MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
contact-patch-test.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11//
12// -----------------------------------------
13// Tribol Miniapp: Mortar contact patch test
14// -----------------------------------------
15//
16// This miniapp depends on Tribol, an open source contact mechanics library
17// available at https://github.com/LLNL/Tribol.
18//
19// This miniapp uses Tribol's mortar method to solve a contact patch test.
20// Tribol has native support for MFEM data structures (ParMesh, ParGridFunction,
21// HypreParMatrix, etc.) which simplifies including contact support in
22// MFEM-based solid mechanics codes. Note the mesh file two_hex.mesh must be in
23// your path for the miniapp to execute correctly. This mesh file contains two
24// cubes occupying [0,1]^3 and [0,1]x[0,1]x[0.99,1.99]. By default, the miniapp
25// will uniformly refine the mesh twice, then split it across MPI ranks. An
26// elasticity bilinear form will be created over the volume mesh and mortar
27// contact constraints will be formed along the z=1 and z=0.99 surfaces of the
28// blocks.
29//
30// Given the elasticity stiffness matrix and the gap constraints and constraint
31// derivatives from Tribol, the miniapp will form and solve a linear system of
32// equations for updated displacements and pressure Lagrange multipliers
33// enforcing the contact constraints. Finally, it will verify force equilibrium
34// and that the gap constraints are satisfied and save output in VisIt format.
35//
36// Command line options:
37// - -r, --refine: number of uniform refinements of the mesh (default: 2)
38// - -l, --lambda: Lame parameter lambda (default: 50)
39// - -m, --mu: Lame parameter mu (default: 50)
40//
41// Compile with: see README.md
42//
43// Sample runs: mpirun -n 2 contact-patch-test
44// mpirun -n 2 contact-patch-test -r 3
45
46#include "mfem.hpp"
47
48#include "axom/slic.hpp"
49
50#include "tribol/interface/tribol.hpp"
51#include "tribol/interface/mfem_tribol.hpp"
52
53// Define MPI_REAL_T
54#if defined(MFEM_USE_DOUBLE)
55#define MPI_REAL_T MPI_DOUBLE
56#else
57#error "Tribol requires MFEM built with double precision!"
58#endif
59
60
61int main(int argc, char *argv[])
62{
63 // Initialize MPI
65
66 // Initialize logging with axom::slic
67 axom::slic::SimpleLogger logger;
68 axom::slic::setIsRoot(mfem::Mpi::Root());
69
70 // Define command line options
71 int ref_levels = 2; // number of times to uniformly refine the serial mesh
72 mfem::real_t lambda = 50.0; // Lame parameter lambda
73 mfem::real_t mu = 50.0; // Lame parameter mu (shear modulus)
74
75 // Parse command line options
76 mfem::OptionsParser args(argc, argv);
77 args.AddOption(&ref_levels, "-r", "--refine",
78 "Number of times to refine the mesh uniformly.");
79 args.AddOption(&lambda, "-l", "--lambda",
80 "Lame parameter lambda.");
81 args.AddOption(&mu, "-m", "--mu",
82 "Lame parameter mu (shear modulus).");
83 args.Parse();
84 if (!args.Good())
85 {
86 if (mfem::Mpi::Root())
87 {
88 args.PrintUsage(std::cout);
89 }
90 return EXIT_FAILURE;
91 }
92 if (mfem::Mpi::Root())
93 {
94 args.PrintOptions(std::cout);
95 }
96
97 // Fixed options
98 // two block mesh; bottom block = [0,1]^3 and top block = [0,1]x[0,1]x[0.99,1.99]
99 std::string mesh_file = "two-hex.mesh";
100 // Problem dimension (NOTE: Tribol's mortar only works in 3D)
101 constexpr int dim = 3;
102 // FE polynomial degree (NOTE: only 1 works for now)
103 constexpr int order = 1;
104 // z=1 plane of bottom block (contact plane)
105 std::set<int> mortar_attrs({4});
106 // z=0.99 plane of top block (contact plane)
107 std::set<int> nonmortar_attrs({5});
108 // per-dimension sets of boundary attributes with homogeneous Dirichlet BCs.
109 // allows transverse deformation of the blocks while precluding rigid body
110 // rotations/translations.
111 std::vector<std::set<int>> fixed_attrs(dim);
112 fixed_attrs[0] = {1}; // x=0 plane of both blocks
113 fixed_attrs[1] = {2}; // y=0 plane of both blocks
114 fixed_attrs[2] = {3, 6}; // 3: z=0 plane of bottom block; 6: z=1.99 plane of top block
115
116 // Read the mesh, refine, and create a mfem::ParMesh
117 mfem::Mesh serial_mesh(mesh_file);
118 for (int i = 0; i < ref_levels; ++i)
119 {
120 serial_mesh.UniformRefinement();
121 }
122 mfem::ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
123 serial_mesh.Clear();
124
125 MFEM_ASSERT(dim == mesh.Dimension(),
126 "This miniapp must be run with the supplied two-hex.mesh file.");
127
128 // Create an H1 finite element space on the mesh for displacements/forces
129 mfem::H1_FECollection fec(order, dim);
130 mfem::ParFiniteElementSpace fespace(&mesh, &fec, dim);
131 auto n_displacement_dofs = fespace.GlobalTrueVSize();
132 if (mfem::Mpi::Root())
133 {
134 std::cout << "Number of displacement unknowns: " << n_displacement_dofs <<
135 std::endl;
136 }
137
138 // Create coordinate and displacement grid functions
139 mfem::ParGridFunction coords(&fespace);
140 mesh.SetNodalGridFunction(&coords);
141 mfem::ParGridFunction displacement(&fespace);
142 displacement = 0.0;
143
144 // Find true dofs with homogeneous Dirichlet BCs
145 mfem::Array<int> ess_tdof_list;
146 {
147 mfem::Array<int> ess_vdof_marker(fespace.GetVSize());
148 ess_vdof_marker = 0;
149 for (int i = 0; i < dim; ++i)
150 {
151 mfem::Array<int> ess_bdr(mesh.bdr_attributes.Max());
152 ess_bdr = 0;
153 for (auto xfixed_attr : fixed_attrs[i])
154 {
155 ess_bdr[xfixed_attr-1] = 1;
156 }
157 mfem::Array<int> new_ess_vdof_marker;
158 fespace.GetEssentialVDofs(ess_bdr, new_ess_vdof_marker, i);
159 for (int j = 0; j < new_ess_vdof_marker.Size(); ++j)
160 {
161 ess_vdof_marker[j] = ess_vdof_marker[j] || new_ess_vdof_marker[j];
162 }
163 }
164 mfem::Array<int> ess_tdof_marker;
165 fespace.GetRestrictionMatrix()->BooleanMult(ess_vdof_marker, ess_tdof_marker);
166 mfem::FiniteElementSpace::MarkerToList(ess_tdof_marker, ess_tdof_list);
167 }
168
169 // Create small deformation linear elastic bilinear form
170 mfem::ParBilinearForm a(&fespace);
171 mfem::ConstantCoefficient lambda_coeff(lambda);
173 a.AddDomainIntegrator(new mfem::ElasticityIntegrator(lambda_coeff, mu_coeff));
174
175 // Compute elasticity contribution to tangent stiffness matrix
176 a.Assemble();
177 std::unique_ptr<mfem::HypreParMatrix> A(new mfem::HypreParMatrix);
178 a.FormSystemMatrix(ess_tdof_list, *A);
179
180 // #1: Initialize Tribol contact library
181 tribol::initialize(dim, MPI_COMM_WORLD);
182
183 // #2: Create a Tribol coupling scheme: defines contact surfaces and enforcement
184 int coupling_scheme_id = 0;
185 // NOTE: While there is a single mfem ParMesh for this problem, Tribol
186 // defines a mortar and a nonmortar contact mesh, each with a unique mesh ID.
187 // The Tribol mesh IDs for each contact surface are defined here.
188 int mesh1_id = 0;
189 int mesh2_id = 1;
190 tribol::registerMfemCouplingScheme(
191 coupling_scheme_id, mesh1_id, mesh2_id,
192 mesh, coords, mortar_attrs, nonmortar_attrs,
193 tribol::SURFACE_TO_SURFACE,
194 tribol::NO_CASE,
195 tribol::SINGLE_MORTAR,
196 tribol::FRICTIONLESS,
197 tribol::LAGRANGE_MULTIPLIER,
198 tribol::BINNING_GRID
199 );
200
201 // #3: Set additional options/access pressure grid function on contact surfaces
202 // Access Tribol's pressure grid function (on the contact surface). The
203 // pressure ParGridFunction is created upon calling
204 // registerMfemCouplingScheme(). It's lifetime coincides with the lifetime of
205 // the coupling scheme, so the host code can reference and update it as
206 // needed.
207 auto& pressure = tribol::getMfemPressure(coupling_scheme_id);
208 if (mfem::Mpi::Root())
209 {
210 std::cout << "Number of pressure unknowns: " <<
211 pressure.ParFESpace()->GlobalTrueVSize() << std::endl;
212 }
213
214 // Set Tribol options for Lagrange multiplier enforcement
215 tribol::setLagrangeMultiplierOptions(
216 coupling_scheme_id,
217 tribol::ImplicitEvalMode::MORTAR_RESIDUAL_JACOBIAN
218 );
219
220 // #4: Update contact mesh decomposition so the on-rank Tribol meshes
221 // coincide with the current configuration of the mesh. This must be called
222 // before tribol::update().
223 tribol::updateMfemParallelDecomposition();
224
225 // #5: Update contact gaps, forces, and tangent stiffness contributions
226 int cycle = 1; // pseudo cycle
227 mfem::real_t t = 1.0; // pseudo time
228 mfem::real_t dt = 1.0; // pseudo dt
229 tribol::update(cycle, t, dt);
230
231 // #6a: Return contact contribution to the tangent stiffness matrix as a
232 // block operator. See documentation for getMfemBlockJacobian() for block
233 // definitions.
234 auto A_blk = tribol::getMfemBlockJacobian(coupling_scheme_id);
235 // Add elasticity contribution to the block operator returned by
236 // getMfemBlockJacobian(). Note the (0,0) block Tribol returns is empty, so
237 // we can simply set the (0,0) block to A.
238 A_blk->SetBlock(0, 0, A.release());
239
240 // Convert block operator to a single HypreParMatrix
242 for (int i{0}; i < 2; ++i)
243 {
244 for (int j{0}; j < 2; ++j)
245 {
246 if (A_blk->GetBlock(i, j).Height() != 0 && A_blk->GetBlock(i, j).Width() != 0)
247 {
248 hypre_blocks(i, j) =
249 dynamic_cast<const mfem::HypreParMatrix*>(&A_blk->GetBlock(i, j));
250 }
251 else
252 {
253 hypre_blocks(i, j) = nullptr;
254 }
255 }
256 }
257 auto A_hyprePar = std::unique_ptr<mfem::HypreParMatrix>(
259 );
260
261 // Create block RHS vector holding forces and gaps at tdofs
262 mfem::BlockVector B_blk(A_blk->RowOffsets());
263 B_blk = 0.0;
264
265 // Fill with initial nodal gaps.
266 // Note forces from contact are currently zero since pressure is zero prior
267 // to first solve.
268 mfem::Vector gap;
269 // #6b: Return computed gap constraints on the contact surfaces
270 tribol::getMfemGap(coupling_scheme_id, gap); // gap on ldofs
271 auto& P_submesh = *pressure.ParFESpace()->GetProlongationMatrix();
272 auto& gap_true = B_blk.GetBlock(1); // gap tdof vectorParFESpace()
273 // gap is a dual vector, so (gap tdof vector) = P^T * (gap ldof vector)
274 P_submesh.MultTranspose(gap, gap_true);
275
276 // Create block solution vector holding displacements and pressures at tdofs
277 mfem::BlockVector X_blk(A_blk->ColOffsets());
278 X_blk = 0.0;
279
280 // Solve for displacements and pressures. Direct solvers such as STRUMPACK
281 // are a good way to solve the saddle point system, but MINRES can be used if
282 // the MFEM build doesn't have access to STRUMPACK.
283#ifdef MFEM_USE_STRUMPACK
284 std::unique_ptr<mfem::STRUMPACKRowLocMatrix> A_strumpk(new
285 mfem::STRUMPACKRowLocMatrix(*A_hyprePar));
286 mfem::STRUMPACKSolver solver(MPI_COMM_WORLD);
287 solver.SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
288 solver.SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
289 solver.SetPrintFactorStatistics(true);
290 solver.SetPrintSolveStatistics(true);
291 solver.SetOperator(*A_strumpk);
292#else
293 mfem::MINRESSolver solver(MPI_COMM_WORLD);
294 solver.SetRelTol(1.0e-12);
295 solver.SetMaxIter(2000);
296 solver.SetPrintLevel(3);
297 solver.SetOperator(*A_hyprePar);
298 mfem::HypreDiagScale prec(*A_hyprePar);
300 solver.SetPreconditioner(prec);
301#endif
302 solver.Mult(B_blk, X_blk);
303
304 // Update displacement and coords grid functions
305 auto& displacement_true = X_blk.GetBlock(0);
306 fespace.GetProlongationMatrix()->Mult(displacement_true, displacement);
307 displacement.Neg();
308 coords += displacement;
309
310 // Update the pressure grid function
311 auto& pressure_true = X_blk.GetBlock(1);
312 P_submesh.Mult(pressure_true, pressure);
313
314 // Verify the forces are in equilibrium, i.e. f_int = A*u = -f_contact = -B^T*p
315 // This should be true if the solver converges.
316 mfem::Vector f_int_true(fespace.GetTrueVSize());
317 f_int_true = 0.0;
318 mfem::Vector f_contact_true(f_int_true);
319 A_blk->GetBlock(0, 0).Mult(displacement_true, f_int_true);
320 A_blk->GetBlock(0, 1).Mult(pressure_true, f_contact_true);
321 mfem::Vector resid_true(f_int_true);
322 resid_true += f_contact_true;
323 for (int i{0}; i < ess_tdof_list.Size(); ++i)
324 {
325 resid_true[ess_tdof_list[i]] = 0.0;
326 }
327 auto resid_linf = resid_true.Normlinf();
328 if (mfem::Mpi::Root())
329 {
330 MPI_Reduce(MPI_IN_PLACE, &resid_linf, 1, MPI_REAL_T, MPI_MAX, 0,
331 MPI_COMM_WORLD);
332 std::cout << "|| force residual ||_(infty) = " << resid_linf << std::endl;
333 }
334 else
335 {
336 MPI_Reduce(&resid_linf, &resid_linf, 1, MPI_REAL_T, MPI_MAX, 0, MPI_COMM_WORLD);
337 }
338
339 // Verify the gap is closed by the displacements, i.e. B*u = gap
340 // This should be true if the solver converges.
341 mfem::Vector gap_resid_true(gap_true.Size());
342 gap_resid_true = 0.0;
343 A_blk->GetBlock(1, 0).Mult(displacement_true, gap_resid_true);
344 gap_resid_true -= gap_true;
345 auto gap_resid_linf = gap_resid_true.Normlinf();
346 if (mfem::Mpi::Root())
347 {
348 MPI_Reduce(MPI_IN_PLACE, &gap_resid_linf, 1, MPI_REAL_T, MPI_MAX, 0,
349 MPI_COMM_WORLD);
350 std::cout << "|| gap residual ||_(infty) = " << gap_resid_linf << std::endl;
351 }
352 else
353 {
354 MPI_Reduce(&gap_resid_linf, &gap_resid_linf, 1, MPI_REAL_T, MPI_MAX, 0,
355 MPI_COMM_WORLD);
356 }
357
358 // Update the Tribol mesh based on deformed configuration
359 // NOTE: This is done to update the contact submesh based on the current
360 // deformed shape of the mesh.
361 tribol::updateMfemParallelDecomposition();
362
363 // Save data in VisIt format
364 mfem::VisItDataCollection visit_vol_dc("contact-patch-test-volume", &mesh);
365 visit_vol_dc.RegisterField("coordinates", &coords);
366 visit_vol_dc.RegisterField("displacement", &displacement);
367 visit_vol_dc.Save();
368 mfem::VisItDataCollection visit_surf_dc("contact-patch-test-surface",
369 pressure.ParFESpace()->GetMesh());
370 visit_surf_dc.RegisterField("pressure", &pressure);
371 visit_surf_dc.Save();
372
373 // #7: Tribol cleanup: deletes coupling schemes and clears associated memory
374 tribol::finalize();
375
376 return 0;
377}
Dynamic 2D array using row-major layout.
Definition array.hpp:372
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
int Size() const
Return the logical size of the array.
Definition array.hpp:144
A class to handle Vectors in a block fashion.
Vector & GetBlock(int i)
Get the i-th vector in the block.
A coefficient that is constant across space and time.
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
Definition fespace.cpp:648
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Jacobi preconditioner in hypre.
Definition hypre.hpp:1481
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
@ IGNORE_HYPRE_ERRORS
Ignore hypre errors (see e.g. HypreADS)
Definition hypre.hpp:1167
void SetErrorMode(ErrorMode err_mode) const
Set the behavior for treating hypre errors, see the ErrorMode enum. The default mode in the base clas...
Definition hypre.hpp:1241
MINRES method.
Definition solvers.hpp:628
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition mesh.cpp:6200
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Class for parallel bilinear form.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
int GetTrueVSize() const override
Return the number of local vector true dofs.
Definition pfespace.hpp:289
const Operator * GetProlongationMatrix() const override
The returned Operator is owned by the FiniteElementSpace.
void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_dofs, int component=-1) const override
Determine the boundary degrees of freedom.
const SparseMatrix * GetRestrictionMatrix() const override
Get the R matrix which restricts a local dof vector to true dof vector.
Definition pfespace.hpp:389
Class for parallel grid function.
Definition pgridfunc.hpp:33
Class for parallel meshes.
Definition pmesh.hpp:34
void SetOperator(const Operator &op)
Set the operator/matrix.
void SetMaxIter(int max_it)
Set the maximum number of iterations for iterative solvers.
void SetPrintFactorStatistics(bool print_stat)
Set up verbose printing during the factor step.
void SetKrylovSolver(strumpack::KrylovSolver method)
Set the Krylov solver method to use.
void Mult(const Vector &x, Vector &y) const
Factor and solve the linear system .
void SetRelTol(double rtol)
Set the relative tolerance for interative solvers.
void SetPrintSolveStatistics(bool print_stat)
Set up verbose printing during the solve step.
void SetReorderingStrategy(strumpack::ReorderingStrategy method)
Set matrix reordering strategy.
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
Vector data type.
Definition vector.hpp:80
void Neg()
(*this) = -(*this)
Definition vector.cpp:300
real_t Normlinf() const
Returns the l_infinity norm of the vector.
Definition vector.cpp:850
Data collection with VisIt I/O routines.
virtual void Save() override
Save the collection and a VisIt root file.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
int dim
Definition ex24.cpp:53
real_t mu
Definition ex25.cpp:140
int main()
real_t a
Definition lissajous.cpp:41
HypreParMatrix * HypreParMatrixFromBlocks(Array2D< const HypreParMatrix * > &blocks, Array2D< real_t > *blockCoeff)
Returns a merged hypre matrix constructed from hypre matrix blocks.
Definition hypre.cpp:3127
float real_t
Definition config.hpp:43
RefCoord t[3]