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