libMesh
Classes | Functions
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

int main ( int  argc,
char **  argv 
)

Definition at line 418 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::LibMeshInit::comm(), LinearElasticity::compute_stresses(), libMesh::default_solver_package(), dim, libMesh::MeshBase::element_ptr_range(), 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(), side, libMesh::LinearImplicitSystem::solve(), and libMesh::ExodusII_IO::write_discontinuous_exodusII().

419 {
420  // Initialize libMesh and any dependent libraries
421  LibMeshInit init (argc, argv);
422 
423  // This example requires a linear solver package.
424  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
425  "--enable-petsc, --enable-trilinos, or --enable-eigen");
426 
427  // Initialize the cantilever mesh
428  const unsigned int dim = 3;
429 
430  // Make sure libMesh was compiled for 3D
431  libmesh_example_requires(dim == LIBMESH_DIM, "3D support");
432 
433  // Create a 3D mesh distributed across the default MPI communicator.
434  Mesh mesh(init.comm(), dim);
436  32,
437  8,
438  4,
439  0., 1.*x_scaling,
440  0., 0.3,
441  0., 0.1,
442  HEX8);
443 
444  // Print information about the mesh to the screen.
445  mesh.print_info();
446 
447  // Let's add some node and edge boundary conditions
448  // Each processor should know about each boundary condition it can
449  // see, so we loop over all elements, not just local elements.
450  for (const auto & elem : mesh.element_ptr_range())
451  {
452  unsigned int
453  side_max_x = 0, side_min_y = 0,
454  side_max_y = 0, side_max_z = 0;
455 
456  bool
457  found_side_max_x = false, found_side_max_y = false,
458  found_side_min_y = false, found_side_max_z = false;
459 
460  for (auto side : elem->side_index_range())
461  {
462  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_X))
463  {
464  side_max_x = side;
465  found_side_max_x = true;
466  }
467 
468  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MIN_Y))
469  {
470  side_min_y = side;
471  found_side_min_y = true;
472  }
473 
474  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Y))
475  {
476  side_max_y = side;
477  found_side_max_y = true;
478  }
479 
480  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Z))
481  {
482  side_max_z = side;
483  found_side_max_z = true;
484  }
485  }
486 
487  // If elem has sides on boundaries
488  // BOUNDARY_ID_MAX_X, BOUNDARY_ID_MAX_Y, BOUNDARY_ID_MAX_Z
489  // then let's set a node boundary condition
490  if (found_side_max_x && found_side_max_y && found_side_max_z)
491  for (auto n : elem->node_index_range())
492  if (elem->is_node_on_side(n, side_max_x) &&
493  elem->is_node_on_side(n, side_max_y) &&
494  elem->is_node_on_side(n, side_max_z))
495  mesh.get_boundary_info().add_node(elem->node_ptr(n), NODE_BOUNDARY_ID);
496 
497 
498  // If elem has sides on boundaries
499  // BOUNDARY_ID_MAX_X and BOUNDARY_ID_MIN_Y
500  // then let's set an edge boundary condition
501  if (found_side_max_x && found_side_min_y)
502  for (auto e : elem->edge_index_range())
503  if (elem->is_edge_on_side(e, side_max_x) &&
504  elem->is_edge_on_side(e, side_min_y))
505  mesh.get_boundary_info().add_edge(elem, e, EDGE_BOUNDARY_ID);
506  }
507 
508  // Create an equation systems object.
509  EquationSystems equation_systems (mesh);
510 
511  // Declare the system and its variables.
512  // Create a system named "Elasticity"
513  LinearImplicitSystem & system =
514  equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
515 
516 #ifdef LIBMESH_HAVE_PETSC
517  // Attach a SolverConfiguration object to system.linear_solver
518  PetscLinearSolver<Number> * petsc_linear_solver =
519  cast_ptr<PetscLinearSolver<Number>*>(system.get_linear_solver());
520  libmesh_assert(petsc_linear_solver);
521  PetscSolverConfiguration petsc_solver_config(*petsc_linear_solver);
522  petsc_linear_solver->set_solver_configuration(petsc_solver_config);
523 #endif
524 
525  // Add three displacement variables, u and v, to the system
526  unsigned int u_var = system.add_variable("u", FIRST, LAGRANGE);
527  unsigned int v_var = system.add_variable("v", FIRST, LAGRANGE);
528  unsigned int w_var = system.add_variable("w", FIRST, LAGRANGE);
529 
530  LinearElasticity le(equation_systems);
531  system.attach_assemble_object(le);
532 
533  std::set<boundary_id_type> boundary_ids;
534  boundary_ids.insert(BOUNDARY_ID_MIN_X);
535  boundary_ids.insert(NODE_BOUNDARY_ID);
536  boundary_ids.insert(EDGE_BOUNDARY_ID);
537 
538  // Create a vector storing the variable numbers which the BC applies to
539  std::vector<unsigned int> variables;
540  variables.push_back(u_var);
541  variables.push_back(v_var);
542  variables.push_back(w_var);
543 
544  // Create a ZeroFunction to initialize dirichlet_bc
545  ZeroFunction<> zf;
546 
547  // Most DirichletBoundary users will want to supply a "locally
548  // indexed" functor
549  DirichletBoundary dirichlet_bc(boundary_ids, variables, zf,
551 
552  // We must add the Dirichlet boundary condition _before_
553  // we call equation_systems.init()
554  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
555 
556  // Also, initialize an ExplicitSystem to store stresses
557  ExplicitSystem & stress_system =
558  equation_systems.add_system<ExplicitSystem> ("StressSystem");
559 
560  stress_system.add_variable("sigma_00", FIRST, L2_LAGRANGE);
561  stress_system.add_variable("sigma_01", FIRST, L2_LAGRANGE);
562  stress_system.add_variable("sigma_02", FIRST, L2_LAGRANGE);
563  stress_system.add_variable("sigma_11", FIRST, L2_LAGRANGE);
564  stress_system.add_variable("sigma_12", FIRST, L2_LAGRANGE);
565  stress_system.add_variable("sigma_22", FIRST, L2_LAGRANGE);
566  stress_system.add_variable("vonMises", FIRST, L2_LAGRANGE);
567 
568  // Initialize the data structures for the equation system.
569  equation_systems.init();
570 
571  // Print information about the system to the screen.
572  equation_systems.print_info();
573 
574  // Solve the system
575  system.solve();
576 
577  // Post-process the solution to compute the stresses
578  le.compute_stresses();
579 
580  // Plot the solution
581 #ifdef LIBMESH_HAVE_EXODUS_API
582 
583  // Use single precision in this case (reduces the size of the exodus file)
584  ExodusII_IO exo_io(mesh, /*single_precision=*/true);
585  exo_io.write_discontinuous_exodusII("displacement_and_stress.exo", equation_systems);
586 
587 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
588 
589  // All done.
590  return 0;
591 }
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
This is the EquationSystems class.
ConstFunction that simply returns 0.
Definition: zero_function.h:35
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
This class provides a specific system class.
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
unsigned short int side
Definition: xdr_io.C:49
MeshBase & mesh
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
libmesh_assert(j)
SolverPackage default_solver_package()
Definition: libmesh.C:995
void add_node(const Node *node, const boundary_id_type id)
Add Node node with boundary id id to the boundary information data structures.
virtual LinearSolver< Number > * get_linear_solver() const libmesh_override
const DofMap & get_dof_map() const
Definition: system.h:2030
virtual SimpleRange< element_iterator > element_ptr_range()=0
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
void attach_assemble_object(Assembly &assemble)
Register a user object to use in assembling the system matrix and RHS.
Definition: system.C:1814
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
The ExplicitSystem provides only "right hand side" storage, which should be sufficient for solving mo...
void add_edge(const dof_id_type elem, const unsigned short int edge, const boundary_id_type id)
Add edge edge of element number elem with boundary id id to the boundary information data structure...