1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 #include "libmesh/enum_solver_package.h"
59
60 // Local includes
61 #include "linear_elasticity_with_contact.h"
62
63 using namespace libMesh;
64
main(int argc,char ** argv)65 int main (int argc, char ** argv)
66 {
67 LibMeshInit init (argc, argv);
68
69 // This example uses an ExodusII input file
70 #ifndef LIBMESH_HAVE_EXODUS_API
71 libmesh_example_requires(false, "--enable-exodus");
72 #endif
73
74 // We use a 3D domain.
75 libmesh_example_requires(LIBMESH_DIM > 2, "--disable-1D-only --disable-2D-only");
76
77 // We use Dirichlet boundary conditions here
78 #ifndef LIBMESH_ENABLE_DIRICHLET
79 libmesh_example_requires(false, "--enable-dirichlet");
80 #endif
81
82 // This example requires the PETSc nonlinear solvers
83 libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
84
85 GetPot infile("systems_of_equations_ex8.in");
86 const std::string approx_order = infile("approx_order", "FIRST");
87 const std::string fe_family = infile("fe_family", "LAGRANGE");
88
89 const Real young_modulus = infile("Young_modulus", 1.0);
90 const Real poisson_ratio = infile("poisson_ratio", 0.3);
91
92 const Real nonlinear_abs_tol = infile("nonlinear_abs_tol", 1.e-8);
93 const Real nonlinear_rel_tol = infile("nonlinear_rel_tol", 1.e-8);
94 const unsigned int nonlinear_max_its = infile("nonlinear_max_its", 50);
95 const Real contact_penalty = infile("contact_penalty", 1.e2);
96 const Real gap_function_tol = infile("gap_function_tol", 1.e-8);
97
98 // This example code has not been written to cope with a distributed mesh
99 ReplicatedMesh mesh(init.comm());
100 mesh.read("systems_of_equations_ex8.exo");
101
102 mesh.print_info();
103
104 EquationSystems equation_systems (mesh);
105
106 NonlinearImplicitSystem & system =
107 equation_systems.add_system<NonlinearImplicitSystem> ("NonlinearElasticity");
108
109 LinearElasticityWithContact le(system, contact_penalty);
110
111 unsigned int u_var =
112 system.add_variable("u",
113 Utility::string_to_enum<Order> (approx_order),
114 Utility::string_to_enum<FEFamily>(fe_family));
115
116 unsigned int v_var =
117 system.add_variable("v",
118 Utility::string_to_enum<Order> (approx_order),
119 Utility::string_to_enum<FEFamily>(fe_family));
120
121 unsigned int w_var =
122 system.add_variable("w",
123 Utility::string_to_enum<Order> (approx_order),
124 Utility::string_to_enum<FEFamily>(fe_family));
125
126 // Also, initialize an ExplicitSystem to store stresses
127 ExplicitSystem & stress_system =
128 equation_systems.add_system<ExplicitSystem> ("StressSystem");
129
130 stress_system.add_variable("sigma_00", CONSTANT, MONOMIAL);
131 stress_system.add_variable("sigma_01", CONSTANT, MONOMIAL);
132 stress_system.add_variable("sigma_02", CONSTANT, MONOMIAL);
133 stress_system.add_variable("sigma_11", CONSTANT, MONOMIAL);
134 stress_system.add_variable("sigma_12", CONSTANT, MONOMIAL);
135 stress_system.add_variable("sigma_22", CONSTANT, MONOMIAL);
136 stress_system.add_variable("vonMises", CONSTANT, MONOMIAL);
137
138 equation_systems.parameters.set<Real> ("nonlinear solver absolute residual tolerance") = nonlinear_abs_tol;
139 equation_systems.parameters.set<Real> ("nonlinear solver relative residual tolerance") = nonlinear_rel_tol;
140 equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum iterations") = nonlinear_max_its;
141
142 system.nonlinear_solver->residual_and_jacobian_object = ≤
143
144 equation_systems.parameters.set<Real>("young_modulus") = young_modulus;
145 equation_systems.parameters.set<Real>("poisson_ratio") = poisson_ratio;
146
147 #ifdef LIBMESH_ENABLE_DIRICHLET
148 // Attach Dirichlet boundary conditions
149 {
150 std::set<boundary_id_type> clamped_boundaries;
151 clamped_boundaries.insert(MIN_Z_BOUNDARY);
152
153 std::vector<unsigned int> uvw;
154 uvw.push_back(u_var);
155 uvw.push_back(v_var);
156 uvw.push_back(w_var);
157
158 ZeroFunction<Number> zero;
159
160 // Most DirichletBoundary users will want to supply a "locally
161 // indexed" functor
162 system.get_dof_map().add_dirichlet_boundary
163 (DirichletBoundary (clamped_boundaries, uvw, zero,
164 LOCAL_VARIABLE_ORDER));
165 }
166 {
167 std::set<boundary_id_type> clamped_boundaries;
168 clamped_boundaries.insert(MAX_Z_BOUNDARY);
169
170 std::vector<unsigned int> uv;
171 uv.push_back(u_var);
172 uv.push_back(v_var);
173
174 ZeroFunction<Number> zero;
175
176 system.get_dof_map().add_dirichlet_boundary
177 (DirichletBoundary (clamped_boundaries, uv, zero,
178 LOCAL_VARIABLE_ORDER));
179 }
180 {
181 std::set<boundary_id_type> clamped_boundaries;
182 clamped_boundaries.insert(MAX_Z_BOUNDARY);
183
184 std::vector<unsigned int> w;
185 w.push_back(w_var);
186
187 ConstFunction<Number> neg_one(-1.);
188
189 system.get_dof_map().add_dirichlet_boundary
190 (DirichletBoundary (clamped_boundaries, w, neg_one,
191 LOCAL_VARIABLE_ORDER));
192 }
193 #else
194 libmesh_ignore(u_var, v_var, w_var);
195 #endif // LIBMESH_ENABLE_DIRICHLET
196
197 le.initialize_contact_load_paths();
198
199 libMesh::out << "Mesh before adding edge connectors" << std::endl;
200 mesh.print_info();
201 le.add_contact_edge_elements();
202
203 libMesh::out << "Mesh after adding edge connectors" << std::endl;
204 mesh.print_info();
205
206 equation_systems.init();
207 equation_systems.print_info();
208
209 libMesh::out << "Contact penalty: " << contact_penalty << std::endl << std::endl;
210
211 Real current_max_gap_function = std::numeric_limits<Real>::max();
212
213 const unsigned int max_outer_iterations = 10;
214 for (unsigned int outer_iteration = 0;
215 outer_iteration != max_outer_iterations; ++outer_iteration)
216 {
217 if (current_max_gap_function <= gap_function_tol)
218 {
219 libMesh::out << "Outer loop converged" << std::endl;
220 break;
221 }
222
223 libMesh::out << "Starting outer iteration " << outer_iteration << std::endl;
224
225 // Perform inner iteration (i.e. Newton's method loop)
226 system.solve();
227 system.nonlinear_solver->print_converged_reason();
228
229 // Perform augmented Lagrangian update
230 le.update_lambdas();
231
232 std::pair<Real, Real> least_and_max_gap_function = le.get_least_and_max_gap_function();
233 Real least_gap_fn = least_and_max_gap_function.first;
234 Real max_gap_fn = least_and_max_gap_function.second;
235
236 libMesh::out << "Finished outer iteration, least gap function: "
237 << least_gap_fn
238 << ", max gap function: "
239 << max_gap_fn
240 << std::endl
241 << std::endl;
242
243 current_max_gap_function = std::max(std::abs(least_gap_fn), std::abs(max_gap_fn));
244 }
245
246 libMesh::out << "Computing stresses..." << std::endl;
247
248 le.compute_stresses();
249
250 std::stringstream filename;
251 filename << "solution.exo";
252 ExodusII_IO (mesh).write_equation_systems(filename.str(),
253 equation_systems);
254
255 return 0;
256 }
257