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