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