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

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 509 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(), 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::libmesh_ignore(), 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::System::prefer_hash_table_matrix_assembly(), libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::Real, libMesh::Parameters::set(), libMesh::NonlinearImplicitSystem::solve(), libMesh::ExodusII_IO::write_equation_systems(), and libMesh::zero.

{
  LibMeshInit init (argc, argv);

  // This example requires the PETSc nonlinear solvers
  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");

  // We use a 3D domain.
  libmesh_example_requires(LIBMESH_DIM > 2, "--disable-1D-only --disable-2D-only");

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

  GetPot infile("systems_of_equations_ex7.in");

  infile.parse_command_line(argc,argv);

  const Real x_length = infile("x_length", 0.);
  const Real y_length = infile("y_length", 0.);
  const Real z_length = infile("z_length", 0.);
  const int n_elem_x = infile("n_elem_x", 0);
  const int n_elem_y = infile("n_elem_y", 0);
  const int n_elem_z = infile("n_elem_z", 0);
  const std::string approx_order = infile("approx_order", "FIRST");
  const std::string fe_family = infile("fe_family", "LAGRANGE");

  const Real young_modulus = infile("Young_modulus", 1.0);
  const Real poisson_ratio = infile("poisson_ratio", 0.3);
  const Real forcing_magnitude = infile("forcing_magnitude", 0.001);

  const Real nonlinear_abs_tol = infile("nonlinear_abs_tol", 1.e-8);
  const Real nonlinear_rel_tol = infile("nonlinear_rel_tol", 1.e-8);
  const unsigned int nonlinear_max_its = infile("nonlinear_max_its", 50);

  const unsigned int n_solves = infile("n_solves", 10);
  const Real force_scaling = infile("force_scaling", 5.0);

  Mesh mesh(init.comm());

  MeshTools::Generation::build_cube(mesh,
                                    n_elem_x,
                                    n_elem_y,
                                    n_elem_z,
                                    0., x_length,
                                    0., y_length,
                                    0., z_length,
                                    HEX27);

  mesh.print_info();

  EquationSystems equation_systems (mesh);
  LargeDeformationElasticity lde(equation_systems);

  NonlinearImplicitSystem & system =
    equation_systems.add_system<NonlinearImplicitSystem> ("NonlinearElasticity");
  system.prefer_hash_table_matrix_assembly(true);

  unsigned int u_var =
    system.add_variable("u",
                        Utility::string_to_enum<Order>   (approx_order),
                        Utility::string_to_enum<FEFamily>(fe_family));

  unsigned int v_var =
    system.add_variable("v",
                        Utility::string_to_enum<Order>   (approx_order),
                        Utility::string_to_enum<FEFamily>(fe_family));

  unsigned int w_var =
    system.add_variable("w",
                        Utility::string_to_enum<Order>   (approx_order),
                        Utility::string_to_enum<FEFamily>(fe_family));

  // Also, initialize an ExplicitSystem to store stresses
  ExplicitSystem & stress_system =
    equation_systems.add_system<ExplicitSystem> ("StressSystem");
  stress_system.add_variable("sigma_00", CONSTANT, MONOMIAL);
  stress_system.add_variable("sigma_01", CONSTANT, MONOMIAL);
  stress_system.add_variable("sigma_02", CONSTANT, MONOMIAL);
  stress_system.add_variable("sigma_11", CONSTANT, MONOMIAL);
  stress_system.add_variable("sigma_12", CONSTANT, MONOMIAL);
  stress_system.add_variable("sigma_22", CONSTANT, MONOMIAL);

  equation_systems.parameters.set<Real>         ("nonlinear solver absolute residual tolerance") = nonlinear_abs_tol;
  equation_systems.parameters.set<Real>         ("nonlinear solver relative residual tolerance") = nonlinear_rel_tol;
  equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum iterations")          = nonlinear_max_its;

  system.nonlinear_solver->residual_object = &lde;
  system.nonlinear_solver->jacobian_object = &lde;

  equation_systems.parameters.set<Real>("young_modulus") = young_modulus;
  equation_systems.parameters.set<Real>("poisson_ratio") = poisson_ratio;
  equation_systems.parameters.set<Real>("forcing_magnitude") = forcing_magnitude;

#ifdef LIBMESH_ENABLE_DIRICHLET
  // Attach Dirichlet (Clamped) boundary conditions
  ZeroFunction<Number> zero;

  // Most DirichletBoundary users will want to supply a "locally
  // indexed" functor
  system.get_dof_map().add_dirichlet_boundary
    (DirichletBoundary ({BOUNDARY_ID_MIN_X}, {u_var, v_var, w_var},
                        zero, LOCAL_VARIABLE_ORDER));
#endif // LIBMESH_ENABLE_DIRICHLET
  libmesh_ignore(u_var, v_var, w_var);

  equation_systems.init();
  equation_systems.print_info();

  // Provide a loop here so that we can do a sequence of solves
  // where solve n gives a good starting guess for solve n+1.
  // This "continuation" approach is helpful for solving for
  // large values of "forcing_magnitude".
  // Set n_solves and force_scaling in nonlinear_elasticity.in.
  for (unsigned int count=0; count<n_solves; count++)
    {
      Real previous_forcing_magnitude = equation_systems.parameters.get<Real>("forcing_magnitude");
      equation_systems.parameters.set<Real>("forcing_magnitude") = previous_forcing_magnitude*force_scaling;

      libMesh::out << "Performing solve "
                   << count
                   << ", forcing_magnitude: "
                   << equation_systems.parameters.get<Real>("forcing_magnitude")
                   << std::endl;

      system.solve();

      libMesh::out << "System solved at nonlinear iteration "
                   << system.n_nonlinear_iterations()
                   << " , final nonlinear residual norm: "
                   << system.final_nonlinear_residual()
                   << std::endl
                   << std::endl;

      libMesh::out << "Computing stresses..." << std::endl;

      lde.compute_stresses();

#ifdef LIBMESH_HAVE_EXODUS_API
      std::stringstream filename;
      filename << "solution_" << count << ".exo";
      ExodusII_IO (mesh).write_equation_systems(filename.str(), equation_systems);
#endif
    }

  return 0;
}