systems_of_equations_ex6.C File Reference

Go to the source code of this file.

Classes
class	PetscSolverConfiguration
class	LinearElasticity

Functions
int	main (int argc, char **argv)

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 414 of file systems_of_equations_ex6.C.

References libMesh::DofMap::add_dirichlet_boundary(), libMesh::BoundaryInfo::add_edge(), libMesh::BoundaryInfo::add_node(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::System::attach_assemble_object(), libMesh::MeshTools::Generation::build_cube(), libMesh::command_line_next(), LinearElasticity::compute_stresses(), libMesh::default_solver_package(), dim, libMesh::FIRST, libMesh::MeshBase::get_boundary_info(), libMesh::System::get_dof_map(), libMesh::LinearImplicitSystem::get_linear_solver(), libMesh::BoundaryInfo::has_boundary_id(), libMesh::HEX8, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, libMesh::L2_LAGRANGE, libMesh::LAGRANGE, libMesh::libmesh_assert(), libMesh::LOCAL_VARIABLE_ORDER, mesh, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::LinearImplicitSystem::solve(), and libMesh::ExodusII_IO::write_discontinuous_exodusII().

{
  // Initialize libMesh and any dependent libraries
  LibMeshInit init (argc, argv);

  // This example requires a linear solver package.
  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
                           "--enable-petsc, --enable-trilinos, or --enable-eigen");

  // Initialize the cantilever mesh
  const unsigned int dim = 3;

  // Make sure libMesh was compiled for 3D
  libmesh_example_requires(dim == LIBMESH_DIM, "3D support");

  // We use Dirichlet boundary conditions here
#ifndef LIBMESH_ENABLE_DIRICHLET
  libmesh_example_requires(false, "--enable-dirichlet");
#endif

  // Get the mesh size from the command line.
  const int nx = libMesh::command_line_next("-nx", 32),
            ny = libMesh::command_line_next("-ny", 8),
            nz = libMesh::command_line_next("-nz", 4);

  // Create a 3D mesh distributed across the default MPI communicator.
  Mesh mesh(init.comm(), dim);
  MeshTools::Generation::build_cube (mesh,
                                     nx,
                                     ny,
                                     nz,
                                     0., 1.*x_scaling,
                                     0., 0.3,
                                     0., 0.1,
                                     HEX8);

  // Print information about the mesh to the screen.
  mesh.print_info();

  // Let's add some node and edge boundary conditions
  // Each processor should know about each boundary condition it can
  // see, so we loop over all elements, not just local elements.
  for (const auto & elem : mesh.element_ptr_range())
    {
      unsigned int
        side_max_x = 0, side_min_y = 0,
        side_max_y = 0, side_max_z = 0;

      bool
        found_side_max_x = false, found_side_max_y = false,
        found_side_min_y = false, found_side_max_z = false;

      for (auto side : elem->side_index_range())
        {
          if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_X))
            {
              side_max_x = side;
              found_side_max_x = true;
            }

          if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MIN_Y))
            {
              side_min_y = side;
              found_side_min_y = true;
            }

          if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Y))
            {
              side_max_y = side;
              found_side_max_y = true;
            }

          if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Z))
            {
              side_max_z = side;
              found_side_max_z = true;
            }
        }

      // If elem has sides on boundaries
      // BOUNDARY_ID_MAX_X, BOUNDARY_ID_MAX_Y, BOUNDARY_ID_MAX_Z
      // then let's set a node boundary condition
      if (found_side_max_x && found_side_max_y && found_side_max_z)
        for (auto n : elem->node_index_range())
          if (elem->is_node_on_side(n, side_max_x) &&
              elem->is_node_on_side(n, side_max_y) &&
              elem->is_node_on_side(n, side_max_z))
            mesh.get_boundary_info().add_node(elem->node_ptr(n), NODE_BOUNDARY_ID);


      // If elem has sides on boundaries
      // BOUNDARY_ID_MAX_X and BOUNDARY_ID_MIN_Y
      // then let's set an edge boundary condition
      if (found_side_max_x && found_side_min_y)
        for (auto e : elem->edge_index_range())
          if (elem->is_edge_on_side(e, side_max_x) &&
              elem->is_edge_on_side(e, side_min_y))
            mesh.get_boundary_info().add_edge(elem, e, EDGE_BOUNDARY_ID);
    }

  // Create an equation systems object.
  EquationSystems equation_systems (mesh);

  // Declare the system and its variables.
  // Create a system named "Elasticity"
  LinearImplicitSystem & system =
    equation_systems.add_system<LinearImplicitSystem> ("Elasticity");

#ifdef LIBMESH_HAVE_PETSC
  // Attach a SolverConfiguration object to system.linear_solver
  PetscLinearSolver<Number> * petsc_linear_solver =
    cast_ptr<PetscLinearSolver<Number>*>(system.get_linear_solver());
  libmesh_assert(petsc_linear_solver);
  PetscSolverConfiguration petsc_solver_config(*petsc_linear_solver);
  petsc_linear_solver->set_solver_configuration(petsc_solver_config);
#endif

  LinearElasticity le(equation_systems);
  system.attach_assemble_object(le);

#ifdef LIBMESH_ENABLE_DIRICHLET
  // Add three displacement variables, u and v, to the system
  unsigned int u_var = system.add_variable("u", FIRST, LAGRANGE);
  unsigned int v_var = system.add_variable("v", FIRST, LAGRANGE);
  unsigned int w_var = system.add_variable("w", FIRST, LAGRANGE);

  // Create a ZeroFunction to initialize dirichlet_bc
  ZeroFunction<> zf;

  // Most DirichletBoundary users will want to supply a "locally
  // indexed" functor
  DirichletBoundary dirichlet_bc({BOUNDARY_ID_MIN_X, NODE_BOUNDARY_ID,
                                 EDGE_BOUNDARY_ID},
                                 {u_var, v_var, w_var}, zf,
                                 LOCAL_VARIABLE_ORDER);

  // We must add the Dirichlet boundary condition _before_
  // we call equation_systems.init()
  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
#endif // LIBMESH_ENABLE_DIRICHLET

  // Also, initialize an ExplicitSystem to store stresses
  ExplicitSystem & stress_system =
    equation_systems.add_system<ExplicitSystem> ("StressSystem");

  stress_system.add_variable("sigma_00", FIRST, L2_LAGRANGE);
  stress_system.add_variable("sigma_01", FIRST, L2_LAGRANGE);
  stress_system.add_variable("sigma_02", FIRST, L2_LAGRANGE);
  stress_system.add_variable("sigma_11", FIRST, L2_LAGRANGE);
  stress_system.add_variable("sigma_12", FIRST, L2_LAGRANGE);
  stress_system.add_variable("sigma_22", FIRST, L2_LAGRANGE);
  stress_system.add_variable("vonMises", FIRST, L2_LAGRANGE);

  // Initialize the data structures for the equation system.
  equation_systems.init();

  // Print information about the system to the screen.
  equation_systems.print_info();

  // Solve the system
  system.solve();

  // Post-process the solution to compute the stresses
  le.compute_stresses();

  // Plot the solution
#ifdef LIBMESH_HAVE_EXODUS_API

  // Use single precision in this case (reduces the size of the exodus file)
  ExodusII_IO exo_io(mesh, /*single_precision=*/true);
  exo_io.write_discontinuous_exodusII("displacement_and_stress.exo", equation_systems);

#endif // #ifdef LIBMESH_HAVE_EXODUS_API

  // All done.
  return 0;
}