libMesh
Functions
systems_of_equations_ex8.C File Reference

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 64 of file systems_of_equations_ex8.C.

References std::abs(), libMesh::DofMap::add_dirichlet_boundary(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::LibMeshInit::comm(), libMesh::CONSTANT, libMesh::default_solver_package(), libMesh::System::get_dof_map(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::LOCAL_VARIABLE_ORDER, std::max(), 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::MeshOutput< MT >::write_equation_systems(), and libMesh::zero.

65 {
66  LibMeshInit init (argc, argv);
67 
68  // This example uses an ExodusII input file
69 #ifndef LIBMESH_HAVE_EXODUS_API
70  libmesh_example_requires(false, "--enable-exodus");
71 #endif
72 
73  // We use a 3D domain.
74  libmesh_example_requires(LIBMESH_DIM > 2, "--disable-1D-only --disable-2D-only");
75 
76  // This example requires the PETSc nonlinear solvers
77  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
78 
79  GetPot infile("systems_of_equations_ex8.in");
80  const std::string approx_order = infile("approx_order", "FIRST");
81  const std::string fe_family = infile("fe_family", "LAGRANGE");
82 
83  const Real young_modulus = infile("Young_modulus", 1.0);
84  const Real poisson_ratio = infile("poisson_ratio", 0.3);
85 
86  const Real nonlinear_abs_tol = infile("nonlinear_abs_tol", 1.e-8);
87  const Real nonlinear_rel_tol = infile("nonlinear_rel_tol", 1.e-8);
88  const unsigned int nonlinear_max_its = infile("nonlinear_max_its", 50);
89  const Real contact_penalty = infile("contact_penalty", 1.e2);
90  const Real gap_function_tol = infile("gap_function_tol", 1.e-8);
91 
92  // This example code has not been written to cope with a distributed mesh
93  ReplicatedMesh mesh(init.comm());
94  mesh.read("systems_of_equations_ex8.exo");
95 
96  mesh.print_info();
97 
98  EquationSystems equation_systems (mesh);
99 
100  NonlinearImplicitSystem & system =
101  equation_systems.add_system<NonlinearImplicitSystem> ("NonlinearElasticity");
102 
103  LinearElasticityWithContact le(system, contact_penalty);
104 
105  unsigned int u_var =
106  system.add_variable("u",
107  Utility::string_to_enum<Order> (approx_order),
108  Utility::string_to_enum<FEFamily>(fe_family));
109 
110  unsigned int v_var =
111  system.add_variable("v",
112  Utility::string_to_enum<Order> (approx_order),
113  Utility::string_to_enum<FEFamily>(fe_family));
114 
115  unsigned int w_var =
116  system.add_variable("w",
117  Utility::string_to_enum<Order> (approx_order),
118  Utility::string_to_enum<FEFamily>(fe_family));
119 
120  // Also, initialize an ExplicitSystem to store stresses
121  ExplicitSystem & stress_system =
122  equation_systems.add_system<ExplicitSystem> ("StressSystem");
123 
124  stress_system.add_variable("sigma_00", CONSTANT, MONOMIAL);
125  stress_system.add_variable("sigma_01", CONSTANT, MONOMIAL);
126  stress_system.add_variable("sigma_02", CONSTANT, MONOMIAL);
127  stress_system.add_variable("sigma_11", CONSTANT, MONOMIAL);
128  stress_system.add_variable("sigma_12", CONSTANT, MONOMIAL);
129  stress_system.add_variable("sigma_22", CONSTANT, MONOMIAL);
130  stress_system.add_variable("vonMises", CONSTANT, MONOMIAL);
131 
132  equation_systems.parameters.set<Real> ("nonlinear solver absolute residual tolerance") = nonlinear_abs_tol;
133  equation_systems.parameters.set<Real> ("nonlinear solver relative residual tolerance") = nonlinear_rel_tol;
134  equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum iterations") = nonlinear_max_its;
135 
136  system.nonlinear_solver->residual_and_jacobian_object = &le;
137 
138  equation_systems.parameters.set<Real>("young_modulus") = young_modulus;
139  equation_systems.parameters.set<Real>("poisson_ratio") = poisson_ratio;
140 
141  // Attach Dirichlet boundary conditions
142  {
143  std::set<boundary_id_type> clamped_boundaries;
144  clamped_boundaries.insert(MIN_Z_BOUNDARY);
145 
146  std::vector<unsigned int> uvw;
147  uvw.push_back(u_var);
148  uvw.push_back(v_var);
149  uvw.push_back(w_var);
150 
152 
153  // Most DirichletBoundary users will want to supply a "locally
154  // indexed" functor
156  (DirichletBoundary (clamped_boundaries, uvw, zero,
158  }
159  {
160  std::set<boundary_id_type> clamped_boundaries;
161  clamped_boundaries.insert(MAX_Z_BOUNDARY);
162 
163  std::vector<unsigned int> uv;
164  uv.push_back(u_var);
165  uv.push_back(v_var);
166 
168 
170  (DirichletBoundary (clamped_boundaries, uv, zero,
172  }
173  {
174  std::set<boundary_id_type> clamped_boundaries;
175  clamped_boundaries.insert(MAX_Z_BOUNDARY);
176 
177  std::vector<unsigned int> w;
178  w.push_back(w_var);
179 
180  ConstFunction<Number> neg_one(-1.);
181 
183  (DirichletBoundary (clamped_boundaries, w, neg_one,
185  }
186 
187  le.initialize_contact_load_paths();
188 
189  libMesh::out << "Mesh before adding edge connectors" << std::endl;
190  mesh.print_info();
191  le.add_contact_edge_elements();
192 
193  libMesh::out << "Mesh after adding edge connectors" << std::endl;
194  mesh.print_info();
195 
196  equation_systems.init();
197  equation_systems.print_info();
198 
199  libMesh::out << "Contact penalty: " << contact_penalty << std::endl << std::endl;
200 
201  Real current_max_gap_function = std::numeric_limits<Real>::max();
202 
203  const unsigned int max_outer_iterations = 10;
204  for (unsigned int outer_iteration = 0;
205  outer_iteration != max_outer_iterations; ++outer_iteration)
206  {
207  if (current_max_gap_function <= gap_function_tol)
208  {
209  libMesh::out << "Outer loop converged" << std::endl;
210  break;
211  }
212 
213  libMesh::out << "Starting outer iteration " << outer_iteration << std::endl;
214 
215  // Perform inner iteration (i.e. Newton's method loop)
216  system.solve();
217  system.nonlinear_solver->print_converged_reason();
218 
219  // Perform augmented Lagrangian update
220  le.update_lambdas();
221 
222  std::pair<Real, Real> least_and_max_gap_function = le.get_least_and_max_gap_function();
223  Real least_gap_fn = least_and_max_gap_function.first;
224  Real max_gap_fn = least_and_max_gap_function.second;
225 
226  libMesh::out << "Finished outer iteration, least gap function: "
227  << least_gap_fn
228  << ", max gap function: "
229  << max_gap_fn
230  << std::endl
231  << std::endl;
232 
233  current_max_gap_function = std::max(std::abs(least_gap_fn), std::abs(max_gap_fn));
234  }
235 
236  libMesh::out << "Computing stresses..." << std::endl;
237 
238  le.compute_stresses();
239 
240  std::stringstream filename;
241  filename << "solution.exo";
242  ExodusII_IO (mesh).write_equation_systems(filename.str(),
243  equation_systems);
244 
245  return 0;
246 }
double abs(double a)
This is the EquationSystems class.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
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 encapsulate all functionality required for assembling and solving a linear elastic model w...
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
long double max(long double a, double b)
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
Function that returns a single value that never changes.
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
The ExplicitSystem provides only "right hand side" storage, which should be sufficient for solving mo...