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