libMesh
Classes | Functions
systems_of_equations_ex7.C File Reference

Go to the source code of this file.

Classes

class  LargeDeformationElasticity
 

Functions

int main (int argc, char **argv)
 

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 523 of file systems_of_equations_ex7.C.

References libMesh::DofMap::add_dirichlet_boundary(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::MeshTools::Generation::build_cube(), libMesh::LibMeshInit::comm(), LargeDeformationElasticity::compute_stresses(), libMesh::CONSTANT, libMesh::default_solver_package(), libMesh::NonlinearImplicitSystem::final_nonlinear_residual(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::HEX27, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::LOCAL_VARIABLE_ORDER, mesh, libMesh::MONOMIAL, libMesh::NonlinearImplicitSystem::n_nonlinear_iterations(), libMesh::NonlinearImplicitSystem::nonlinear_solver, libMesh::out, libMesh::EquationSystems::parameters, libMesh::PETSC_SOLVERS, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::Real, libMesh::Parameters::set(), libMesh::NonlinearImplicitSystem::solve(), libMesh::MeshOutput< MT >::write_equation_systems(), and libMesh::zero.

524 {
525  LibMeshInit init (argc, argv);
526 
527  // This example requires the PETSc nonlinear solvers
528  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
529 
530  // We use a 3D domain.
531  libmesh_example_requires(LIBMESH_DIM > 2, "--disable-1D-only --disable-2D-only");
532 
533  GetPot infile("systems_of_equations_ex7.in");
534  const Real x_length = infile("x_length", 0.);
535  const Real y_length = infile("y_length", 0.);
536  const Real z_length = infile("z_length", 0.);
537  const Real n_elem_x = infile("n_elem_x", 0);
538  const Real n_elem_y = infile("n_elem_y", 0);
539  const Real n_elem_z = infile("n_elem_z", 0);
540  const std::string approx_order = infile("approx_order", "FIRST");
541  const std::string fe_family = infile("fe_family", "LAGRANGE");
542 
543  const Real young_modulus = infile("Young_modulus", 1.0);
544  const Real poisson_ratio = infile("poisson_ratio", 0.3);
545  const Real forcing_magnitude = infile("forcing_magnitude", 0.001);
546 
547  const Real nonlinear_abs_tol = infile("nonlinear_abs_tol", 1.e-8);
548  const Real nonlinear_rel_tol = infile("nonlinear_rel_tol", 1.e-8);
549  const unsigned int nonlinear_max_its = infile("nonlinear_max_its", 50);
550 
551  const unsigned int n_solves = infile("n_solves", 10);
552  const Real force_scaling = infile("force_scaling", 5.0);
553 
554  Mesh mesh(init.comm());
555 
557  n_elem_x,
558  n_elem_y,
559  n_elem_z,
560  0., x_length,
561  0., y_length,
562  0., z_length,
563  HEX27);
564 
565  mesh.print_info();
566 
567  EquationSystems equation_systems (mesh);
568  LargeDeformationElasticity lde(equation_systems);
569 
570  NonlinearImplicitSystem & system =
571  equation_systems.add_system<NonlinearImplicitSystem> ("NonlinearElasticity");
572 
573  unsigned int u_var =
574  system.add_variable("u",
575  Utility::string_to_enum<Order> (approx_order),
576  Utility::string_to_enum<FEFamily>(fe_family));
577 
578  unsigned int v_var =
579  system.add_variable("v",
580  Utility::string_to_enum<Order> (approx_order),
581  Utility::string_to_enum<FEFamily>(fe_family));
582 
583  unsigned int w_var =
584  system.add_variable("w",
585  Utility::string_to_enum<Order> (approx_order),
586  Utility::string_to_enum<FEFamily>(fe_family));
587 
588  // Also, initialize an ExplicitSystem to store stresses
589  ExplicitSystem & stress_system =
590  equation_systems.add_system<ExplicitSystem> ("StressSystem");
591  stress_system.add_variable("sigma_00", CONSTANT, MONOMIAL);
592  stress_system.add_variable("sigma_01", CONSTANT, MONOMIAL);
593  stress_system.add_variable("sigma_02", CONSTANT, MONOMIAL);
594  stress_system.add_variable("sigma_11", CONSTANT, MONOMIAL);
595  stress_system.add_variable("sigma_12", CONSTANT, MONOMIAL);
596  stress_system.add_variable("sigma_22", CONSTANT, MONOMIAL);
597 
598  equation_systems.parameters.set<Real> ("nonlinear solver absolute residual tolerance") = nonlinear_abs_tol;
599  equation_systems.parameters.set<Real> ("nonlinear solver relative residual tolerance") = nonlinear_rel_tol;
600  equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum iterations") = nonlinear_max_its;
601 
602  system.nonlinear_solver->residual_object = &lde;
603  system.nonlinear_solver->jacobian_object = &lde;
604 
605  equation_systems.parameters.set<Real>("young_modulus") = young_modulus;
606  equation_systems.parameters.set<Real>("poisson_ratio") = poisson_ratio;
607  equation_systems.parameters.set<Real>("forcing_magnitude") = forcing_magnitude;
608 
609  // Attach Dirichlet boundary conditions
610  std::set<boundary_id_type> clamped_boundaries;
611  clamped_boundaries.insert(BOUNDARY_ID_MIN_X);
612 
613  std::vector<unsigned int> uvw;
614  uvw.push_back(u_var);
615  uvw.push_back(v_var);
616  uvw.push_back(w_var);
617 
619 
620  // Most DirichletBoundary users will want to supply a "locally
621  // indexed" functor
623  (DirichletBoundary (clamped_boundaries, uvw, zero,
625 
626  equation_systems.init();
627  equation_systems.print_info();
628 
629  // Provide a loop here so that we can do a sequence of solves
630  // where solve n gives a good starting guess for solve n+1.
631  // This "continuation" approach is helpful for solving for
632  // large values of "forcing_magnitude".
633  // Set n_solves and force_scaling in nonlinear_elasticity.in.
634  for (unsigned int count=0; count<n_solves; count++)
635  {
636  Real previous_forcing_magnitude = equation_systems.parameters.get<Real>("forcing_magnitude");
637  equation_systems.parameters.set<Real>("forcing_magnitude") = previous_forcing_magnitude*force_scaling;
638 
639  libMesh::out << "Performing solve "
640  << count
641  << ", forcing_magnitude: "
642  << equation_systems.parameters.get<Real>("forcing_magnitude")
643  << std::endl;
644 
645  system.solve();
646 
647  libMesh::out << "System solved at nonlinear iteration "
648  << system.n_nonlinear_iterations()
649  << " , final nonlinear residual norm: "
650  << system.final_nonlinear_residual()
651  << std::endl
652  << std::endl;
653 
654  libMesh::out << "Computing stresses..." << std::endl;
655 
656  lde.compute_stresses();
657 
658 #ifdef LIBMESH_HAVE_EXODUS_API
659  std::stringstream filename;
660  filename << "solution_" << count << ".exo";
661  ExodusII_IO (mesh).write_equation_systems(filename.str(), equation_systems);
662 #endif
663  }
664 
665  return 0;
666 }
This is the EquationSystems class.
ConstFunction that simply returns 0.
Definition: zero_function.h:35
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
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
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
const Number zero
.
Definition: libmesh.h:178
SolverPackage default_solver_package()
Definition: libmesh.C:995
UniquePtr< NonlinearSolver< Number > > nonlinear_solver
The NonlinearSolver defines the default interface used to solve the nonlinear_implicit system...
const DofMap & get_dof_map() const
Definition: system.h:2030
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
This class provides a specific system class.
virtual void solve() libmesh_override
Assembles & solves the nonlinear system R(x) = 0.
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
streamT * get()
Rather than implement every ostream/ios/ios_base function, we&#39;ll be lazy and make esoteric uses go th...
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
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...