libMesh
systems_of_equations_ex8.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // <h1> Systems Example 8 - "Small-sliding" contact. </h1>
21 // \author David Knezevic
22 // \date 2015
23 //
24 // In this example, we consider a linear elastic model with contact. We restrict ourselves
25 // to considering "small sliding", which means that we set the contact "load transfer" between
26 // surfaces up front, based on the undeformed mesh, and retain that load transfer throughout
27 // the solve. Even though we consider linear elasticity here, this is a nonlinear problem due
28 // to the contact condition.
29 //
30 // The contact condition is imposed using the augmented Lagrangian method, e.g. see
31 // Simo & Laursen (1992). For the sake of simplicity, in this example we assume that contact
32 // nodes are perfectly aligned (this assumption can be eliminated relatively easily).
33 //
34 // The mesh in this example consists of two disconnected cylinders. We add edge elements into
35 // the mesh in order to ensure correct parallel communication of data on contact surfaces,
36 // and also so that we do not have to manually augment the sparsity pattern.
37 
38 // We impose a displacement boundary condition of -1 in the z-direction on the top cylinder
39 // in order to impose contact.
40 
41 // C++ include files that we need
42 #include <iostream>
43 
44 // libMesh includes
45 #include "libmesh/libmesh.h"
46 #include "libmesh/replicated_mesh.h"
47 #include "libmesh/exodusII_io.h"
48 #include "libmesh/equation_systems.h"
49 #include "libmesh/dof_map.h"
50 #include "libmesh/numeric_vector.h"
51 #include "libmesh/getpot.h"
52 #include "libmesh/dirichlet_boundaries.h"
53 #include "libmesh/string_to_enum.h"
54 #include "libmesh/zero_function.h"
55 #include "libmesh/nonlinear_solver.h"
56 #include "libmesh/nonlinear_implicit_system.h"
57 #include "libmesh/petsc_macro.h"
58 
59 // Local includes
61 
62 using namespace libMesh;
63 
64 int main (int argc, char ** argv)
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
int main(int argc, char **argv)
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
The libMesh namespace provides an interface to certain functionality in the library.
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 print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
This class provides a specific system class.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
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
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
Parameters parameters
Data structure holding arbitrary parameters.
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 init()
Initialize all the systems.
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...