Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "ContactAction.h"
11 :
12 : #include "Factory.h"
13 : #include "FEProblem.h"
14 : #include "Conversion.h"
15 : #include "AddVariableAction.h"
16 : #include "MortarConstraintBase.h"
17 : #include "MortarConsumerInterface.h"
18 : #include "NonlinearSystemBase.h"
19 : #include "Parser.h"
20 : #include "AugmentedLagrangianContactProblem.h"
21 :
22 : #include "NanoflannMeshAdaptor.h"
23 : #include "PointListAdaptor.h"
24 :
25 : #include <set>
26 : #include <algorithm>
27 : #include <unordered_map>
28 : #include <limits>
29 :
30 : #include "libmesh/petsc_nonlinear_solver.h"
31 : #include "libmesh/string_to_enum.h"
32 :
33 : // Make newer nanoflann API compatible with older nanoflann versions
34 : #if NANOFLANN_VERSION < 0x150
35 : namespace nanoflann
36 : {
37 : typedef SearchParams SearchParameters;
38 :
39 : template <typename T, typename U>
40 : using ResultItem = std::pair<T, U>;
41 : }
42 : #endif
43 :
44 : using NodeBoundaryIDInfo = std::pair<const Node *, BoundaryID>;
45 :
46 : // Counter for naming mortar auxiliary kernels
47 : static unsigned int contact_mortar_auxkernel_counter = 0;
48 :
49 : // Counter for naming auxiliary kernels
50 : static unsigned int contact_auxkernel_counter = 0;
51 :
52 : // Counter for naming nodal area user objects
53 : static unsigned int contact_userobject_counter = 0;
54 :
55 : // Counter for distinct contact action objects
56 : static unsigned int contact_action_counter = 0;
57 :
58 : // For mortar subdomains
59 : registerMooseAction("ContactApp", ContactAction, "append_mesh_generator");
60 : registerMooseAction("ContactApp", ContactAction, "add_aux_variable");
61 : // For mortar Lagrange multiplier
62 : registerMooseAction("ContactApp", ContactAction, "add_contact_aux_variable");
63 : registerMooseAction("ContactApp", ContactAction, "add_mortar_variable");
64 : registerMooseAction("ContactApp", ContactAction, "add_aux_kernel");
65 : // For mortar constraint
66 : registerMooseAction("ContactApp", ContactAction, "add_constraint");
67 : registerMooseAction("ContactApp", ContactAction, "output_penetration_info_vars");
68 : registerMooseAction("ContactApp", ContactAction, "add_user_object");
69 : // For automatic generation of contact pairs
70 : registerMooseAction("ContactApp", ContactAction, "post_mesh_prepared");
71 :
72 : InputParameters
73 1593 : ContactAction::validParams()
74 : {
75 1593 : InputParameters params = Action::validParams();
76 1593 : params += ContactAction::commonParameters();
77 :
78 3186 : params.addParam<std::vector<BoundaryName>>(
79 : "primary", "The list of boundary IDs referring to primary sidesets");
80 3186 : params.addParam<std::vector<BoundaryName>>(
81 : "secondary", "The list of boundary IDs referring to secondary sidesets");
82 3186 : params.addParam<std::vector<BoundaryName>>(
83 : "automatic_pairing_boundaries",
84 : {},
85 : "List of boundary IDs for sidesets that are automatically paired with any other boundary in "
86 : "this list having a centroid-to-centroid distance less than the value specified in the "
87 : "'automatic_pairing_distance' parameter. ");
88 3186 : params.addRangeCheckedParam<Real>(
89 : "automatic_pairing_distance",
90 : "automatic_pairing_distance>=0",
91 : "The maximum distance the centroids of the boundaries provided in the "
92 : "'automatic_pairing_boundaries' parameter can be to generate a contact pair automatically. "
93 : "Due to numerical error in the determination of the centroids, it is encouraged that "
94 : "the user adds a tolerance to this distance (e.g. extra 10%) to make sure no suitable "
95 : "contact pair is missed. If the 'automatic_pairing_method = NODE' option is chosen instead, "
96 : "this distance is recommended to be set to at least twice the minimum distance between "
97 : "nodes of boundaries to be paired.");
98 3186 : params.addDeprecatedParam<MeshGeneratorName>(
99 : "mesh",
100 : "The mesh generator for mortar method",
101 : "This parameter is not used anymore and can simply be removed");
102 3186 : params.addParam<VariableName>("secondary_gap_offset",
103 : "Offset to gap distance from secondary side");
104 3186 : params.addParam<VariableName>("mapped_primary_gap_offset",
105 : "Offset to gap distance mapped from primary side");
106 3186 : params.addParam<std::vector<VariableName>>(
107 : "displacements",
108 : {},
109 : "The displacements appropriate for the simulation geometry and coordinate system");
110 3186 : params.addParam<Real>(
111 : "penalty",
112 3186 : 1e8,
113 : "The penalty to apply. This can vary depending on the stiffness of your materials");
114 3186 : params.addParam<Real>(
115 : "penalty_friction",
116 3186 : 1e8,
117 : "The penalty factor to apply in mortar penalty frictional constraints. It is applied to the "
118 : "tangential accumulated slip to build the frictional force");
119 4779 : params.addRangeCheckedParam<Real>(
120 : "penalty_multiplier",
121 3186 : 1.0,
122 : "penalty_multiplier > 0",
123 : "The growth factor for the penalty applied at the end of each augmented "
124 : "Lagrange update iteration (a value larger than one, e.g., 10, tends to speed up "
125 : "convergence.)");
126 4779 : params.addRangeCheckedParam<Real>(
127 : "penalty_multiplier_friction",
128 3186 : 1.0,
129 : "penalty_multiplier_friction > 0",
130 : "The penalty growth factor between augmented Lagrange "
131 : "iterations for penalizing relative slip distance if the node is under stick conditions.(a "
132 : "value larger than one, e.g., 10, tends to speed up convergence.)");
133 3186 : params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
134 3186 : params.addParam<Real>("tension_release",
135 3186 : 0.0,
136 : "Tension release threshold. A node in contact "
137 : "will not be released if its tensile load is below "
138 : "this value. No tension release if negative.");
139 3186 : params.addParam<MooseEnum>("model", ContactAction::getModelEnum(), "The contact model to use");
140 3186 : params.addParam<Real>("tangential_tolerance",
141 : "Tangential distance to extend edges of contact surfaces");
142 3186 : params.addParam<Real>("capture_tolerance",
143 3186 : 0.0,
144 : "Normal distance from surface within which nodes are captured. This "
145 : "parameter is used for node-face and mortar formulations.");
146 3186 : params.addParam<Real>(
147 : "normal_smoothing_distance",
148 : "Distance from edge in parametric coordinates over which to smooth contact normal");
149 :
150 3186 : params.addParam<bool>("normalize_penalty",
151 3186 : false,
152 : "Whether to normalize the penalty parameter with the nodal area.");
153 3186 : params.addParam<bool>(
154 : "primary_secondary_jacobian",
155 3186 : true,
156 : "Whether to include Jacobian entries coupling primary and secondary nodes.");
157 3186 : params.addParam<bool>(
158 : "ghost_whole_interface",
159 3186 : false,
160 : "Whether to geometrically and algebraically ghost the entire primary side of the interface "
161 : "for node-face contact constraints.");
162 3186 : params.addParam<Real>("al_penetration_tolerance",
163 : "The tolerance of the penetration for augmented Lagrangian method.");
164 3186 : params.addParam<Real>("al_incremental_slip_tolerance",
165 : "The tolerance of the incremental slip for augmented Lagrangian method.");
166 4779 : params.addRangeCheckedParam<Real>(
167 : "max_penalty_multiplier",
168 3186 : 1.0e3,
169 : "max_penalty_multiplier >= 1.0",
170 : "Maximum multiplier applied to penalty factors when adaptivity is used in an augmented "
171 : "Lagrange setting. The penalty factor supplied by the user is used as a reference to "
172 : "determine its maximum. If this multiplier is too large, the condition number of the system "
173 : "to be solved may be negatively impacted.");
174 3186 : MooseEnum adaptivity_penalty_normal("SIMPLE BUSSETTA", "SIMPLE");
175 3186 : adaptivity_penalty_normal.addDocumentation(
176 : "SIMPLE", "Keep multiplying by the penalty multiplier between AL iterations");
177 3186 : adaptivity_penalty_normal.addDocumentation(
178 : "BUSSETTA",
179 : "Modify the penalty using an algorithm from Bussetta et al, 2012, Comput Mech 49:259-275 "
180 : "between AL iterations.");
181 3186 : params.addParam<MooseEnum>(
182 : "adaptivity_penalty_normal",
183 : adaptivity_penalty_normal,
184 : "The augmented Lagrange update strategy used on the normal penalty coefficient.");
185 3186 : MooseEnum adaptivity_penalty_friction("SIMPLE FRICTION_LIMIT", "FRICTION_LIMIT");
186 3186 : adaptivity_penalty_friction.addDocumentation(
187 : "SIMPLE", "Keep multiplying by the frictional penalty multiplier between AL iterations");
188 3186 : adaptivity_penalty_friction.addDocumentation(
189 : "FRICTION_LIMIT",
190 : "This strategy will be guided by the Coulomb limit and be less reliant on the initial "
191 : "penalty factor provided by the user.");
192 3186 : params.addParam<MooseEnum>(
193 : "adaptivity_penalty_friction",
194 : adaptivity_penalty_friction,
195 : "The augmented Lagrange update strategy used on the frictional penalty coefficient.");
196 3186 : params.addParam<Real>("al_frictional_force_tolerance",
197 : "The tolerance of the frictional force for augmented Lagrangian method.");
198 3186 : params.addParam<Real>(
199 : "c_normal",
200 3186 : 1e6,
201 : "Parameter for balancing the size of the gap and contact pressure for a mortar formulation. "
202 : "This purely numerical "
203 : "parameter affects convergence behavior and, in general, should be larger for stiffer "
204 : "materials. It is recommended that the user tries out various orders of magnitude for this "
205 : "parameter if the default value generates poor contact convergence.");
206 3186 : params.addParam<Real>(
207 3186 : "c_tangential", 1, "Numerical parameter for nonlinear mortar frictional constraints");
208 3186 : params.addParam<bool>("ping_pong_protection",
209 3186 : false,
210 : "Whether to protect against ping-ponging, e.g. the oscillation of the "
211 : "secondary node between two "
212 : "different primary faces, by tying the secondary node to the "
213 : "edge between the involved primary faces");
214 3186 : params.addParam<Real>(
215 : "normal_lm_scaling",
216 3186 : 1.,
217 : "Scaling factor to apply to the normal LM variable for a mortar formulation");
218 3186 : params.addParam<Real>(
219 : "tangential_lm_scaling",
220 3186 : 1.,
221 : "Scaling factor to apply to the tangential LM variable for a mortar formulation");
222 3186 : MooseEnum lm_space(getContactLMSpaceOptions(), "MATCH_DISPLACEMENT");
223 3186 : lm_space.addDocumentation(
224 : "MATCH_DISPLACEMENT",
225 : "Use the same finite element order as the displacement variables for generated mortar "
226 : "Lagrange multiplier variables.");
227 3186 : lm_space.addDocumentation(
228 : "LINEAR",
229 : "Use first-order LAGRANGE generated mortar Lagrange multiplier variables, independent "
230 : "of the displacement variable order.");
231 3186 : params.addParam<MooseEnum>(
232 : "lm_space",
233 : lm_space,
234 : "Finite element space for mortar Lagrange multiplier variables generated by the "
235 : "contact action. This parameter only applies to the 'mortar' contact formulation.");
236 3186 : params.addParam<bool>(
237 : "normalize_c",
238 3186 : false,
239 : "Whether to normalize c by weighting function norm for mortar contact. When unnormalized "
240 : "the value of c effectively depends on element size since in the constraint we compare nodal "
241 : "Lagrange Multiplier values to integrated gap values (LM nodal value is independent of "
242 : "element size, where integrated values are dependent on element size).");
243 1593 : params.addClassDescription("Sets up all objects needed for mechanical contact enforcement");
244 3186 : params.addParam<bool>(
245 : "use_dual",
246 : "Whether to use the dual mortar approach within a mortar formulation. It is defaulted to "
247 : "true for "
248 : "weighted quantity approach, and to false for the legacy approach. To avoid instabilities "
249 : "in the solution and obtain the full benefits of a variational enforcement,"
250 : "use of dual mortar with weighted constraints is strongly recommended. This "
251 : "input is only intended for advanced users.");
252 3186 : params.addParam<bool>(
253 : "correct_edge_dropping",
254 3186 : false,
255 : "Whether to enable correct edge dropping treatment for mortar constraints. When disabled "
256 : "any Lagrange Multiplier degree of freedom on a secondary element without full primary "
257 : "contributions will be set (strongly) to 0.");
258 1593 : params += MortarConsumerInterface::triangulationParams();
259 3186 : params.addParam<bool>(
260 : "generate_mortar_mesh",
261 3186 : true,
262 : "Whether to generate the mortar mesh from the action. Typically this will be the case, but "
263 : "one may also want to reuse an existing lower-dimensional mesh prior to a restart.");
264 3186 : params.addParam<MooseEnum>("automatic_pairing_method",
265 3186 : ContactAction::getProximityMethod(),
266 : "The proximity method used for automatic pairing of boundaries.");
267 3186 : params.addParam<bool>(
268 : "mortar_dynamics",
269 3186 : false,
270 : "Whether to use constraints that account for the persistency condition, giving rise to "
271 : "smoother normal contact pressure evolution. This flag should only be set to yes for dynamic "
272 : "simulations using the Newmark-beta numerical integrator");
273 3186 : params.addParam<Real>(
274 : "newmark_beta",
275 3186 : 0.25,
276 : "Newmark-beta beta parameter for its inclusion in the weighted gap update formula");
277 3186 : params.addParam<Real>(
278 : "newmark_gamma",
279 3186 : 0.5,
280 : "Newmark-beta gamma parameter for its inclusion in the weighted gap update formula");
281 3186 : params.addCoupledVar("wear_depth",
282 : "The name of the mortar auxiliary variable that is used to modify the "
283 : "weighted gap definition");
284 3186 : params.addParam<std::vector<TagName>>(
285 : "extra_vector_tags",
286 : "The tag names for extra vectors that residual data should be saved into");
287 3186 : params.addParam<std::vector<TagName>>(
288 : "absolute_value_vector_tags",
289 : "The tags for the vectors this residual object should fill with the "
290 : "absolute value of the residual contribution");
291 3186 : params.addParam<bool>(
292 : "use_petrov_galerkin",
293 3186 : false,
294 : "Whether to use the Petrov-Galerkin approach for the mortar-based constraints. If set to "
295 : "true, we use the standard basis as the test function and dual basis as "
296 : "the shape function for the interpolation of the Lagrange multiplier variable.");
297 3186 : params.addParam<bool>(
298 : "debug_mesh",
299 3186 : false,
300 : "Whether we are going to enable mortar segment mesh debug information. An exodus"
301 : "file will be generated if the user sets this flag to true");
302 1593 : const auto mortar_constraint_params = MortarConstraintBase::validParams();
303 3186 : params.transferParam<MooseEnum>(mortar_constraint_params, "segment_quadrature");
304 3186 : params.transferParam<Real>(mortar_constraint_params, "minimum_projection_angle");
305 3186 : params.transferParam<MooseEnum>(mortar_constraint_params, "mortar_3d_subpatch_plane");
306 3186 : params.transferParam<MooseEnum>(mortar_constraint_params, "mortar_3d_qp_mapping");
307 :
308 : // Contact surface definition
309 3186 : params.addParamNamesToGroup("primary secondary displacements", "Contact Surface Definition");
310 : // Automatic pairing
311 3186 : params.addParamNamesToGroup(
312 : "automatic_pairing_boundaries automatic_pairing_distance automatic_pairing_method",
313 : "Automatic Contact Pair Generation");
314 : // Contact formulation and model
315 3186 : params.addParamNamesToGroup("formulation model", "Contact Formulation");
316 : // Penalty parameters
317 3186 : params.addParamNamesToGroup(
318 : "penalty penalty_friction penalty_multiplier penalty_multiplier_friction "
319 : "max_penalty_multiplier normalize_penalty",
320 : "Penalty Parameters");
321 : // Augmented Lagrange settings
322 3186 : params.addParamNamesToGroup(
323 : "al_penetration_tolerance al_incremental_slip_tolerance al_frictional_force_tolerance "
324 : "adaptivity_penalty_normal adaptivity_penalty_friction",
325 : "Augmented Lagrange");
326 : // Friction
327 3186 : params.addParamNamesToGroup("friction_coefficient tension_release", "Friction");
328 : // Mortar-specific parameters
329 3186 : params.addParamNamesToGroup("c_normal c_tangential normal_lm_scaling tangential_lm_scaling "
330 : "lm_space "
331 : "use_dual correct_edge_dropping normalize_c use_petrov_galerkin "
332 : "generate_mortar_mesh segment_quadrature minimum_projection_angle "
333 : "mortar_3d_subpatch_plane mortar_3d_qp_mapping wear_depth debug_mesh",
334 : "Mortar");
335 : // Mortar dynamics (Newmark-beta)
336 3186 : params.addParamNamesToGroup("mortar_dynamics newmark_beta newmark_gamma", "Mortar Dynamics");
337 : // Gap and tolerance settings
338 3186 : params.addParamNamesToGroup(
339 : "secondary_gap_offset mapped_primary_gap_offset capture_tolerance "
340 : "tangential_tolerance normal_smoothing_distance normal_smoothing_method",
341 : "Gap and Tolerance");
342 : // Jacobian and solver options
343 3186 : params.addParamNamesToGroup("primary_secondary_jacobian ping_pong_protection", "Solver Options");
344 : // Interface ghosting
345 3186 : params.addParamNamesToGroup("ghost_whole_interface", "Interface Ghosting");
346 : // Residual vector tags
347 3186 : params.addParamNamesToGroup("extra_vector_tags absolute_value_vector_tags", "Residual Tags");
348 :
349 1593 : return params;
350 1593 : }
351 :
352 1593 : ContactAction::ContactAction(const InputParameters & params)
353 : : Action(params),
354 3186 : _boundary_pairs(getParam<BoundaryName, BoundaryName>("primary", "secondary")),
355 3186 : _model(getParam<MooseEnum>("model").getEnum<ContactModel>()),
356 3186 : _formulation(getParam<MooseEnum>("formulation").getEnum<ContactFormulation>()),
357 3186 : _lm_space(getParam<MooseEnum>("lm_space").getEnum<ContactLMSpace>()),
358 3186 : _generate_mortar_mesh(getParam<bool>("generate_mortar_mesh")),
359 4779 : _mortar_dynamics(getParam<bool>("mortar_dynamics"))
360 : {
361 : // Check for automatic selection of contact pairs.
362 3186 : if (getParam<std::vector<BoundaryName>>("automatic_pairing_boundaries").size() > 1)
363 : _automatic_pairing_boundaries =
364 84 : getParam<std::vector<BoundaryName>>("automatic_pairing_boundaries");
365 :
366 1677 : if (_automatic_pairing_boundaries.size() > 0 && !isParamValid("automatic_pairing_distance"))
367 0 : paramError("automatic_pairing_distance",
368 : "For automatic selection of contact pairs (for particular geometries) in contact "
369 : "action, 'automatic_pairing_distance' needs to be provided.");
370 :
371 1677 : if (_automatic_pairing_boundaries.size() > 0 && !isParamValid("automatic_pairing_method"))
372 0 : paramError("automatic_pairing_distance",
373 : "For automatic selection of contact pairs (for particular geometries) in contact "
374 : "action, 'automatic_pairing_method' needs to be provided.");
375 :
376 1593 : if (_automatic_pairing_boundaries.size() > 0 && _boundary_pairs.size() != 0)
377 0 : paramError("automatic_pairing_boundaries",
378 : "If a boundary list is provided, primary and secondary surfaces will be identified "
379 : "automatically. Therefore, one cannot provide an automatic pairing boundary list "
380 : "and primary/secondary lists.");
381 1593 : else if (_automatic_pairing_boundaries.size() == 0 && _boundary_pairs.size() == 0)
382 2 : paramError("primary",
383 : "'primary' and 'secondary' surfaces or a list of boundaries for automatic pair "
384 : "generation need to be provided.");
385 :
386 : // End of checks for automatic selection of contact pairs.
387 :
388 1591 : if (_boundary_pairs.size() != 1 && _formulation == ContactFormulation::MORTAR)
389 0 : paramError("formulation", "When using mortar, a vector of contact pairs cannot be used");
390 :
391 3182 : if ((_formulation == ContactFormulation::MORTAR ||
392 1591 : _formulation == ContactFormulation::MORTAR_PENALTY) &&
393 2006 : params.isParamSetByUser("ghost_whole_interface"))
394 4 : paramError("ghost_whole_interface",
395 : "The 'ghost_whole_interface' parameter is only supported for node-face contact "
396 : "formulations. Mortar contact always geometrically and algebraically ghosts the "
397 : "interface.");
398 :
399 1587 : if (_formulation == ContactFormulation::TANGENTIAL_PENALTY && _model != ContactModel::COULOMB)
400 0 : paramError("formulation",
401 : "The 'tangential_penalty' formulation can only be used with the 'coulomb' model");
402 :
403 1587 : if (_formulation == ContactFormulation::MORTAR_PENALTY)
404 : {
405 : // Use dual basis functions for contact traction interpolation
406 114 : if (isParamValid("use_dual"))
407 40 : _use_dual = getParam<bool>("use_dual");
408 : else
409 37 : _use_dual = true;
410 :
411 57 : if (_model == ContactModel::GLUED)
412 0 : paramError("model", "The 'mortar_penalty' formulation does not support glued contact");
413 :
414 114 : if (getParam<bool>("mortar_dynamics"))
415 0 : paramError("mortar_dynamics",
416 : "The 'mortar_penalty' formulation does not support implicit dynamic simulations");
417 :
418 114 : if (getParam<bool>("use_petrov_galerkin"))
419 0 : paramError("use_petrov_galerkin",
420 : "The 'mortar_penalty' formulation does not support usage of the Petrov-Galerkin "
421 : "flag. The default (use_dual = true) behavior is such that contact tractions are "
422 : "interpolated with dual bases whereas mortar or weighted contact quantities are "
423 : "interpolated with Lagrange shape functions.");
424 : }
425 :
426 1587 : if (_formulation == ContactFormulation::MORTAR)
427 : {
428 354 : if (_model == ContactModel::GLUED)
429 0 : paramError("model", "The 'mortar' formulation does not support glued contact (yet)");
430 :
431 : // use dual basis function for Lagrange multipliers?
432 708 : if (isParamValid("use_dual"))
433 82 : _use_dual = getParam<bool>("use_dual");
434 : else
435 313 : _use_dual = true;
436 :
437 708 : if (!getParam<bool>("mortar_dynamics"))
438 : {
439 666 : if (params.isParamSetByUser("newmark_beta"))
440 2 : paramError("newmark_beta", "newmark_beta can only be used with the mortar_dynamics option");
441 :
442 662 : if (params.isParamSetByUser("newmark_gamma"))
443 2 : paramError("newmark_gamma",
444 : "newmark_gamma can only be used with the mortar_dynamics option");
445 : }
446 :
447 700 : if (isParamSetByUser("penalty"))
448 0 : paramError("penalty",
449 : "The 'penalty' parameter is not used for the 'mortar' formulation which instead "
450 : "uses Lagrange multipliers");
451 : }
452 : else
453 : {
454 2466 : if (params.isParamSetByUser("correct_edge_dropping"))
455 0 : paramError(
456 : "correct_edge_dropping",
457 : "The 'correct_edge_dropping' option can only be used with the 'mortar' formulation "
458 : "(weighted)");
459 2466 : else if (params.isParamSetByUser("triangulation") &&
460 0 : _formulation != ContactFormulation::MORTAR_PENALTY)
461 0 : paramError("triangulation",
462 : "The 'triangulation' option can only be used with mortar-based formulations.");
463 2466 : else if (params.isParamSetByUser("triangulate_triangles") &&
464 0 : _formulation != ContactFormulation::MORTAR_PENALTY)
465 0 : paramError("triangulate_triangles",
466 : "The 'triangulate_triangles' option can only be used with mortar-based "
467 : "formulations.");
468 2466 : else if (params.isParamSetByUser("minimum_projection_angle") &&
469 6 : _formulation != ContactFormulation::MORTAR_PENALTY)
470 0 : paramError("minimum_projection_angle",
471 : "The 'minimum_projection_angle' option can only be used with mortar-based "
472 : "formulations.");
473 2466 : else if (params.isParamSetByUser("mortar_3d_subpatch_plane") &&
474 8 : _formulation != ContactFormulation::MORTAR_PENALTY)
475 2 : paramError("mortar_3d_subpatch_plane",
476 : "The 'mortar_3d_subpatch_plane' option can only be used with mortar-based "
477 : "formulations.");
478 2462 : else if (params.isParamSetByUser("mortar_3d_qp_mapping") &&
479 9 : _formulation != ContactFormulation::MORTAR_PENALTY)
480 2 : paramError("mortar_3d_qp_mapping",
481 : "The 'mortar_3d_qp_mapping' option can only be used with mortar-based "
482 : "formulations.");
483 2458 : else if (params.isParamSetByUser("use_dual") &&
484 20 : _formulation != ContactFormulation::MORTAR_PENALTY)
485 0 : paramError("use_dual",
486 : "The 'use_dual' option can only be used with the 'mortar' formulation");
487 2458 : else if (params.isParamSetByUser("c_normal"))
488 0 : paramError("c_normal",
489 : "The 'c_normal' option can only be used with the 'mortar' formulation");
490 2458 : else if (params.isParamSetByUser("c_tangential"))
491 0 : paramError("c_tangential",
492 : "The 'c_tangential' option can only be used with the 'mortar' formulation");
493 2458 : else if (params.isParamSetByUser("mortar_dynamics"))
494 0 : paramError("mortar_dynamics",
495 : "The 'mortar_dynamics' constraint option can only be used with the 'mortar' "
496 : "formulation and in dynamic simulations using Newmark-beta");
497 2458 : else if (params.isParamSetByUser("segment_quadrature"))
498 0 : paramError("segment_quadrature",
499 : "The 'segment_quadrature' option can only be used with the "
500 : "'mortar' formulation.");
501 2458 : else if (params.isParamSetByUser("lm_space"))
502 4 : paramError("lm_space",
503 : "The 'lm_space' option can only be used with the 'mortar' formulation.");
504 : }
505 :
506 1575 : if (_formulation == ContactFormulation::RANFS)
507 : {
508 18 : if (isParamValid("secondary_gap_offset"))
509 0 : paramError("secondary_gap_offset",
510 : "The 'secondary_gap_offset' option can only be used with the "
511 : "'MechanicalContactConstraint'");
512 18 : if (isParamValid("mapped_primary_gap_offset"))
513 0 : paramError("mapped_primary_gap_offset",
514 : "The 'mapped_primary_gap_offset' option can only be used with the "
515 : "'MechanicalContactConstraint'");
516 : }
517 3132 : else if (getParam<bool>("ping_pong_protection"))
518 4 : paramError("ping_pong_protection",
519 : "The 'ping_pong_protection' option can only be used with the 'ranfs' formulation");
520 :
521 : // Remove repeated pairs from input file.
522 1575 : removeRepeatedPairs();
523 1587 : }
524 :
525 : void
526 1603 : ContactAction::removeRepeatedPairs()
527 : {
528 1603 : if (_boundary_pairs.size() == 0 && _automatic_pairing_boundaries.size() == 0)
529 0 : paramError(
530 : "primary",
531 : "Number of contact pairs in the contact action is zero. Please revise your input file.");
532 :
533 : // Remove repeated interactions
534 : std::vector<std::pair<BoundaryName, BoundaryName>> lean_boundary_pairs;
535 :
536 10527 : for (const auto & [primary, secondary] : _boundary_pairs)
537 : {
538 : // Structured bindings are not capturable (primary_copy, secondary_copy)
539 8924 : auto it = std::find_if(lean_boundary_pairs.begin(),
540 : lean_boundary_pairs.end(),
541 17848 : [&, primary_copy = primary, secondary_copy = secondary](
542 : const std::pair<BoundaryName, BoundaryName> & lean_pair)
543 : {
544 14646 : const bool match_one = lean_pair.second == secondary_copy &&
545 9689 : lean_pair.first == primary_copy;
546 14646 : const bool match_two = lean_pair.second == primary_copy &&
547 9 : lean_pair.first == secondary_copy;
548 14646 : const bool exist = match_one || match_two;
549 14646 : return exist;
550 : });
551 :
552 8924 : if (it == lean_boundary_pairs.end())
553 1682 : lean_boundary_pairs.emplace_back(primary, secondary);
554 : else
555 7242 : mooseInfo("Contact pair ",
556 : primary,
557 : "--",
558 : secondary,
559 : " has been removed from the contact interaction list due to "
560 : "duplicates in the input file.");
561 : }
562 :
563 1603 : _boundary_pairs = lean_boundary_pairs;
564 1603 : }
565 :
566 : void
567 14125 : ContactAction::act()
568 : {
569 : // proform problem checks/corrections once during the first feasible task
570 14125 : if (_current_task == "add_contact_aux_variable")
571 : {
572 3146 : if (!_problem->getDisplacedProblem())
573 2 : mooseError(
574 : "Contact requires updated coordinates. Use the 'displacements = ...' parameter in the "
575 : "Mesh block.");
576 :
577 : // It is risky to apply this optimization to contact problems
578 : // since the problem configuration may be changed during Jacobian
579 : // evaluation. We therefore turn it off for all contact problems so that
580 : // PETSc-3.8.4 or higher will have the same behavior as PETSc-3.8.3.
581 1571 : if (!_problem->isSNESMFReuseBaseSetbyUser())
582 : _problem->setSNESMFReuseBase(false, false);
583 : }
584 :
585 14123 : if (_formulation == ContactFormulation::MORTAR ||
586 : _formulation == ContactFormulation::MORTAR_PENALTY)
587 3621 : addMortarContact();
588 : else
589 10502 : addNodeFaceContact();
590 :
591 14117 : if (_current_task == "add_aux_kernel")
592 : {
593 3126 : if (!_problem->getDisplacedProblem())
594 0 : mooseError("Contact requires updated coordinates. Use the 'displacements = ...' line in the "
595 : "Mesh block.");
596 :
597 : // Create auxiliary kernels for each contact pairs
598 3233 : for (const auto & contact_pair : _boundary_pairs)
599 : {
600 : const auto & [primary_name, secondary_name] = contact_pair;
601 1670 : if ((_formulation != ContactFormulation::MORTAR) &&
602 : (_formulation != ContactFormulation::MORTAR_PENALTY))
603 : {
604 1273 : InputParameters params = _factory.getValidParams("PenetrationAux");
605 1273 : params.applyParameters(parameters(),
606 : {"secondary_gap_offset", "mapped_primary_gap_offset", "order"});
607 :
608 : std::vector<VariableName> displacements =
609 2546 : getParam<std::vector<VariableName>>("displacements");
610 1273 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
611 1273 : .system()
612 1273 : .variable_type(displacements[0])
613 : .order.get_order();
614 :
615 2546 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
616 5092 : params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_LINEAR};
617 3819 : params.set<std::vector<BoundaryName>>("boundary") = {secondary_name};
618 2546 : params.set<BoundaryName>("paired_boundary") = primary_name;
619 2546 : params.set<AuxVariableName>("variable") = "penetration";
620 2546 : if (isParamValid("secondary_gap_offset"))
621 12 : params.set<std::vector<VariableName>>("secondary_gap_offset") = {
622 24 : getParam<VariableName>("secondary_gap_offset")};
623 2546 : if (isParamValid("mapped_primary_gap_offset"))
624 12 : params.set<std::vector<VariableName>>("mapped_primary_gap_offset") = {
625 24 : getParam<VariableName>("mapped_primary_gap_offset")};
626 1273 : params.set<bool>("use_displaced_mesh") = true;
627 2546 : std::string name = _name + "_contact_" + Moose::stringify(contact_auxkernel_counter++);
628 :
629 2546 : _problem->addAuxKernel("PenetrationAux", name, params);
630 1273 : }
631 : else
632 : {
633 : const auto type = "MortarUserObjectAux";
634 794 : InputParameters params = _factory.getValidParams(type);
635 1191 : params.set<std::vector<BoundaryName>>("boundary") = {secondary_name};
636 794 : params.set<AuxVariableName>("variable") = "gap";
637 397 : params.set<bool>("use_displaced_mesh") = true; // Unnecessary as this object only operates
638 : // on nodes, but we'll do it for consistency
639 794 : params.set<MooseEnum>("contact_quantity") = "normal_gap";
640 : const auto & [primary_id, secondary_id, uo_name] =
641 397 : libmesh_map_find(_bnd_pair_to_mortar_info, contact_pair);
642 794 : params.set<UserObjectName>("user_object") = uo_name;
643 794 : std::string name = _name + "_contact_gap_" + std::to_string(primary_id) + "_" +
644 794 : std::to_string(secondary_id);
645 :
646 794 : _problem->addAuxKernel(type, name, params);
647 397 : }
648 : }
649 :
650 1563 : addContactPressureAuxKernel();
651 :
652 3126 : const unsigned int ndisp = getParam<std::vector<VariableName>>("displacements").size();
653 :
654 : // Add MortarFrictionalPressureVectorAux
655 1563 : if (_formulation == ContactFormulation::MORTAR && _model == ContactModel::COULOMB && ndisp > 2)
656 : {
657 : {
658 28 : InputParameters params = _factory.getValidParams("MortarFrictionalPressureVectorAux");
659 :
660 28 : params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
661 28 : params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
662 42 : params.set<std::vector<BoundaryName>>("boundary") = {_boundary_pairs[0].second};
663 42 : params.set<ExecFlagEnum>("execute_on", true) = {EXEC_NONLINEAR};
664 :
665 14 : std::string action_name = MooseUtils::shortName(name());
666 14 : const std::string tangential_lagrange_multiplier_name = action_name + "_tangential_lm";
667 : const std::string tangential_lagrange_multiplier_3d_name =
668 14 : action_name + "_tangential_3d_lm";
669 :
670 28 : params.set<std::vector<VariableName>>("tangent_one") = {
671 56 : tangential_lagrange_multiplier_name};
672 28 : params.set<std::vector<VariableName>>("tangent_two") = {
673 56 : tangential_lagrange_multiplier_3d_name};
674 :
675 56 : std::vector<std::string> disp_components({"x", "y", "z"});
676 : unsigned component_index = 0;
677 :
678 : // Loop over three displacements
679 56 : for (const auto & disp_component : disp_components)
680 : {
681 126 : params.set<AuxVariableName>("variable") = _name + "_tangent_" + disp_component;
682 42 : params.set<unsigned int>("component") = component_index;
683 :
684 84 : std::string name = _name + "_mortar_frictional_pressure_" + disp_component + "_" +
685 84 : Moose::stringify(contact_mortar_auxkernel_counter++);
686 :
687 42 : _problem->addAuxKernel("MortarFrictionalPressureVectorAux", name, params);
688 42 : component_index++;
689 : }
690 28 : }
691 : }
692 : }
693 :
694 14117 : if (_current_task == "add_contact_aux_variable")
695 : {
696 3142 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
697 1571 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
698 1571 : .system()
699 1571 : .variable_type(displacements[0])
700 : .order.get_order();
701 : const auto mortar_lm_order =
702 1571 : _lm_space == ContactLMSpace::LINEAR ? static_cast<int>(FIRST) : order;
703 1571 : std::unique_ptr<InputParameters> current_params;
704 : const auto create_aux_var_params =
705 4308 : [this, order, mortar_lm_order, ¤t_params]() -> InputParameters &
706 : {
707 4308 : current_params = std::make_unique<InputParameters>(_factory.getValidParams("MooseVariable"));
708 : // Node/face and mortar-penalty contact aux variables continue to follow the displacement
709 : // order. Mortar LM contact aux variables live on the same contact surface as the generated LM
710 : // field, so they use the selected generated LM space.
711 4308 : const auto aux_order = _formulation == ContactFormulation::MORTAR ? mortar_lm_order : order;
712 8616 : current_params->set<MooseEnum>("order") =
713 8616 : Utility::enum_to_string<Order>(OrderWrapper{aux_order});
714 8616 : current_params->set<MooseEnum>("family") = "LAGRANGE";
715 4308 : return *current_params;
716 1571 : };
717 :
718 1571 : if ((_formulation != ContactFormulation::MORTAR) &&
719 : (_formulation != ContactFormulation::MORTAR_PENALTY))
720 : {
721 : // Add penetration aux variable
722 2332 : _problem->addAuxVariable("MooseVariable", "penetration", create_aux_var_params());
723 : // Add nodal area aux variable
724 2332 : _problem->addAuxVariable("MooseVariable", "nodal_area", create_aux_var_params());
725 : }
726 : else
727 810 : _problem->addAuxVariable("MooseVariable", "gap", create_aux_var_params());
728 :
729 : // Add contact pressure aux variable
730 3142 : _problem->addAuxVariable("MooseVariable", "contact_pressure", create_aux_var_params());
731 :
732 3142 : const unsigned int ndisp = getParam<std::vector<VariableName>>("displacements").size();
733 :
734 : // Add MortarFrictionalPressureVectorAux variables
735 1571 : if (_formulation == ContactFormulation::MORTAR && _model == ContactModel::COULOMB && ndisp > 2)
736 : {
737 : {
738 56 : std::vector<std::string> disp_components({"x", "y", "z"});
739 : // Loop over three displacements
740 56 : for (const auto & disp_component : disp_components)
741 : {
742 42 : auto var_params = _factory.getValidParams("MooseVariable");
743 84 : var_params.set<MooseEnum>("order") =
744 84 : Utility::enum_to_string<Order>(OrderWrapper{mortar_lm_order});
745 84 : var_params.set<MooseEnum>("family") = "LAGRANGE";
746 :
747 126 : _problem->addAuxVariable(
748 42 : "MooseVariable", _name + "_tangent_" + disp_component, var_params);
749 42 : }
750 14 : }
751 : }
752 1571 : }
753 :
754 14117 : if (_current_task == "add_user_object" && (_formulation != ContactFormulation::MORTAR) &&
755 : (_formulation != ContactFormulation::MORTAR_PENALTY))
756 : {
757 1166 : auto var_params = _factory.getValidParams("NodalArea");
758 :
759 : // Get secondary_boundary_vector from possibly updated set from the
760 : // ContactAction constructor cleanup
761 1166 : const auto actions = _awh.getActions<ContactAction>();
762 :
763 : std::vector<BoundaryName> secondary_boundary_vector;
764 2440 : for (const auto * const action : actions)
765 2655 : for (const auto j : index_range(action->_boundary_pairs))
766 1381 : secondary_boundary_vector.push_back(action->_boundary_pairs[j].second);
767 :
768 2332 : var_params.set<std::vector<BoundaryName>>("boundary") = secondary_boundary_vector;
769 3498 : var_params.set<std::vector<VariableName>>("variable") = {"nodal_area"};
770 :
771 : mooseAssert(_problem, "Problem pointer is NULL");
772 4664 : var_params.set<ExecFlagEnum>("execute_on", true) = {EXEC_INITIAL, EXEC_TIMESTEP_BEGIN};
773 1166 : var_params.set<bool>("use_displaced_mesh") = true;
774 :
775 3498 : _problem->addUserObject("NodalArea",
776 1166 : "nodal_area_object_" + Moose::stringify(contact_userobject_counter++),
777 : var_params);
778 1166 : }
779 16650 : }
780 :
781 : void
782 1563 : ContactAction::addContactPressureAuxKernel()
783 : {
784 : // Increment counter for contact action objects
785 1563 : contact_action_counter++;
786 :
787 1563 : if ((_formulation != ContactFormulation::MORTAR) &&
788 : (_formulation != ContactFormulation::MORTAR_PENALTY))
789 : {
790 : // Add ContactPressureAux: Only one object for all contact pairs
791 1166 : const auto actions = _awh.getActions<ContactAction>();
792 :
793 : // Add auxiliary kernel if we are the last contact action object.
794 1166 : if (contact_action_counter == actions.size())
795 : {
796 : std::vector<BoundaryName> boundary_vector;
797 : std::vector<BoundaryName> pair_boundary_vector;
798 :
799 2296 : for (const auto * const action : actions)
800 2439 : for (const auto j : index_range(action->_boundary_pairs))
801 : {
802 1273 : boundary_vector.push_back(action->_boundary_pairs[j].second);
803 1273 : pair_boundary_vector.push_back(action->_boundary_pairs[j].first);
804 : }
805 :
806 1130 : InputParameters params = _factory.getValidParams("ContactPressureAux");
807 1130 : params.applyParameters(parameters(), {"order"});
808 :
809 : std::vector<VariableName> displacements =
810 2260 : getParam<std::vector<VariableName>>("displacements");
811 1130 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
812 1130 : .system()
813 1130 : .variable_type(displacements[0])
814 : .order.get_order();
815 :
816 2260 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
817 1130 : params.set<std::vector<BoundaryName>>("boundary") = boundary_vector;
818 2260 : params.set<std::vector<BoundaryName>>("paired_boundary") = pair_boundary_vector;
819 2260 : params.set<AuxVariableName>("variable") = "contact_pressure";
820 2260 : params.addRequiredCoupledVar("nodal_area", "The nodal area");
821 3390 : params.set<std::vector<VariableName>>("nodal_area") = {"nodal_area"};
822 1130 : params.set<bool>("use_displaced_mesh") = true;
823 :
824 1130 : std::string name = _name + "_contact_pressure";
825 2260 : params.set<ExecFlagEnum>("execute_on",
826 6780 : true) = {EXEC_NONLINEAR, EXEC_TIMESTEP_END, EXEC_TIMESTEP_BEGIN};
827 2260 : _problem->addAuxKernel("ContactPressureAux", name, params);
828 1130 : }
829 1166 : }
830 : else
831 794 : for (const auto & contact_pair : _boundary_pairs)
832 : {
833 : const auto & [_, secondary_name] = contact_pair;
834 : const auto type = "MortarUserObjectAux";
835 794 : InputParameters params = _factory.getValidParams(type);
836 1191 : params.set<std::vector<BoundaryName>>("boundary") = {secondary_name};
837 794 : params.set<AuxVariableName>("variable") = "contact_pressure";
838 397 : params.set<bool>("use_displaced_mesh") = true; // Unecessary as this object only operates on
839 : // nodes, but we'll do it for consistency
840 794 : params.set<MooseEnum>("contact_quantity") = "normal_pressure";
841 : const auto & [primary_id, secondary_id, uo_name] =
842 397 : libmesh_map_find(_bnd_pair_to_mortar_info, contact_pair);
843 794 : params.set<UserObjectName>("user_object") = uo_name;
844 794 : const std::string name = _name + "_contact_pressure" + std::to_string(primary_id) + "_" +
845 794 : std::to_string(secondary_id);
846 :
847 794 : _problem->addAuxKernel(type, name, params);
848 397 : }
849 2693 : }
850 :
851 : void
852 4701 : ContactAction::addRelationshipManagers(Moose::RelationshipManagerType input_rm_type)
853 : {
854 4701 : if (_formulation == ContactFormulation::MORTAR ||
855 : _formulation == ContactFormulation::MORTAR_PENALTY)
856 : {
857 1199 : auto params = MortarConstraintBase::validParams();
858 1199 : params.set<bool>("use_displaced_mesh") = true;
859 1199 : std::string action_name = MooseUtils::shortName(name());
860 1199 : const std::string primary_subdomain_name = action_name + "_primary_subdomain";
861 1199 : const std::string secondary_subdomain_name = action_name + "_secondary_subdomain";
862 2398 : params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
863 2398 : params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
864 2398 : params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
865 2398 : params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
866 2398 : params.set<bool>("use_petrov_galerkin") = getParam<bool>("use_petrov_galerkin");
867 2398 : params.set<Real>("minimum_projection_angle") = getParam<Real>("minimum_projection_angle");
868 2398 : params.set<MooseEnum>("mortar_3d_subpatch_plane") =
869 2398 : getParam<MooseEnum>("mortar_3d_subpatch_plane");
870 : const bool augmented_penalty =
871 1199 : _formulation == ContactFormulation::MORTAR_PENALTY &&
872 157 : dynamic_cast<AugmentedLagrangianContactProblemInterface *>(_problem.get());
873 1288 : params.set<bool>("ghost_point_neighbors") = !_mortar_dynamics && !augmented_penalty;
874 1199 : addRelationshipManagers(input_rm_type, params);
875 1199 : }
876 : else
877 : {
878 : const std::string constraint_type = _formulation == ContactFormulation::RANFS
879 : ? "RANFSNormalMechanicalContact"
880 6977 : : "MechanicalContactConstraint";
881 :
882 7262 : for (const auto & contact_pair : _boundary_pairs)
883 : {
884 3760 : auto params = _factory.getValidParams(constraint_type);
885 3760 : params.set<bool>("use_displaced_mesh") = true;
886 7520 : params.set<bool>("ghost_whole_interface") = getParam<bool>("ghost_whole_interface");
887 3760 : params.set<BoundaryName>("primary") = contact_pair.first;
888 3760 : params.set<BoundaryName>("secondary") = contact_pair.second;
889 3760 : addRelationshipManagers(input_rm_type, params);
890 3760 : }
891 : }
892 4701 : }
893 :
894 : void
895 3621 : ContactAction::addMortarContact()
896 : {
897 3621 : std::string action_name = MooseUtils::shortName(name());
898 :
899 10863 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
900 3621 : const unsigned int ndisp = displacements.size();
901 :
902 : // Definitions for mortar contact.
903 3621 : const std::string primary_subdomain_name = action_name + "_primary_subdomain";
904 3621 : const std::string secondary_subdomain_name = action_name + "_secondary_subdomain";
905 3621 : const std::string normal_lagrange_multiplier_name = action_name + "_normal_lm";
906 3621 : const std::string tangential_lagrange_multiplier_name = action_name + "_tangential_lm";
907 3621 : const std::string tangential_lagrange_multiplier_3d_name = action_name + "_tangential_3d_lm";
908 3621 : const std::string auxiliary_lagrange_multiplier_name = action_name + "_aux_lm";
909 :
910 3621 : if (_current_task == "append_mesh_generator")
911 : {
912 : // Don't do mesh generators when recovering or when the user has requested for us not to
913 : // (presumably because the lower-dimensional blocks are already in the mesh due to manual
914 : // addition or because we are restarting)
915 405 : if (!(_app.isRecovering() && _app.isUltimateMaster()) && !_app.useMasterMesh() &&
916 378 : _generate_mortar_mesh)
917 : {
918 371 : const MeshGeneratorName primary_name = primary_subdomain_name + "_generator";
919 371 : const MeshGeneratorName secondary_name = secondary_subdomain_name + "_generator";
920 :
921 371 : auto primary_params = _factory.getValidParams("LowerDBlockFromSidesetGenerator");
922 742 : auto secondary_params = _factory.getValidParams("LowerDBlockFromSidesetGenerator");
923 :
924 742 : primary_params.set<SubdomainName>("new_block_name") = primary_subdomain_name;
925 742 : secondary_params.set<SubdomainName>("new_block_name") = secondary_subdomain_name;
926 :
927 1113 : primary_params.set<std::vector<BoundaryName>>("sidesets") = {_boundary_pairs[0].first};
928 1113 : secondary_params.set<std::vector<BoundaryName>>("sidesets") = {_boundary_pairs[0].second};
929 :
930 742 : _app.appendMeshGenerator("LowerDBlockFromSidesetGenerator", primary_name, primary_params);
931 742 : _app.appendMeshGenerator("LowerDBlockFromSidesetGenerator", secondary_name, secondary_params);
932 371 : }
933 : }
934 :
935 : // Add the lagrange multiplier on the secondary subdomain.
936 : const auto addLagrangeMultiplier =
937 533 : [this, &secondary_subdomain_name, &displacements](const std::string & variable_name,
938 : const Real scaling_factor,
939 : const bool add_aux_lm,
940 : const bool penalty_traction) //
941 : {
942 533 : InputParameters params = _factory.getValidParams("MooseVariableBase");
943 :
944 : // Allow the user to select "weighted" constraints and standard bases (use_dual = false) or
945 : // "legacy" constraints and dual bases (use_dual = true). Unless it's for testing purposes,
946 : // this combination isn't recommended
947 533 : if (!add_aux_lm || penalty_traction)
948 519 : params.set<bool>("use_dual") = _use_dual;
949 :
950 : mooseAssert(_problem->systemBaseNonlinear(/*nl_sys_num=*/0).hasVariable(displacements[0]),
951 : "Displacement variable is missing");
952 : const auto primal_type =
953 533 : _problem->systemBaseNonlinear(/*nl_sys_num=*/0).system().variable_type(displacements[0]);
954 :
955 : // The lm_space option is only valid for the mortar Lagrange multiplier formulation. Mortar
956 : // penalty traction variables continue to use the displacement order.
957 : const int lm_order =
958 493 : _formulation == ContactFormulation::MORTAR && _lm_space == ContactLMSpace::LINEAR
959 533 : ? static_cast<int>(FIRST)
960 : : primal_type.order.get_order();
961 :
962 533 : if (primal_type.family == LAGRANGE)
963 : {
964 1066 : params.set<MooseEnum>("family") = Utility::enum_to_string<FEFamily>(primal_type.family);
965 1066 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{lm_order});
966 : }
967 : else
968 0 : mooseError("Invalid bases for mortar contact.");
969 :
970 1599 : params.set<std::vector<SubdomainName>>("block") = {secondary_subdomain_name};
971 533 : if (!(add_aux_lm || penalty_traction))
972 958 : params.set<std::vector<Real>>("scaling") = {scaling_factor};
973 :
974 533 : auto fe_type = AddVariableAction::feType(params);
975 533 : auto var_type = AddVariableAction::variableType(fe_type);
976 533 : if (add_aux_lm || penalty_traction)
977 54 : _problem->addAuxVariable(var_type, variable_name, params);
978 : else
979 479 : _problem->addVariable(var_type, variable_name, params);
980 1066 : };
981 :
982 3621 : if (_current_task == "add_mortar_variable" && _formulation == ContactFormulation::MORTAR)
983 : {
984 350 : addLagrangeMultiplier(
985 350 : normal_lagrange_multiplier_name, getParam<Real>("normal_lm_scaling"), false, false);
986 :
987 350 : if (_model == ContactModel::COULOMB)
988 : {
989 115 : addLagrangeMultiplier(tangential_lagrange_multiplier_name,
990 115 : getParam<Real>("tangential_lm_scaling"),
991 : false,
992 : false);
993 115 : if (ndisp > 2)
994 14 : addLagrangeMultiplier(tangential_lagrange_multiplier_3d_name,
995 28 : getParam<Real>("tangential_lm_scaling"),
996 : false,
997 : false);
998 : }
999 :
1000 700 : if (getParam<bool>("use_petrov_galerkin"))
1001 14 : addLagrangeMultiplier(auxiliary_lagrange_multiplier_name, 1.0, true, false);
1002 : }
1003 3271 : else if (_current_task == "add_mortar_variable" &&
1004 55 : _formulation == ContactFormulation::MORTAR_PENALTY)
1005 : {
1006 55 : if (_use_dual)
1007 40 : addLagrangeMultiplier(auxiliary_lagrange_multiplier_name, 1.0, false, true);
1008 : }
1009 :
1010 3621 : if (_current_task == "add_user_object")
1011 : {
1012 403 : const auto register_mortar_uo_name = [this](const auto & bnd_pair, const auto & uo_prefix)
1013 : {
1014 : const auto & [primary_name, secondary_name] = bnd_pair;
1015 403 : const auto primary_id = _mesh->getBoundaryID(primary_name);
1016 403 : const auto secondary_id = _mesh->getBoundaryID(secondary_name);
1017 403 : const auto uo_name = uo_prefix + name();
1018 403 : _bnd_pair_to_mortar_info.emplace(bnd_pair, MortarInfo{primary_id, secondary_id, uo_name});
1019 403 : return uo_name;
1020 808 : };
1021 :
1022 : // check if the correct problem class is selected if AL parameters are provided
1023 405 : if (_formulation == ContactFormulation::MORTAR_PENALTY &&
1024 55 : !dynamic_cast<AugmentedLagrangianContactProblemInterface *>(_problem.get()))
1025 : {
1026 : const std::vector<std::string> params = {"penalty_multiplier",
1027 : "penalty_multiplier_friction",
1028 : "al_penetration_tolerance",
1029 : "al_incremental_slip_tolerance",
1030 42 : "al_frictional_force_tolerance"};
1031 242 : for (const auto & param : params)
1032 202 : if (parameters().isParamSetByUser(param))
1033 2 : paramError(param,
1034 : "Augmented Lagrange parameter was specified, but the selected problem type "
1035 : "does not support Augmented Lagrange iterations.");
1036 40 : }
1037 :
1038 403 : if (_model != ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR)
1039 : {
1040 470 : auto uo_params = _factory.getValidParams("LMWeightedGapUserObject");
1041 :
1042 470 : uo_params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
1043 470 : uo_params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
1044 470 : uo_params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
1045 470 : uo_params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
1046 705 : uo_params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
1047 705 : uo_params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
1048 235 : if (ndisp > 2)
1049 216 : uo_params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
1050 235 : uo_params.set<bool>("use_displaced_mesh") = true;
1051 705 : uo_params.set<std::vector<VariableName>>("lm_variable") = {normal_lagrange_multiplier_name};
1052 235 : uo_params.applySpecificParameters(parameters(),
1053 : {"correct_edge_dropping",
1054 : "triangulation",
1055 : "triangulate_triangles",
1056 : "minimum_projection_angle",
1057 : "mortar_3d_subpatch_plane",
1058 : "mortar_3d_qp_mapping",
1059 : "use_petrov_galerkin",
1060 : "debug_mesh"});
1061 470 : if (getParam<bool>("use_petrov_galerkin"))
1062 21 : uo_params.set<std::vector<VariableName>>("aux_lm") = {auxiliary_lagrange_multiplier_name};
1063 :
1064 468 : _problem->addUserObject("LMWeightedGapUserObject",
1065 235 : register_mortar_uo_name(_boundary_pairs[0], "lm_weightedgap_object_"),
1066 : uo_params);
1067 233 : }
1068 168 : else if (_model == ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR)
1069 : {
1070 230 : auto uo_params = _factory.getValidParams("LMWeightedVelocitiesUserObject");
1071 230 : uo_params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
1072 230 : uo_params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
1073 230 : uo_params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
1074 230 : uo_params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
1075 345 : uo_params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
1076 345 : uo_params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
1077 115 : if (ndisp > 2)
1078 42 : uo_params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
1079 :
1080 115 : uo_params.set<VariableName>("secondary_variable") = displacements[0];
1081 115 : uo_params.set<bool>("use_displaced_mesh") = true;
1082 230 : uo_params.set<std::vector<VariableName>>("lm_variable_normal") = {
1083 460 : normal_lagrange_multiplier_name};
1084 230 : uo_params.set<std::vector<VariableName>>("lm_variable_tangential_one") = {
1085 460 : tangential_lagrange_multiplier_name};
1086 115 : if (ndisp > 2)
1087 28 : uo_params.set<std::vector<VariableName>>("lm_variable_tangential_two") = {
1088 56 : tangential_lagrange_multiplier_3d_name};
1089 115 : uo_params.applySpecificParameters(parameters(),
1090 : {"correct_edge_dropping",
1091 : "triangulation",
1092 : "triangulate_triangles",
1093 : "minimum_projection_angle",
1094 : "mortar_3d_subpatch_plane",
1095 : "mortar_3d_qp_mapping",
1096 : "use_petrov_galerkin",
1097 : "debug_mesh"});
1098 230 : if (getParam<bool>("use_petrov_galerkin"))
1099 21 : uo_params.set<std::vector<VariableName>>("aux_lm") = {auxiliary_lagrange_multiplier_name};
1100 :
1101 115 : const auto uo_name = _problem->addUserObject(
1102 : "LMWeightedVelocitiesUserObject",
1103 115 : register_mortar_uo_name(_boundary_pairs[0], "lm_weightedvelocities_object_"),
1104 230 : uo_params);
1105 115 : }
1106 :
1107 401 : if (_model != ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR_PENALTY)
1108 : {
1109 56 : auto uo_params = _factory.getValidParams("PenaltyWeightedGapUserObject");
1110 :
1111 56 : uo_params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
1112 56 : uo_params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
1113 56 : uo_params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
1114 56 : uo_params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
1115 84 : uo_params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
1116 84 : uo_params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
1117 :
1118 : // AL parameters
1119 28 : uo_params.applySpecificParameters(parameters(),
1120 : {"correct_edge_dropping",
1121 : "triangulation",
1122 : "triangulate_triangles",
1123 : "minimum_projection_angle",
1124 : "mortar_3d_subpatch_plane",
1125 : "mortar_3d_qp_mapping",
1126 : "penalty",
1127 : "debug_mesh",
1128 : "max_penalty_multiplier",
1129 : "adaptivity_penalty_normal"});
1130 :
1131 56 : if (isParamValid("al_penetration_tolerance"))
1132 0 : uo_params.set<Real>("penetration_tolerance") = getParam<Real>("al_penetration_tolerance");
1133 56 : if (isParamValid("penalty_multiplier"))
1134 56 : uo_params.set<Real>("penalty_multiplier") = getParam<Real>("penalty_multiplier");
1135 : // In the contact action, we force the physical value of the normal gap, which also normalizes
1136 : // the penalty factor with the "area" around the node
1137 28 : uo_params.set<bool>("use_physical_gap") = true;
1138 :
1139 28 : if (_use_dual)
1140 60 : uo_params.set<std::vector<VariableName>>("aux_lm") = {auxiliary_lagrange_multiplier_name};
1141 :
1142 28 : if (ndisp > 2)
1143 39 : uo_params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
1144 28 : uo_params.set<bool>("use_displaced_mesh") = true;
1145 :
1146 54 : _problem->addUserObject(
1147 : "PenaltyWeightedGapUserObject",
1148 28 : register_mortar_uo_name(_boundary_pairs[0], "penalty_weightedgap_object_"),
1149 : uo_params);
1150 26 : _problem->haveADObjects(true);
1151 26 : }
1152 373 : else if (_model == ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR_PENALTY)
1153 : {
1154 50 : auto uo_params = _factory.getValidParams("PenaltyFrictionUserObject");
1155 50 : uo_params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
1156 50 : uo_params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
1157 50 : uo_params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
1158 50 : uo_params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
1159 75 : uo_params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
1160 50 : uo_params.set<bool>("correct_edge_dropping") = getParam<bool>("correct_edge_dropping");
1161 75 : uo_params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
1162 25 : if (ndisp > 2)
1163 15 : uo_params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
1164 :
1165 25 : uo_params.set<VariableName>("secondary_variable") = displacements[0];
1166 25 : uo_params.set<bool>("use_displaced_mesh") = true;
1167 50 : uo_params.set<Real>("friction_coefficient") = getParam<Real>("friction_coefficient");
1168 50 : uo_params.set<Real>("penalty") = getParam<Real>("penalty");
1169 50 : uo_params.set<Real>("penalty_friction") = getParam<Real>("penalty_friction");
1170 :
1171 : // AL parameters
1172 50 : uo_params.set<Real>("max_penalty_multiplier") = getParam<Real>("max_penalty_multiplier");
1173 50 : uo_params.set<MooseEnum>("adaptivity_penalty_normal") =
1174 50 : getParam<MooseEnum>("adaptivity_penalty_normal");
1175 50 : uo_params.set<MooseEnum>("adaptivity_penalty_friction") =
1176 50 : getParam<MooseEnum>("adaptivity_penalty_friction");
1177 50 : if (isParamValid("al_penetration_tolerance"))
1178 26 : uo_params.set<Real>("penetration_tolerance") = getParam<Real>("al_penetration_tolerance");
1179 50 : if (isParamValid("penalty_multiplier"))
1180 50 : uo_params.set<Real>("penalty_multiplier") = getParam<Real>("penalty_multiplier");
1181 50 : if (isParamValid("penalty_multiplier_friction"))
1182 25 : uo_params.set<Real>("penalty_multiplier_friction") =
1183 75 : getParam<Real>("penalty_multiplier_friction");
1184 :
1185 50 : if (isParamValid("al_incremental_slip_tolerance"))
1186 26 : uo_params.set<Real>("slip_tolerance") = getParam<Real>("al_incremental_slip_tolerance");
1187 : // In the contact action, we force the physical value of the normal gap, which also normalizes
1188 : // the penalty factor with the "area" around the node
1189 25 : uo_params.set<bool>("use_physical_gap") = true;
1190 :
1191 25 : if (_use_dual)
1192 54 : uo_params.set<std::vector<VariableName>>("aux_lm") = {auxiliary_lagrange_multiplier_name};
1193 :
1194 25 : uo_params.applySpecificParameters(parameters(),
1195 : {"triangulation",
1196 : "triangulate_triangles",
1197 : "minimum_projection_angle",
1198 : "mortar_3d_subpatch_plane",
1199 : "mortar_3d_qp_mapping",
1200 : "friction_coefficient",
1201 : "penalty",
1202 : "penalty_friction"});
1203 :
1204 50 : _problem->addUserObject(
1205 : "PenaltyFrictionUserObject",
1206 25 : register_mortar_uo_name(_boundary_pairs[0], "penalty_friction_object_"),
1207 : uo_params);
1208 25 : _problem->haveADObjects(true);
1209 25 : }
1210 : }
1211 :
1212 3615 : if (_current_task == "add_constraint")
1213 : {
1214 : // Prepare problem for enforcement with Lagrange multipliers
1215 397 : if (_model != ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR)
1216 : {
1217 : std::string mortar_constraint_name;
1218 :
1219 231 : if (!_mortar_dynamics)
1220 : mortar_constraint_name = "ComputeWeightedGapLMMechanicalContact";
1221 : else
1222 : mortar_constraint_name = "ComputeDynamicWeightedGapLMMechanicalContact";
1223 :
1224 231 : InputParameters params = _factory.getValidParams(mortar_constraint_name);
1225 231 : if (_mortar_dynamics)
1226 7 : params.applySpecificParameters(
1227 : parameters(), {"newmark_beta", "newmark_gamma", "capture_tolerance", "wear_depth"});
1228 :
1229 : else // We need user objects for quasistatic constraints
1230 896 : params.set<UserObjectName>("weighted_gap_uo") = "lm_weightedgap_object_" + name();
1231 :
1232 462 : params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
1233 462 : params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
1234 462 : params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
1235 462 : params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
1236 462 : params.set<NonlinearVariableName>("variable") = normal_lagrange_multiplier_name;
1237 693 : params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
1238 462 : params.set<Real>("c") = getParam<Real>("c_normal");
1239 :
1240 231 : if (ndisp > 1)
1241 693 : params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
1242 231 : if (ndisp > 2)
1243 204 : params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
1244 :
1245 231 : params.set<bool>("use_displaced_mesh") = true;
1246 :
1247 231 : params.applySpecificParameters(parameters(),
1248 : {"correct_edge_dropping",
1249 : "triangulation",
1250 : "triangulate_triangles",
1251 : "minimum_projection_angle",
1252 : "mortar_3d_subpatch_plane",
1253 : "mortar_3d_qp_mapping",
1254 : "normalize_c",
1255 : "extra_vector_tags",
1256 : "absolute_value_vector_tags",
1257 : "debug_mesh"});
1258 :
1259 231 : _problem->addConstraint(
1260 231 : mortar_constraint_name, action_name + "_normal_lm_weighted_gap", params);
1261 231 : _problem->haveADObjects(true);
1262 462 : }
1263 : // Add the tangential and normal Lagrange's multiplier constraints on the secondary boundary.
1264 166 : else if (_model == ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR)
1265 : {
1266 : std::string mortar_constraint_name;
1267 :
1268 115 : if (!_mortar_dynamics)
1269 : mortar_constraint_name = "ComputeFrictionalForceLMMechanicalContact";
1270 : else
1271 : mortar_constraint_name = "ComputeDynamicFrictionalForceLMMechanicalContact";
1272 :
1273 115 : InputParameters params = _factory.getValidParams(mortar_constraint_name);
1274 115 : if (_mortar_dynamics)
1275 14 : params.applySpecificParameters(
1276 : parameters(), {"newmark_beta", "newmark_gamma", "capture_tolerance", "wear_depth"});
1277 : else
1278 : { // We need user objects for quasistatic constraints
1279 404 : params.set<UserObjectName>("weighted_gap_uo") = "lm_weightedvelocities_object_" + name();
1280 202 : params.set<UserObjectName>("weighted_velocities_uo") =
1281 202 : "lm_weightedvelocities_object_" + name();
1282 : }
1283 :
1284 230 : params.set<bool>("correct_edge_dropping") = getParam<bool>("correct_edge_dropping");
1285 230 : params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
1286 230 : params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
1287 230 : params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
1288 230 : params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
1289 115 : params.set<bool>("use_displaced_mesh") = true;
1290 230 : params.set<Real>("c_t") = getParam<Real>("c_tangential");
1291 230 : params.set<Real>("c") = getParam<Real>("c_normal");
1292 230 : params.set<bool>("normalize_c") = getParam<bool>("normalize_c");
1293 115 : params.set<bool>("compute_primal_residuals") = false;
1294 :
1295 345 : params.set<MooseEnum>("segment_quadrature") = getParam<MooseEnum>("segment_quadrature");
1296 :
1297 345 : params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
1298 :
1299 115 : if (ndisp > 1)
1300 345 : params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
1301 115 : if (ndisp > 2)
1302 42 : params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
1303 :
1304 230 : params.set<NonlinearVariableName>("variable") = normal_lagrange_multiplier_name;
1305 345 : params.set<std::vector<VariableName>>("friction_lm") = {tangential_lagrange_multiplier_name};
1306 :
1307 115 : if (ndisp > 2)
1308 28 : params.set<std::vector<VariableName>>("friction_lm_dir") = {
1309 56 : tangential_lagrange_multiplier_3d_name};
1310 :
1311 230 : params.set<Real>("mu") = getParam<Real>("friction_coefficient");
1312 115 : params.applySpecificParameters(parameters(),
1313 : {"triangulation",
1314 : "triangulate_triangles",
1315 : "minimum_projection_angle",
1316 : "mortar_3d_subpatch_plane",
1317 : "mortar_3d_qp_mapping",
1318 : "extra_vector_tags",
1319 : "absolute_value_vector_tags",
1320 : "debug_mesh"});
1321 :
1322 115 : _problem->addConstraint(mortar_constraint_name, action_name + "_tangential_lm", params);
1323 115 : _problem->haveADObjects(true);
1324 115 : }
1325 :
1326 : const auto addMechanicalContactConstraints =
1327 556 : [this, &primary_subdomain_name, &secondary_subdomain_name, &displacements](
1328 : const std::string & variable_name,
1329 : const std::string & constraint_prefix,
1330 : const std::string & constraint_type,
1331 : const bool is_additional_frictional_constraint,
1332 : const bool is_normal_constraint)
1333 : {
1334 556 : InputParameters params = _factory.getValidParams(constraint_type);
1335 :
1336 1112 : params.set<bool>("correct_edge_dropping") = getParam<bool>("correct_edge_dropping");
1337 556 : params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
1338 556 : params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
1339 1112 : params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
1340 1112 : params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
1341 :
1342 556 : if (_formulation == ContactFormulation::MORTAR)
1343 950 : params.set<NonlinearVariableName>("variable") = variable_name;
1344 :
1345 1668 : params.set<MooseEnum>("segment_quadrature") = getParam<MooseEnum>("segment_quadrature");
1346 556 : params.set<bool>("use_displaced_mesh") = true;
1347 556 : params.set<bool>("compute_lm_residuals") = false;
1348 :
1349 : // Additional displacement residual for frictional problem
1350 : // The second frictional LM acts on a perpendicular direction.
1351 556 : if (is_additional_frictional_constraint)
1352 38 : params.set<MooseEnum>("direction") = "direction_2";
1353 556 : params.applySpecificParameters(parameters(),
1354 : {"triangulation",
1355 : "triangulate_triangles",
1356 : "minimum_projection_angle",
1357 : "mortar_3d_subpatch_plane",
1358 : "mortar_3d_qp_mapping",
1359 : "extra_vector_tags",
1360 : "absolute_value_vector_tags",
1361 : "debug_mesh"});
1362 :
1363 1804 : for (unsigned int i = 0; i < displacements.size(); ++i)
1364 : {
1365 1248 : std::string constraint_name = constraint_prefix + Moose::stringify(i);
1366 :
1367 1248 : params.set<VariableName>("secondary_variable") = displacements[i];
1368 1248 : params.set<MooseEnum>("component") = i;
1369 :
1370 1248 : if (is_normal_constraint && _model != ContactModel::COULOMB &&
1371 593 : _formulation == ContactFormulation::MORTAR)
1372 2120 : params.set<UserObjectName>("weighted_gap_uo") = "lm_weightedgap_object_" + name();
1373 362 : else if (is_normal_constraint && _model == ContactModel::COULOMB &&
1374 : _formulation == ContactFormulation::MORTAR)
1375 976 : params.set<UserObjectName>("weighted_gap_uo") = "lm_weightedvelocities_object_" + name();
1376 474 : else if (_formulation == ContactFormulation::MORTAR)
1377 572 : params.set<UserObjectName>("weighted_velocities_uo") =
1378 572 : "lm_weightedvelocities_object_" + name();
1379 188 : else if (is_normal_constraint && _model != ContactModel::COULOMB &&
1380 : _formulation == ContactFormulation::MORTAR_PENALTY)
1381 252 : params.set<UserObjectName>("weighted_gap_uo") = "penalty_weightedgap_object_" + name();
1382 55 : else if (is_normal_constraint && _model == ContactModel::COULOMB &&
1383 : _formulation == ContactFormulation::MORTAR_PENALTY)
1384 220 : params.set<UserObjectName>("weighted_gap_uo") = "penalty_friction_object_" + name();
1385 70 : else if (_formulation == ContactFormulation::MORTAR_PENALTY)
1386 140 : params.set<UserObjectName>("weighted_velocities_uo") =
1387 140 : "penalty_friction_object_" + name();
1388 :
1389 1248 : _problem->addConstraint(constraint_type, constraint_name, params);
1390 : }
1391 556 : _problem->haveADObjects(true);
1392 556 : };
1393 :
1394 : // Add mortar mechanical contact constraint objects for primal variables
1395 794 : addMechanicalContactConstraints(normal_lagrange_multiplier_name,
1396 794 : action_name + "_normal_constraint_",
1397 : "NormalMortarMechanicalContact",
1398 : /* is_additional_frictional_constraint = */ false,
1399 : /* is_normal_constraint = */ true);
1400 :
1401 397 : if (_model == ContactModel::COULOMB)
1402 : {
1403 280 : addMechanicalContactConstraints(tangential_lagrange_multiplier_name,
1404 280 : action_name + "_tangential_constraint_",
1405 : "TangentialMortarMechanicalContact",
1406 : /* is_additional_frictional_constraint = */ false,
1407 : /* is_normal_constraint = */ false);
1408 140 : if (ndisp > 2)
1409 57 : addMechanicalContactConstraints(tangential_lagrange_multiplier_3d_name,
1410 38 : action_name + "_tangential_constraint_3d_",
1411 : "TangentialMortarMechanicalContact",
1412 : /* is_additional_frictional_constraint = */ true,
1413 : /* is_normal_constraint = */ false);
1414 : }
1415 : }
1416 7230 : }
1417 :
1418 : void
1419 10502 : ContactAction::addNodeFaceContact()
1420 : {
1421 10502 : if (_current_task == "post_mesh_prepared" && _automatic_pairing_boundaries.size() > 0)
1422 : {
1423 56 : if (getParam<MooseEnum>("automatic_pairing_method").getEnum<ProximityMethod>() ==
1424 : ProximityMethod::NODE)
1425 21 : createSidesetsFromNodeProximity();
1426 14 : else if (getParam<MooseEnum>("automatic_pairing_method").getEnum<ProximityMethod>() ==
1427 : ProximityMethod::CENTROID)
1428 7 : createSidesetPairsFromGeometry();
1429 : }
1430 :
1431 10502 : if (_current_task != "add_constraint")
1432 9336 : return;
1433 :
1434 1166 : std::string action_name = MooseUtils::shortName(name());
1435 3498 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
1436 1166 : const unsigned int ndisp = displacements.size();
1437 :
1438 : std::string constraint_type;
1439 :
1440 1166 : if (_formulation == ContactFormulation::RANFS)
1441 : constraint_type = "RANFSNormalMechanicalContact";
1442 : else
1443 : constraint_type = "MechanicalContactConstraint";
1444 :
1445 1166 : InputParameters params = _factory.getValidParams(constraint_type);
1446 :
1447 1166 : params.applyParameters(parameters(),
1448 : {"displacements",
1449 : "secondary_gap_offset",
1450 : "mapped_primary_gap_offset",
1451 : "primary",
1452 : "secondary"});
1453 :
1454 1166 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
1455 1166 : .system()
1456 1166 : .variable_type(displacements[0])
1457 : .order.get_order();
1458 :
1459 1166 : params.set<std::vector<VariableName>>("displacements") = displacements;
1460 1166 : params.set<bool>("use_displaced_mesh") = true;
1461 2332 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
1462 :
1463 2439 : for (const auto & contact_pair : _boundary_pairs)
1464 : {
1465 1273 : if (_formulation != ContactFormulation::RANFS)
1466 : {
1467 3792 : params.set<std::vector<VariableName>>("nodal_area") = {"nodal_area"};
1468 1264 : params.set<BoundaryName>("boundary") = contact_pair.first;
1469 2528 : if (isParamValid("secondary_gap_offset"))
1470 12 : params.set<std::vector<VariableName>>("secondary_gap_offset") = {
1471 24 : getParam<VariableName>("secondary_gap_offset")};
1472 2528 : if (isParamValid("mapped_primary_gap_offset"))
1473 12 : params.set<std::vector<VariableName>>("mapped_primary_gap_offset") = {
1474 24 : getParam<VariableName>("mapped_primary_gap_offset")};
1475 : }
1476 :
1477 4082 : for (unsigned int i = 0; i < ndisp; ++i)
1478 : {
1479 8427 : std::string name = action_name + "_constraint_" + Moose::stringify(contact_pair, "_") + "_" +
1480 2809 : Moose::stringify(i);
1481 :
1482 2809 : if (_formulation == ContactFormulation::RANFS)
1483 36 : params.set<MooseEnum>("component") = i;
1484 : else
1485 2791 : params.set<unsigned int>("component") = i;
1486 :
1487 2809 : params.set<BoundaryName>("primary") = contact_pair.first;
1488 2809 : params.set<BoundaryName>("secondary") = contact_pair.second;
1489 5618 : params.set<NonlinearVariableName>("variable") = displacements[i];
1490 8427 : params.set<std::vector<VariableName>>("primary_variable") = {displacements[i]};
1491 2809 : params.applySpecificParameters(parameters(),
1492 : {"extra_vector_tags", "absolute_value_vector_tags"});
1493 2809 : _problem->addConstraint(constraint_type, name, params);
1494 : }
1495 : }
1496 2356 : }
1497 :
1498 : // Specialization for PointListAdaptor<MooseMesh::PeriodicNodeInfo>
1499 : // Return node location from NodeBoundaryIDInfo pairs
1500 : template <>
1501 : inline const Point &
1502 : PointListAdaptor<NodeBoundaryIDInfo>::getPoint(const NodeBoundaryIDInfo & item) const
1503 : {
1504 0 : return *(item.first);
1505 : }
1506 :
1507 : void
1508 21 : ContactAction::createSidesetsFromNodeProximity()
1509 : {
1510 21 : mooseInfo("The contact action is reading the list of boundaries and automatically pairs them "
1511 : "if the distance between nodes is less than a specified distance.");
1512 :
1513 21 : if (!_mesh)
1514 0 : mooseError("Failed to obtain mesh for automatically generating contact pairs.");
1515 :
1516 21 : if (!_mesh->getMesh().is_serial())
1517 0 : paramError(
1518 : "automatic_pairing_boundaries",
1519 : "The generation of automatic contact pairs in the contact action requires a serial mesh.");
1520 :
1521 : // Create automatic_pairing_boundaries_id
1522 : std::vector<BoundaryID> _automatic_pairing_boundaries_id;
1523 84 : for (const auto & sideset_name : _automatic_pairing_boundaries)
1524 63 : _automatic_pairing_boundaries_id.emplace_back(_mesh->getBoundaryID(sideset_name));
1525 :
1526 : // Vector of pairs node-boundary id
1527 : std::vector<NodeBoundaryIDInfo> node_boundary_id_vector;
1528 :
1529 : // Data structures to hold the boundary nodes
1530 21 : const ConstBndNodeRange & bnd_nodes = *_mesh->getBoundaryNodeRange();
1531 :
1532 2100 : for (const auto & bnode : bnd_nodes)
1533 : {
1534 2079 : const BoundaryID boundary_id = bnode->_bnd_id;
1535 2079 : const Node * node_ptr = bnode->_node;
1536 :
1537 : // Make sure node is on a boundary chosen for contact mechanics
1538 2079 : auto it = std::find(_automatic_pairing_boundaries_id.begin(),
1539 : _automatic_pairing_boundaries_id.end(),
1540 : boundary_id);
1541 :
1542 2079 : if (it != _automatic_pairing_boundaries_id.end())
1543 819 : node_boundary_id_vector.emplace_back(node_ptr, boundary_id);
1544 : }
1545 :
1546 : // sort by increasing boundary id
1547 21 : std::sort(node_boundary_id_vector.begin(),
1548 : node_boundary_id_vector.end(),
1549 : [](const NodeBoundaryIDInfo & first_pair, const NodeBoundaryIDInfo & second_pair)
1550 3171 : { return first_pair.second < second_pair.second; });
1551 :
1552 : // build kd-tree
1553 : using KDTreeType = nanoflann::KDTreeSingleIndexAdaptor<
1554 : nanoflann::L2_Simple_Adaptor<Real, PointListAdaptor<NodeBoundaryIDInfo>, Real, std::size_t>,
1555 : PointListAdaptor<NodeBoundaryIDInfo>,
1556 : LIBMESH_DIM,
1557 : std::size_t>;
1558 :
1559 : // This parameter can be tuned. Others use '10'
1560 : const unsigned int max_leaf_size = 20;
1561 :
1562 : // Build point list adaptor with all nodes-sidesets pairs for possible mechanical contact
1563 : auto point_list = PointListAdaptor<NodeBoundaryIDInfo>(node_boundary_id_vector.begin(),
1564 21 : node_boundary_id_vector.end());
1565 : auto kd_tree = std::make_unique<KDTreeType>(
1566 21 : LIBMESH_DIM, point_list, nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size));
1567 :
1568 21 : if (!kd_tree)
1569 0 : mooseError("Internal error. KDTree was not properly initialized in the contact action.");
1570 :
1571 21 : kd_tree->buildIndex();
1572 :
1573 : // data structures for kd-tree search
1574 : nanoflann::SearchParameters search_params;
1575 : std::vector<nanoflann::ResultItem<std::size_t, Real>> ret_matches;
1576 :
1577 42 : const auto radius_for_search = getParam<Real>("automatic_pairing_distance");
1578 :
1579 : // For all nodes
1580 840 : for (const auto & pair : node_boundary_id_vector)
1581 : {
1582 : // clear result buffer
1583 819 : ret_matches.clear();
1584 :
1585 : // position where we expect a periodic partner for the current node and boundary
1586 819 : const Point search_point = *pair.first;
1587 :
1588 : // search at the expected point
1589 819 : kd_tree->radiusSearch(
1590 819 : &(search_point)(0), radius_for_search * radius_for_search, ret_matches, search_params);
1591 :
1592 12180 : for (auto & match_pair : ret_matches)
1593 : {
1594 11361 : const auto & match = node_boundary_id_vector[match_pair.first];
1595 :
1596 : //
1597 : // If the proximity node identified belongs to a boundary in the input, add boundary pair
1598 : //
1599 :
1600 : // Make sure node is on a boundary chosen for contact mechanics
1601 11361 : auto it = std::find(_automatic_pairing_boundaries_id.begin(),
1602 : _automatic_pairing_boundaries_id.end(),
1603 11361 : match.second);
1604 :
1605 : // If nodes are on the same boundary, pass.
1606 11361 : if (match.second == pair.second)
1607 : continue;
1608 :
1609 : // At this point we will likely create many repeated pairs because many nodal pairs may
1610 : // fulfill the distance condition imposed by the automatic_pairing_distance user input
1611 : // parameter.
1612 7266 : if (it != _automatic_pairing_boundaries_id.end())
1613 : {
1614 : const auto index_one = cast_int<int>(it - _automatic_pairing_boundaries_id.begin());
1615 7266 : auto it_other = std::find(_automatic_pairing_boundaries_id.begin(),
1616 : _automatic_pairing_boundaries_id.end(),
1617 7266 : pair.second);
1618 :
1619 : mooseAssert(it_other != _automatic_pairing_boundaries_id.end(),
1620 : "Error in contact action. Unable to find boundary ID for node proximity "
1621 : "automatic pairing.");
1622 :
1623 : const auto index_two = cast_int<int>(it_other - _automatic_pairing_boundaries_id.begin());
1624 :
1625 7266 : if (pair.second > match.second)
1626 3633 : _boundary_pairs.push_back(
1627 3633 : {_automatic_pairing_boundaries[index_two], _automatic_pairing_boundaries[index_one]});
1628 : else
1629 3633 : _boundary_pairs.push_back(
1630 3633 : {_automatic_pairing_boundaries[index_one], _automatic_pairing_boundaries[index_two]});
1631 : }
1632 : }
1633 : }
1634 :
1635 : // Let's remove likely repeated pairs
1636 21 : removeRepeatedPairs();
1637 :
1638 21 : mooseInfo(
1639 : "The following boundary pairs were created by the contact action using nodal proximity: ");
1640 63 : for (const auto & [primary, secondary] : _boundary_pairs)
1641 : mooseInfoRepeated(
1642 : "Primary boundary ID: ", primary, " and secondary boundary ID: ", secondary, ".");
1643 21 : }
1644 :
1645 : void
1646 7 : ContactAction::createSidesetPairsFromGeometry()
1647 : {
1648 7 : mooseInfo("The contact action is reading the list of boundaries and automatically pairs them "
1649 : "if their centroids fall within a specified distance of each other.");
1650 :
1651 7 : if (!_mesh)
1652 0 : mooseError("Failed to obtain mesh for automatically generating contact pairs.");
1653 :
1654 7 : if (!_mesh->getMesh().is_serial())
1655 0 : paramError(
1656 : "automatic_pairing_boundaries",
1657 : "The generation of automatic contact pairs in the contact action requires a serial mesh.");
1658 :
1659 : // Compute centers of gravity for each sideset
1660 : std::vector<std::pair<BoundaryName, Point>> automatic_pairing_boundaries_cog;
1661 7 : const auto & sideset_ids = _mesh->meshSidesetIds();
1662 :
1663 7 : const auto & bnd_to_elem_map = _mesh->getBoundariesToActiveSemiLocalElemIds();
1664 :
1665 28 : for (const auto & sideset_name : _automatic_pairing_boundaries)
1666 : {
1667 : // If the sideset provided in the input file isn't in the mesh, error out.
1668 21 : const auto find_set = sideset_ids.find(_mesh->getBoundaryID(sideset_name));
1669 21 : if (find_set == sideset_ids.end())
1670 0 : paramError("automatic_pairing_boundaries",
1671 : sideset_name,
1672 : " is not defined as a sideset in the mesh.");
1673 :
1674 21 : auto dofs_set = bnd_to_elem_map.find(_mesh->getBoundaryID(sideset_name));
1675 :
1676 : // Initialize data for sideset
1677 : Point center_of_gravity(0, 0, 0);
1678 : Real accumulated_sideset_area(0);
1679 :
1680 : // Pointer to lower-dimensional element on the sideset
1681 21 : std::unique_ptr<const Elem> side_ptr;
1682 : const std::unordered_set<dof_id_type> & bnd_elems = dofs_set->second;
1683 :
1684 273 : for (auto elem_id : bnd_elems)
1685 : {
1686 252 : const Elem * elem = _mesh->elemPtr(elem_id);
1687 252 : unsigned int side = _mesh->sideWithBoundaryID(elem, _mesh->getBoundaryID(sideset_name));
1688 :
1689 : // update side_ptr
1690 252 : elem->side_ptr(side_ptr, side);
1691 :
1692 : // area of the (linearized) side
1693 252 : const auto side_area = side_ptr->volume();
1694 :
1695 : // position of the side
1696 252 : const auto side_position = side_ptr->true_centroid();
1697 :
1698 : center_of_gravity += side_position * side_area;
1699 252 : accumulated_sideset_area += side_area;
1700 : }
1701 :
1702 : // Average each element's center of gravity (centroid) with its area
1703 : center_of_gravity /= accumulated_sideset_area;
1704 :
1705 : // Add sideset-cog pair to vector
1706 21 : automatic_pairing_boundaries_cog.emplace_back(sideset_name, center_of_gravity);
1707 21 : }
1708 :
1709 : // Vectors of distances for each pair
1710 : std::vector<std::pair<std::pair<BoundaryName, BoundaryName>, Real>> pairs_distances;
1711 :
1712 : // Assign distances to identify nearby pairs.
1713 21 : for (std::size_t i = 0; i < automatic_pairing_boundaries_cog.size() - 1; i++)
1714 35 : for (std::size_t j = i + 1; j < automatic_pairing_boundaries_cog.size(); j++)
1715 : {
1716 : const Point & distance_vector =
1717 : automatic_pairing_boundaries_cog[i].second - automatic_pairing_boundaries_cog[j].second;
1718 :
1719 21 : if (automatic_pairing_boundaries_cog[i].first != automatic_pairing_boundaries_cog[j].first)
1720 : {
1721 21 : const Real distance = distance_vector.norm();
1722 21 : const std::pair pair = std::make_pair(automatic_pairing_boundaries_cog[i].first,
1723 21 : automatic_pairing_boundaries_cog[j].first);
1724 42 : pairs_distances.emplace_back(std::make_pair(pair, distance));
1725 : }
1726 : }
1727 :
1728 14 : const auto automatic_pairing_distance = getParam<Real>("automatic_pairing_distance");
1729 :
1730 : // Loop over all pairs
1731 : std::vector<std::pair<std::pair<BoundaryName, BoundaryName>, Real>> lean_pairs_distances;
1732 28 : for (const auto & pair_distance : pairs_distances)
1733 21 : if (pair_distance.second <= automatic_pairing_distance)
1734 : {
1735 21 : lean_pairs_distances.emplace_back(pair_distance);
1736 : mooseInfoRepeated("Generating contact pair primary--secondary ",
1737 21 : pair_distance.first.first,
1738 : "--",
1739 21 : pair_distance.first.second,
1740 : ", with a relative distance of ",
1741 21 : pair_distance.second);
1742 : }
1743 :
1744 : // Create the boundary pairs (possibly with repeated pairs depending on user input)
1745 28 : for (const auto & lean_pairs_distance : lean_pairs_distances)
1746 : {
1747 : // Make sure secondary surface's boundary ID is less than primary surface's boundary ID.
1748 : // This is done to ensure some consistency in the boundary matching, which helps in defining
1749 : // auxiliary kernels in the input file.
1750 42 : if (_mesh->getBoundaryID(lean_pairs_distance.first.first) >
1751 21 : _mesh->getBoundaryID(lean_pairs_distance.first.second))
1752 0 : _boundary_pairs.push_back(
1753 : {lean_pairs_distance.first.first, lean_pairs_distance.first.second});
1754 : else
1755 42 : _boundary_pairs.push_back(
1756 : {lean_pairs_distance.first.second, lean_pairs_distance.first.first});
1757 : }
1758 :
1759 : // Let's remove possibly repeated pairs
1760 7 : removeRepeatedPairs();
1761 7 : }
1762 :
1763 : MooseEnum
1764 10867 : ContactAction::getModelEnum()
1765 : {
1766 21734 : return MooseEnum(getContactModelOptions(), "frictionless");
1767 : }
1768 :
1769 : MooseEnum
1770 1593 : ContactAction::getProximityMethod()
1771 : {
1772 3186 : return MooseEnum(getProximityMethodOptions());
1773 : }
1774 :
1775 : MooseEnum
1776 9274 : ContactAction::getFormulationEnum()
1777 : {
1778 18548 : auto formulations = MooseEnum(getContactFormulationOptions(), "kinematic");
1779 :
1780 18548 : formulations.addDocumentation(
1781 : "ranfs",
1782 : "Reduced Active Nonlinear Function Set scheme for node-on-face contact. Provides exact "
1783 : "enforcement without Lagrange multipliers or penalty terms.");
1784 18548 : formulations.addDocumentation(
1785 : "kinematic",
1786 : "Kinematic contact constraint enforcement transfers the internal forces at secondary nodes "
1787 : "to the corresponding primary face for node-on-face contact. Provides exact "
1788 : "enforcement without Lagrange multipliers or penalty terms.");
1789 18548 : formulations.addDocumentation(
1790 : "penalty",
1791 : "Node-on-face penalty based contact constraint enforcement. Interpenetration is penalized. "
1792 : "Enforcement depends on the penalty magnitude. High penalties can introduce ill conditioning "
1793 : "of the system.");
1794 18548 : formulations.addDocumentation("augmented_lagrange",
1795 : "Node-on-face augmented Lagrange penalty based contact constraint "
1796 : "enforcement. Interpenetration is enforced up to a user specified "
1797 : "tolerance, ill-conditioning is generally avoided. Requires an "
1798 : "Augmented Lagrange Problem class to be used in the simulation.");
1799 18548 : formulations.addDocumentation(
1800 : "tangential_penalty",
1801 : "Node-on-face penalty based frictional contact constraint enforcement. Interpenetration and "
1802 : "slip distance for sticking nodes are penalized. Enforcement depends on the penalty "
1803 : "magnitudes. High penalties can introduce ill conditioning of the system.");
1804 18548 : formulations.addDocumentation(
1805 : "mortar",
1806 : "Mortar based contact constraint enforcement using Lagrange multipliers. Provides exact "
1807 : "enforcement and a variationally consistent formulation. Lagrange multipliers introduce a "
1808 : "saddle point character in the system matrix which can have a negative impact on scalability "
1809 : "with iterative solvers");
1810 18548 : formulations.addDocumentation(
1811 : "mortar_penalty",
1812 : "Mortar and penalty based contact constraint enforcement. When using an Augmented Lagrange "
1813 : "Problem class this provides normal (and tangential) contact constratint enforced up to a "
1814 : "user specified tolerances. Without AL the enforcement depends on the penalty magnitudes. "
1815 : "High penalties can introduce ill conditioning of the system.");
1816 :
1817 9274 : return formulations;
1818 0 : }
1819 :
1820 : MooseEnum
1821 0 : ContactAction::getSystemEnum()
1822 : {
1823 0 : return MooseEnum("Constraint", "Constraint");
1824 : }
1825 :
1826 : MooseEnum
1827 9274 : ContactAction::getSmoothingEnum()
1828 : {
1829 18548 : return MooseEnum("edge_based nodal_normal_based", "");
1830 : }
1831 :
1832 : InputParameters
1833 9274 : ContactAction::commonParameters()
1834 : {
1835 9274 : InputParameters params = emptyInputParameters();
1836 :
1837 18548 : params.addParam<MooseEnum>("normal_smoothing_method",
1838 18548 : ContactAction::getSmoothingEnum(),
1839 : "Method to use to smooth normals");
1840 18548 : params.addParam<Real>(
1841 : "normal_smoothing_distance",
1842 : "Distance from edge in parametric coordinates over which to smooth contact normal");
1843 :
1844 18548 : params.addParam<MooseEnum>(
1845 18548 : "formulation", ContactAction::getFormulationEnum(), "The contact formulation");
1846 :
1847 18548 : params.addParam<MooseEnum>("model", ContactAction::getModelEnum(), "The contact model to use");
1848 :
1849 9274 : return params;
1850 0 : }
|
