systems_of_equations_ex8.C File Reference

Go to the source code of this file.

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

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 66 of file systems_of_equations_ex8.C.

References std::abs(), libMesh::DofMap::add_dirichlet_boundary(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::CONSTANT, libMesh::default_solver_package(), libMesh::DofMap::enforce_constraints_exactly(), libMesh::System::get_dof_map(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::libmesh_ignore(), libMesh::LOCAL_VARIABLE_ORDER, mesh, libMesh::MONOMIAL, libMesh::NonlinearImplicitSystem::nonlinear_solver, libMesh::out, libMesh::EquationSystems::parameters, libMesh::PETSC_SOLVERS, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::MeshBase::read(), libMesh::Real, libMesh::Parameters::set(), libMesh::NonlinearImplicitSystem::solve(), libMesh::MeshRefinement::uniformly_refine(), libMesh::System::update(), libMesh::MeshOutput< MT >::write_equation_systems(), and libMesh::zero.

{
  LibMeshInit init (argc, argv);

  // This example uses an ExodusII input file
#ifndef LIBMESH_HAVE_EXODUS_API
  libmesh_example_requires(false, "--enable-exodus");
#endif

  // 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

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

  GetPot infile("systems_of_equations_ex8.in");
  infile.parse_command_line(argc,argv);

  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 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 Real contact_penalty = infile("contact_penalty", 1.e2);
  const Real gap_function_tol = infile("gap_function_tol", 1.e-8);

  // This example code has not been written to cope with a distributed mesh
  ReplicatedMesh mesh(init.comm());
  mesh.read("systems_of_equations_ex8.exo");

  const unsigned int n_refinements = infile("n_refinements", 0);
  // Skip adaptive runs on a non-adaptive libMesh build
#ifndef LIBMESH_ENABLE_AMR
  libmesh_example_requires(n_refinements==0, "--enable-amr");
#else
  MeshRefinement mesh_refinement(mesh);
  mesh_refinement.uniformly_refine(n_refinements);
#endif

  mesh.print_info();

  EquationSystems equation_systems (mesh);

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

  LinearElasticityWithContact le(system, contact_penalty);

  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);
  stress_system.add_variable("vonMises", 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_and_jacobian_object = &le;

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

#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 ({MIN_Z_BOUNDARY}, {u_var, v_var, w_var},
                        zero, LOCAL_VARIABLE_ORDER));

  system.get_dof_map().add_dirichlet_boundary
    (DirichletBoundary ({MAX_Z_BOUNDARY}, {u_var, v_var}, zero,
                        LOCAL_VARIABLE_ORDER));

  ConstFunction<Number> neg_one(-1.);

  system.get_dof_map().add_dirichlet_boundary
    (DirichletBoundary ({MAX_Z_BOUNDARY}, {w_var}, neg_one,
                        LOCAL_VARIABLE_ORDER));
#endif // LIBMESH_ENABLE_DIRICHLET
  libmesh_ignore(u_var, v_var, w_var);

  le.initialize_contact_load_paths();

  libMesh::out << "Mesh before adding edge connectors" << std::endl;
  mesh.print_info();
  le.add_contact_edge_elements();

  libMesh::out << "Mesh after adding edge connectors" << std::endl;
  mesh.print_info();

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

  libMesh::out << "Contact penalty: " << contact_penalty  << std::endl << std::endl;

  Real current_max_gap_function = std::numeric_limits<Real>::max();

  const unsigned int max_outer_iterations = 10;
  for (unsigned int outer_iteration = 0;
       outer_iteration != max_outer_iterations; ++outer_iteration)
    {
      if (current_max_gap_function <= gap_function_tol)
        {
          libMesh::out << "Outer loop converged" << std::endl;
          break;
        }

      libMesh::out << "Starting outer iteration " << outer_iteration << std::endl;

      // Perform inner iteration (i.e. Newton's method loop)
      system.solve();
      system.nonlinear_solver->print_converged_reason();

      // Perform augmented Lagrangian update
      le.update_lambdas();

      std::pair<Real, Real> least_and_max_gap_function = le.get_least_and_max_gap_function();
      Real least_gap_fn = least_and_max_gap_function.first;
      Real max_gap_fn = least_and_max_gap_function.second;

      libMesh::out << "Finished outer iteration, least gap function: "
                   << least_gap_fn
                   << ", max gap function: "
                   << max_gap_fn
                   << std::endl
                   << std::endl;

      current_max_gap_function = std::max(std::abs(least_gap_fn), std::abs(max_gap_fn));
    }

  // Enforcing constraints imposes the non-zero Dirichlet constraints on the solution.
  system.get_dof_map().enforce_constraints_exactly(system);
  system.update();

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

  le.compute_stresses();

  std::stringstream filename;
  filename << "solution.exo";
  ExodusII_IO (mesh).write_equation_systems(filename.str(),
                                            equation_systems);

  return 0;
}