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