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