61int main(
int argc,
char *argv[])
67 axom::slic::SimpleLogger logger;
77 args.
AddOption(&ref_levels,
"-r",
"--refine",
78 "Number of times to refine the mesh uniformly.");
80 "Lame parameter lambda.");
82 "Lame parameter mu (shear modulus).");
99 std::string mesh_file =
"two-hex.mesh";
101 constexpr int dim = 3;
103 constexpr int order = 1;
105 std::set<int> mortar_attrs({4});
107 std::set<int> nonmortar_attrs({5});
111 std::vector<std::set<int>> fixed_attrs(
dim);
112 fixed_attrs[0] = {1};
113 fixed_attrs[1] = {2};
114 fixed_attrs[2] = {3, 6};
118 for (
int i = 0; i < ref_levels; ++i)
126 "This miniapp must be run with the supplied two-hex.mesh file.");
134 std::cout <<
"Number of displacement unknowns: " << n_displacement_dofs <<
149 for (
int i = 0; i <
dim; ++i)
153 for (
auto xfixed_attr : fixed_attrs[i])
155 ess_bdr[xfixed_attr-1] = 1;
159 for (
int j = 0; j < new_ess_vdof_marker.
Size(); ++j)
161 ess_vdof_marker[j] = ess_vdof_marker[j] || new_ess_vdof_marker[j];
178 a.FormSystemMatrix(ess_tdof_list, *A);
181 tribol::initialize(
dim, MPI_COMM_WORLD);
184 int coupling_scheme_id = 0;
190 tribol::registerMfemCouplingScheme(
191 coupling_scheme_id, mesh1_id, mesh2_id,
192 mesh, coords, mortar_attrs, nonmortar_attrs,
193 tribol::SURFACE_TO_SURFACE,
195 tribol::SINGLE_MORTAR,
196 tribol::FRICTIONLESS,
197 tribol::LAGRANGE_MULTIPLIER,
207 auto& pressure = tribol::getMfemPressure(coupling_scheme_id);
210 std::cout <<
"Number of pressure unknowns: " <<
211 pressure.ParFESpace()->GlobalTrueVSize() << std::endl;
215 tribol::setLagrangeMultiplierOptions(
217 tribol::ImplicitEvalMode::MORTAR_RESIDUAL_JACOBIAN
223 tribol::updateMfemParallelDecomposition();
229 tribol::update(cycle,
t, dt);
234 auto A_blk = tribol::getMfemBlockJacobian(coupling_scheme_id);
238 A_blk->SetBlock(0, 0, A.release());
242 for (
int i{0}; i < 2; ++i)
244 for (
int j{0}; j < 2; ++j)
246 if (A_blk->GetBlock(i, j).Height() != 0 && A_blk->GetBlock(i, j).Width() != 0)
253 hypre_blocks(i, j) =
nullptr;
257 auto A_hyprePar = std::unique_ptr<mfem::HypreParMatrix>(
270 tribol::getMfemGap(coupling_scheme_id, gap);
271 auto& P_submesh = *pressure.ParFESpace()->GetProlongationMatrix();
274 P_submesh.MultTranspose(gap, gap_true);
283#ifdef MFEM_USE_STRUMPACK
284 std::unique_ptr<mfem::STRUMPACKRowLocMatrix> A_strumpk(
new
296 solver.SetPrintLevel(3);
300 solver.SetPreconditioner(prec);
302 solver.
Mult(B_blk, X_blk);
305 auto& displacement_true = X_blk.
GetBlock(0);
308 coords += displacement;
311 auto& pressure_true = X_blk.
GetBlock(1);
312 P_submesh.Mult(pressure_true, pressure);
319 A_blk->GetBlock(0, 0).Mult(displacement_true, f_int_true);
320 A_blk->GetBlock(0, 1).Mult(pressure_true, f_contact_true);
322 resid_true += f_contact_true;
323 for (
int i{0}; i < ess_tdof_list.
Size(); ++i)
325 resid_true[ess_tdof_list[i]] = 0.0;
327 auto resid_linf = resid_true.
Normlinf();
330 MPI_Reduce(MPI_IN_PLACE, &resid_linf, 1, MPI_REAL_T, MPI_MAX, 0,
332 std::cout <<
"|| force residual ||_(infty) = " << resid_linf << std::endl;
336 MPI_Reduce(&resid_linf, &resid_linf, 1, MPI_REAL_T, MPI_MAX, 0, MPI_COMM_WORLD);
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();
348 MPI_Reduce(MPI_IN_PLACE, &gap_resid_linf, 1, MPI_REAL_T, MPI_MAX, 0,
350 std::cout <<
"|| gap residual ||_(infty) = " << gap_resid_linf << std::endl;
354 MPI_Reduce(&gap_resid_linf, &gap_resid_linf, 1, MPI_REAL_T, MPI_MAX, 0,
361 tribol::updateMfemParallelDecomposition();
369 pressure.ParFESpace()->GetMesh());
371 visit_surf_dc.
Save();