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 2372 : ContactAction::validParams()
73 : {
74 2372 : InputParameters params = Action::validParams();
75 2372 : params += ContactAction::commonParameters();
76 :
77 4744 : params.addParam<std::vector<BoundaryName>>(
78 : "primary", "The list of boundary IDs referring to primary sidesets");
79 4744 : params.addParam<std::vector<BoundaryName>>(
80 : "secondary", "The list of boundary IDs referring to secondary sidesets");
81 4744 : 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 4744 : 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 4744 : 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 4744 : params.addParam<VariableName>("secondary_gap_offset",
102 : "Offset to gap distance from secondary side");
103 4744 : params.addParam<VariableName>("mapped_primary_gap_offset",
104 : "Offset to gap distance mapped from primary side");
105 4744 : params.addParam<std::vector<VariableName>>(
106 : "displacements",
107 : {},
108 : "The displacements appropriate for the simulation geometry and coordinate system");
109 4744 : params.addParam<Real>(
110 : "penalty",
111 4744 : 1e8,
112 : "The penalty to apply. This can vary depending on the stiffness of your materials");
113 4744 : params.addParam<Real>(
114 : "penalty_friction",
115 4744 : 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 7116 : params.addRangeCheckedParam<Real>(
119 : "penalty_multiplier",
120 4744 : 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 7116 : params.addRangeCheckedParam<Real>(
126 : "penalty_multiplier_friction",
127 4744 : 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 4744 : params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
133 4744 : params.addParam<Real>("tension_release",
134 4744 : 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 4744 : params.addParam<MooseEnum>("model", ContactAction::getModelEnum(), "The contact model to use");
139 4744 : params.addParam<Real>("tangential_tolerance",
140 : "Tangential distance to extend edges of contact surfaces");
141 4744 : params.addParam<Real>("capture_tolerance",
142 4744 : 0.0,
143 : "Normal distance from surface within which nodes are captured. This "
144 : "parameter is used for node-face and mortar formulations.");
145 4744 : params.addParam<Real>(
146 : "normal_smoothing_distance",
147 : "Distance from edge in parametric coordinates over which to smooth contact normal");
148 :
149 4744 : params.addParam<bool>("normalize_penalty",
150 4744 : false,
151 : "Whether to normalize the penalty parameter with the nodal area.");
152 4744 : params.addParam<bool>(
153 : "primary_secondary_jacobian",
154 4744 : true,
155 : "Whether to include Jacobian entries coupling primary and secondary nodes.");
156 4744 : params.addParam<Real>("al_penetration_tolerance",
157 : "The tolerance of the penetration for augmented Lagrangian method.");
158 4744 : params.addParam<Real>("al_incremental_slip_tolerance",
159 : "The tolerance of the incremental slip for augmented Lagrangian method.");
160 7116 : params.addRangeCheckedParam<Real>(
161 : "max_penalty_multiplier",
162 4744 : 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 4744 : MooseEnum adaptivity_penalty_normal("SIMPLE BUSSETTA", "SIMPLE");
169 4744 : adaptivity_penalty_normal.addDocumentation(
170 : "SIMPLE", "Keep multiplying by the penalty multiplier between AL iterations");
171 4744 : 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 4744 : params.addParam<MooseEnum>(
176 : "adaptivity_penalty_normal",
177 : adaptivity_penalty_normal,
178 : "The augmented Lagrange update strategy used on the normal penalty coefficient.");
179 4744 : MooseEnum adaptivity_penalty_friction("SIMPLE FRICTION_LIMIT", "FRICTION_LIMIT");
180 4744 : adaptivity_penalty_friction.addDocumentation(
181 : "SIMPLE", "Keep multiplying by the frictional penalty multiplier between AL iterations");
182 4744 : 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 4744 : params.addParam<MooseEnum>(
187 : "adaptivity_penalty_friction",
188 : adaptivity_penalty_friction,
189 : "The augmented Lagrange update strategy used on the frictional penalty coefficient.");
190 4744 : params.addParam<Real>("al_frictional_force_tolerance",
191 : "The tolerance of the frictional force for augmented Lagrangian method.");
192 4744 : params.addParam<Real>(
193 : "c_normal",
194 4744 : 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 4744 : params.addParam<Real>(
201 4744 : "c_tangential", 1, "Numerical parameter for nonlinear mortar frictional constraints");
202 4744 : params.addParam<bool>("ping_pong_protection",
203 4744 : 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 4744 : params.addParam<Real>(
209 : "normal_lm_scaling",
210 4744 : 1.,
211 : "Scaling factor to apply to the normal LM variable for a mortar formulation");
212 4744 : params.addParam<Real>(
213 : "tangential_lm_scaling",
214 4744 : 1.,
215 : "Scaling factor to apply to the tangential LM variable for a mortar formulation");
216 4744 : params.addParam<bool>(
217 : "normalize_c",
218 4744 : 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 2372 : params.addClassDescription("Sets up all objects needed for mechanical contact enforcement");
224 4744 : 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 4744 : params.addParam<bool>(
233 : "correct_edge_dropping",
234 4744 : 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 4744 : params.addParam<bool>(
239 : "generate_mortar_mesh",
240 4744 : 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 4744 : params.addParam<MooseEnum>("automatic_pairing_method",
244 4744 : ContactAction::getProximityMethod(),
245 : "The proximity method used for automatic pairing of boundaries.");
246 4744 : params.addParam<bool>(
247 : "mortar_dynamics",
248 4744 : 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 4744 : params.addParam<Real>(
253 : "newmark_beta",
254 4744 : 0.25,
255 : "Newmark-beta beta parameter for its inclusion in the weighted gap update formula");
256 4744 : params.addParam<Real>(
257 : "newmark_gamma",
258 4744 : 0.5,
259 : "Newmark-beta gamma parameter for its inclusion in the weighted gap update formula");
260 4744 : params.addCoupledVar("wear_depth",
261 : "The name of the mortar auxiliary variable that is used to modify the "
262 : "weighted gap definition");
263 4744 : params.addParam<std::vector<TagName>>(
264 : "extra_vector_tags",
265 : "The tag names for extra vectors that residual data should be saved into");
266 4744 : 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 4744 : params.addParam<bool>(
271 : "use_petrov_galerkin",
272 4744 : 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 2372 : return params;
277 2372 : }
278 :
279 2372 : ContactAction::ContactAction(const InputParameters & params)
280 : : Action(params),
281 4744 : _boundary_pairs(getParam<BoundaryName, BoundaryName>("primary", "secondary")),
282 4744 : _model(getParam<MooseEnum>("model").getEnum<ContactModel>()),
283 4744 : _formulation(getParam<MooseEnum>("formulation").getEnum<ContactFormulation>()),
284 4744 : _generate_mortar_mesh(getParam<bool>("generate_mortar_mesh")),
285 9488 : _mortar_dynamics(getParam<bool>("mortar_dynamics"))
286 : {
287 : // Check for automatic selection of contact pairs.
288 4744 : if (getParam<std::vector<BoundaryName>>("automatic_pairing_boundaries").size() > 1)
289 : _automatic_pairing_boundaries =
290 120 : getParam<std::vector<BoundaryName>>("automatic_pairing_boundaries");
291 :
292 2492 : 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 2492 : 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 2372 : 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 2372 : 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 2370 : 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 2370 : 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 2370 : if (_formulation == ContactFormulation::MORTAR_PENALTY)
322 : {
323 : // Use dual basis functions for contact traction interpolation
324 148 : if (isParamValid("use_dual"))
325 52 : _use_dual = getParam<bool>("use_dual");
326 : else
327 48 : _use_dual = true;
328 :
329 74 : if (_model == ContactModel::GLUED)
330 0 : paramError("model", "The 'mortar_penalty' formulation does not support glued contact");
331 :
332 148 : if (getParam<bool>("mortar_dynamics"))
333 0 : paramError("mortar_dynamics",
334 : "The 'mortar_penalty' formulation does not support implicit dynamic simulations");
335 :
336 148 : 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 2370 : if (_formulation == ContactFormulation::MORTAR)
345 : {
346 405 : 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 810 : if (isParamValid("use_dual"))
351 76 : _use_dual = getParam<bool>("use_dual");
352 : else
353 367 : _use_dual = true;
354 :
355 810 : if (!getParam<bool>("mortar_dynamics"))
356 : {
357 762 : if (params.isParamSetByUser("newmark_beta"))
358 2 : paramError("newmark_beta", "newmark_beta can only be used with the mortar_dynamics option");
359 :
360 758 : 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 3930 : 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 3930 : else if (params.isParamSetByUser("use_dual") &&
373 26 : _formulation != ContactFormulation::MORTAR_PENALTY)
374 0 : paramError("use_dual",
375 : "The 'use_dual' option can only be used with the 'mortar' formulation");
376 3930 : 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 3930 : 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 3930 : 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 2366 : if (_formulation == ContactFormulation::RANFS)
389 : {
390 36 : 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 36 : 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 4696 : 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 2366 : removeRepeatedPairs();
405 2366 : }
406 :
407 : void
408 2406 : ContactAction::removeRepeatedPairs()
409 : {
410 2406 : 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 15322 : for (const auto & [primary, secondary] : _boundary_pairs)
419 : {
420 : // Structured bindings are not capturable (primary_copy, secondary_copy)
421 12916 : auto it = std::find_if(lean_boundary_pairs.begin(),
422 : lean_boundary_pairs.end(),
423 25832 : [&, primary_copy = primary, secondary_copy = secondary](
424 : const std::pair<BoundaryName, BoundaryName> & lean_pair)
425 : {
426 21000 : const bool match_one = lean_pair.second == secondary_copy &&
427 13862 : lean_pair.first == primary_copy;
428 21000 : const bool match_two = lean_pair.second == primary_copy &&
429 18 : lean_pair.first == secondary_copy;
430 21000 : const bool exist = match_one || match_two;
431 21000 : return exist;
432 : });
433 :
434 12916 : if (it == lean_boundary_pairs.end())
435 2560 : lean_boundary_pairs.emplace_back(primary, secondary);
436 : else
437 10356 : 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 2406 : _boundary_pairs = lean_boundary_pairs;
446 2406 : }
447 :
448 : void
449 21280 : ContactAction::act()
450 : {
451 : // proform problem checks/corrections once during the first feasible task
452 21280 : if (_current_task == "add_contact_aux_variable")
453 : {
454 4732 : 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 2364 : if (!_problem->isSNESMFReuseBaseSetbyUser())
464 : _problem->setSNESMFReuseBase(false, false);
465 : }
466 :
467 21278 : if (_formulation == ContactFormulation::MORTAR ||
468 : _formulation == ContactFormulation::MORTAR_PENALTY)
469 4269 : addMortarContact();
470 : else
471 17009 : addNodeFaceContact();
472 :
473 21276 : if (_current_task == "add_aux_kernel")
474 : { // Add ContactPenetrationAuxAction.
475 4724 : 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 4918 : for (const auto & contact_pair : _boundary_pairs)
480 : {
481 : {
482 2556 : InputParameters params = _factory.getValidParams("PenetrationAux");
483 10224 : params.applyParameters(parameters(),
484 : {"secondary_gap_offset", "mapped_primary_gap_offset", "order"});
485 :
486 : std::vector<VariableName> displacements =
487 2556 : getParam<std::vector<VariableName>>("displacements");
488 2556 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
489 2556 : .system()
490 2556 : .variable_type(displacements[0])
491 : .order.get_order();
492 :
493 5112 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
494 10224 : params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_LINEAR};
495 7668 : params.set<std::vector<BoundaryName>>("boundary") = {contact_pair.second};
496 5112 : params.set<BoundaryName>("paired_boundary") = contact_pair.first;
497 5112 : params.set<AuxVariableName>("variable") = "penetration";
498 5112 : if (isParamValid("secondary_gap_offset"))
499 22 : params.set<std::vector<VariableName>>("secondary_gap_offset") = {
500 44 : getParam<VariableName>("secondary_gap_offset")};
501 5112 : if (isParamValid("mapped_primary_gap_offset"))
502 22 : params.set<std::vector<VariableName>>("mapped_primary_gap_offset") = {
503 44 : getParam<VariableName>("mapped_primary_gap_offset")};
504 2556 : params.set<bool>("use_displaced_mesh") = true;
505 5112 : std::string name = _name + "_contact_" + Moose::stringify(contact_auxkernel_counter++);
506 :
507 5112 : _problem->addAuxKernel("PenetrationAux", name, params);
508 2556 : }
509 : }
510 :
511 2362 : addContactPressureAuxKernel();
512 :
513 4724 : const unsigned int ndisp = getParam<std::vector<VariableName>>("displacements").size();
514 :
515 : // Add MortarFrictionalPressureVectorAux
516 2362 : if (_formulation == ContactFormulation::MORTAR && _model == ContactModel::COULOMB && ndisp > 2)
517 : {
518 : {
519 32 : InputParameters params = _factory.getValidParams("MortarFrictionalPressureVectorAux");
520 :
521 32 : params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
522 32 : params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
523 48 : params.set<std::vector<BoundaryName>>("boundary") = {_boundary_pairs[0].second};
524 48 : params.set<ExecFlagEnum>("execute_on", true) = {EXEC_NONLINEAR};
525 :
526 16 : std::string action_name = MooseUtils::shortName(name());
527 16 : const std::string tangential_lagrange_multiplier_name = action_name + "_tangential_lm";
528 : const std::string tangential_lagrange_multiplier_3d_name =
529 16 : action_name + "_tangential_3d_lm";
530 :
531 32 : params.set<std::vector<VariableName>>("tangent_one") = {
532 64 : tangential_lagrange_multiplier_name};
533 32 : params.set<std::vector<VariableName>>("tangent_two") = {
534 64 : tangential_lagrange_multiplier_3d_name};
535 :
536 80 : std::vector<std::string> disp_components({"x", "y", "z"});
537 : unsigned component_index = 0;
538 :
539 : // Loop over three displacements
540 64 : for (const auto & disp_component : disp_components)
541 : {
542 144 : params.set<AuxVariableName>("variable") = _name + "_tangent_" + disp_component;
543 48 : params.set<unsigned int>("component") = component_index;
544 :
545 96 : std::string name = _name + "_mortar_frictional_pressure_" + disp_component + "_" +
546 96 : Moose::stringify(contact_mortar_auxkernel_counter++);
547 :
548 48 : _problem->addAuxKernel("MortarFrictionalPressureVectorAux", name, params);
549 48 : component_index++;
550 : }
551 32 : }
552 : }
553 : }
554 :
555 21276 : if (_current_task == "add_contact_aux_variable")
556 : {
557 4728 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
558 2364 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
559 2364 : .system()
560 2364 : .variable_type(displacements[0])
561 : .order.get_order();
562 : // Add ContactPenetrationVarAction
563 : {
564 2364 : auto var_params = _factory.getValidParams("MooseVariable");
565 4728 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
566 4728 : var_params.set<MooseEnum>("family") = "LAGRANGE";
567 :
568 4728 : _problem->addAuxVariable("MooseVariable", "penetration", var_params);
569 2364 : }
570 : // Add ContactPressureVarAction
571 : {
572 2364 : auto var_params = _factory.getValidParams("MooseVariable");
573 4728 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
574 4728 : var_params.set<MooseEnum>("family") = "LAGRANGE";
575 :
576 4728 : _problem->addAuxVariable("MooseVariable", "contact_pressure", var_params);
577 2364 : }
578 : // Add nodal area contact variable
579 : {
580 2364 : auto var_params = _factory.getValidParams("MooseVariable");
581 4728 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
582 4728 : var_params.set<MooseEnum>("family") = "LAGRANGE";
583 :
584 4728 : _problem->addAuxVariable("MooseVariable", "nodal_area", var_params);
585 2364 : }
586 :
587 4728 : const unsigned int ndisp = getParam<std::vector<VariableName>>("displacements").size();
588 :
589 : // Add MortarFrictionalPressureVectorAux variables
590 2364 : if (_formulation == ContactFormulation::MORTAR && _model == ContactModel::COULOMB && ndisp > 2)
591 : {
592 : {
593 80 : std::vector<std::string> disp_components({"x", "y", "z"});
594 : // Loop over three displacements
595 64 : for (const auto & disp_component : disp_components)
596 : {
597 48 : auto var_params = _factory.getValidParams("MooseVariable");
598 96 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
599 96 : var_params.set<MooseEnum>("family") = "LAGRANGE";
600 :
601 144 : _problem->addAuxVariable(
602 48 : "MooseVariable", _name + "_tangent_" + disp_component, var_params);
603 48 : }
604 16 : }
605 : }
606 2364 : }
607 :
608 21276 : if (_current_task == "add_user_object")
609 : {
610 2362 : auto var_params = _factory.getValidParams("NodalArea");
611 :
612 : // Get secondary_boundary_vector from possibly updated set from the
613 : // ContactAction constructor cleanup
614 2362 : const auto actions = _awh.getActions<ContactAction>();
615 :
616 : std::vector<BoundaryName> secondary_boundary_vector;
617 5012 : for (const auto * const action : actions)
618 5494 : for (const auto j : index_range(action->_boundary_pairs))
619 2844 : secondary_boundary_vector.push_back(action->_boundary_pairs[j].second);
620 :
621 4724 : var_params.set<std::vector<BoundaryName>>("boundary") = secondary_boundary_vector;
622 7086 : var_params.set<std::vector<VariableName>>("variable") = {"nodal_area"};
623 :
624 : mooseAssert(_problem, "Problem pointer is NULL");
625 9448 : var_params.set<ExecFlagEnum>("execute_on", true) = {EXEC_INITIAL, EXEC_TIMESTEP_BEGIN};
626 2362 : var_params.set<bool>("use_displaced_mesh") = true;
627 :
628 7086 : _problem->addUserObject("NodalArea",
629 2362 : "nodal_area_object_" + Moose::stringify(contact_userobject_counter++),
630 : var_params);
631 4724 : }
632 31430 : }
633 :
634 : void
635 2362 : ContactAction::addContactPressureAuxKernel()
636 : {
637 : // Add ContactPressureAux: Only one object for all contact pairs
638 : // if (_formulation != ContactFormulation::MORTAR)
639 2362 : const auto actions = _awh.getActions<ContactAction>();
640 :
641 : // Increment counter for contact action objects
642 2362 : contact_action_counter++;
643 : // Add auxiliary kernel if we are the last contact action object.
644 2362 : if (contact_action_counter == actions.size())
645 : {
646 : std::vector<BoundaryName> boundary_vector;
647 : std::vector<BoundaryName> pair_boundary_vector;
648 :
649 4616 : for (const auto * const action : actions)
650 4918 : for (const auto j : index_range(action->_boundary_pairs))
651 : {
652 2556 : boundary_vector.push_back(action->_boundary_pairs[j].second);
653 2556 : pair_boundary_vector.push_back(action->_boundary_pairs[j].first);
654 : }
655 :
656 2254 : InputParameters params = _factory.getValidParams("ContactPressureAux");
657 4508 : params.applyParameters(parameters(), {"order"});
658 :
659 2254 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
660 2254 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
661 2254 : .system()
662 2254 : .variable_type(displacements[0])
663 : .order.get_order();
664 :
665 4508 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
666 2254 : params.set<std::vector<BoundaryName>>("boundary") = boundary_vector;
667 4508 : params.set<std::vector<BoundaryName>>("paired_boundary") = pair_boundary_vector;
668 4508 : params.set<AuxVariableName>("variable") = "contact_pressure";
669 4508 : params.addRequiredCoupledVar("nodal_area", "The nodal area");
670 6762 : params.set<std::vector<VariableName>>("nodal_area") = {"nodal_area"};
671 2254 : params.set<bool>("use_displaced_mesh") = true;
672 :
673 2254 : std::string name = _name + "_contact_pressure";
674 4508 : params.set<ExecFlagEnum>("execute_on",
675 13524 : true) = {EXEC_NONLINEAR, EXEC_TIMESTEP_END, EXEC_TIMESTEP_BEGIN};
676 4508 : _problem->addAuxKernel("ContactPressureAux", name, params);
677 2254 : }
678 9124 : }
679 :
680 : void
681 7090 : ContactAction::addRelationshipManagers(Moose::RelationshipManagerType input_rm_type)
682 : {
683 7090 : if (_formulation == ContactFormulation::MORTAR ||
684 : _formulation == ContactFormulation::MORTAR_PENALTY)
685 : {
686 1421 : auto params = MortarConstraintBase::validParams();
687 1421 : params.set<bool>("use_displaced_mesh") = true;
688 1421 : std::string action_name = MooseUtils::shortName(name());
689 1421 : const std::string primary_subdomain_name = action_name + "_primary_subdomain";
690 1421 : const std::string secondary_subdomain_name = action_name + "_secondary_subdomain";
691 2842 : params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
692 2842 : params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
693 2842 : params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
694 2842 : params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
695 2842 : params.set<bool>("use_petrov_galerkin") = getParam<bool>("use_petrov_galerkin");
696 1421 : addRelationshipManagers(input_rm_type, params);
697 1421 : }
698 7090 : }
699 :
700 : void
701 4269 : ContactAction::addMortarContact()
702 : {
703 4269 : std::string action_name = MooseUtils::shortName(name());
704 :
705 12807 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
706 4269 : const unsigned int ndisp = displacements.size();
707 :
708 : // Definitions for mortar contact.
709 4269 : const std::string primary_subdomain_name = action_name + "_primary_subdomain";
710 4269 : const std::string secondary_subdomain_name = action_name + "_secondary_subdomain";
711 4269 : const std::string normal_lagrange_multiplier_name = action_name + "_normal_lm";
712 4269 : const std::string tangential_lagrange_multiplier_name = action_name + "_tangential_lm";
713 4269 : const std::string tangential_lagrange_multiplier_3d_name = action_name + "_tangential_3d_lm";
714 4269 : const std::string auxiliary_lagrange_multiplier_name = action_name + "_aux_lm";
715 :
716 4269 : 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 475 : if (!(_app.isRecovering() && _app.isUltimateMaster()) && !_app.masterMesh() &&
722 369 : _generate_mortar_mesh)
723 : {
724 361 : const MeshGeneratorName primary_name = primary_subdomain_name + "_generator";
725 361 : const MeshGeneratorName secondary_name = secondary_subdomain_name + "_generator";
726 :
727 361 : auto primary_params = _factory.getValidParams("LowerDBlockFromSidesetGenerator");
728 722 : auto secondary_params = _factory.getValidParams("LowerDBlockFromSidesetGenerator");
729 :
730 722 : primary_params.set<SubdomainName>("new_block_name") = primary_subdomain_name;
731 722 : secondary_params.set<SubdomainName>("new_block_name") = secondary_subdomain_name;
732 :
733 1083 : primary_params.set<std::vector<BoundaryName>>("sidesets") = {_boundary_pairs[0].first};
734 1083 : secondary_params.set<std::vector<BoundaryName>>("sidesets") = {_boundary_pairs[0].second};
735 :
736 722 : _app.appendMeshGenerator("LowerDBlockFromSidesetGenerator", primary_name, primary_params);
737 722 : _app.appendMeshGenerator("LowerDBlockFromSidesetGenerator", secondary_name, secondary_params);
738 361 : }
739 : }
740 :
741 : // Add the lagrange multiplier on the secondary subdomain.
742 : const auto addLagrangeMultiplier =
743 679 : [this, &secondary_subdomain_name, &displacements](const std::string & variable_name,
744 : const Real scaling_factor,
745 : const bool add_aux_lm,
746 679 : const bool penalty_traction) //
747 : {
748 679 : 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 679 : if (!add_aux_lm || penalty_traction)
754 663 : 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 1358 : _problem->systemBaseNonlinear(/*nl_sys_num=*/0).system().variable_type(displacements[0]);
760 :
761 : const int lm_order = primal_type.order.get_order();
762 :
763 679 : if (primal_type.family == LAGRANGE)
764 : {
765 1358 : params.set<MooseEnum>("family") = Utility::enum_to_string<FEFamily>(primal_type.family);
766 1358 : 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 2037 : params.set<std::vector<SubdomainName>>("block") = {secondary_subdomain_name};
772 679 : if (!(add_aux_lm || penalty_traction))
773 1230 : params.set<std::vector<Real>>("scaling") = {scaling_factor};
774 :
775 679 : auto fe_type = AddVariableAction::feType(params);
776 679 : auto var_type = AddVariableAction::variableType(fe_type);
777 679 : if (add_aux_lm || penalty_traction)
778 64 : _problem->addAuxVariable(var_type, variable_name, params);
779 : else
780 615 : _problem->addVariable(var_type, variable_name, params);
781 679 : };
782 :
783 4269 : if (_current_task == "add_mortar_variable" && _formulation == ContactFormulation::MORTAR)
784 : {
785 401 : addLagrangeMultiplier(
786 401 : normal_lagrange_multiplier_name, getParam<Real>("normal_lm_scaling"), false, false);
787 :
788 401 : if (_model == ContactModel::COULOMB)
789 : {
790 198 : addLagrangeMultiplier(tangential_lagrange_multiplier_name,
791 198 : getParam<Real>("tangential_lm_scaling"),
792 : false,
793 : false);
794 198 : if (ndisp > 2)
795 16 : addLagrangeMultiplier(tangential_lagrange_multiplier_3d_name,
796 32 : getParam<Real>("tangential_lm_scaling"),
797 : false,
798 : false);
799 : }
800 :
801 802 : if (getParam<bool>("use_petrov_galerkin"))
802 16 : addLagrangeMultiplier(auxiliary_lagrange_multiplier_name, 1.0, true, false);
803 : }
804 3868 : else if (_current_task == "add_mortar_variable" &&
805 74 : _formulation == ContactFormulation::MORTAR_PENALTY)
806 : {
807 74 : if (_use_dual)
808 48 : addLagrangeMultiplier(auxiliary_lagrange_multiplier_name, 1.0, false, true);
809 : }
810 :
811 4269 : if (_current_task == "add_user_object")
812 : {
813 : // check if the correct problem class is selected if AL parameters are provided
814 475 : if (_formulation == ContactFormulation::MORTAR_PENALTY &&
815 74 : !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 228 : "al_frictional_force_tolerance"};
822 218 : for (const auto & param : params)
823 182 : 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 36 : }
828 :
829 473 : if (_model != ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR)
830 : {
831 406 : auto var_params = _factory.getValidParams("LMWeightedGapUserObject");
832 :
833 406 : var_params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
834 406 : var_params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
835 406 : var_params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
836 406 : var_params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
837 609 : var_params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
838 406 : var_params.set<bool>("correct_edge_dropping") = getParam<bool>("correct_edge_dropping");
839 609 : var_params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
840 203 : if (ndisp > 2)
841 48 : var_params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
842 203 : var_params.set<bool>("use_displaced_mesh") = true;
843 609 : var_params.set<std::vector<VariableName>>("lm_variable") = {normal_lagrange_multiplier_name};
844 406 : var_params.set<bool>("use_petrov_galerkin") = getParam<bool>("use_petrov_galerkin");
845 406 : if (getParam<bool>("use_petrov_galerkin"))
846 24 : var_params.set<std::vector<VariableName>>("aux_lm") = {auxiliary_lagrange_multiplier_name};
847 :
848 609 : _problem->addUserObject(
849 203 : "LMWeightedGapUserObject", "lm_weightedgap_object_" + name(), var_params);
850 203 : }
851 270 : else if (_model == ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR)
852 : {
853 396 : auto var_params = _factory.getValidParams("LMWeightedVelocitiesUserObject");
854 396 : var_params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
855 396 : var_params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
856 396 : var_params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
857 396 : var_params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
858 594 : var_params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
859 396 : var_params.set<bool>("correct_edge_dropping") = getParam<bool>("correct_edge_dropping");
860 594 : var_params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
861 198 : if (ndisp > 2)
862 48 : var_params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
863 :
864 198 : var_params.set<VariableName>("secondary_variable") = displacements[0];
865 198 : var_params.set<bool>("use_displaced_mesh") = true;
866 396 : var_params.set<std::vector<VariableName>>("lm_variable_normal") = {
867 792 : normal_lagrange_multiplier_name};
868 396 : var_params.set<std::vector<VariableName>>("lm_variable_tangential_one") = {
869 792 : tangential_lagrange_multiplier_name};
870 198 : if (ndisp > 2)
871 32 : var_params.set<std::vector<VariableName>>("lm_variable_tangential_two") = {
872 64 : tangential_lagrange_multiplier_3d_name};
873 396 : var_params.set<bool>("use_petrov_galerkin") = getParam<bool>("use_petrov_galerkin");
874 396 : if (getParam<bool>("use_petrov_galerkin"))
875 24 : var_params.set<std::vector<VariableName>>("aux_lm") = {auxiliary_lagrange_multiplier_name};
876 :
877 594 : _problem->addUserObject(
878 198 : "LMWeightedVelocitiesUserObject", "lm_weightedvelocities_object_" + name(), var_params);
879 198 : }
880 :
881 473 : if (_model != ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR_PENALTY)
882 : {
883 56 : auto var_params = _factory.getValidParams("PenaltyWeightedGapUserObject");
884 :
885 56 : var_params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
886 56 : var_params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
887 56 : var_params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
888 56 : var_params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
889 84 : var_params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
890 56 : var_params.set<bool>("correct_edge_dropping") = getParam<bool>("correct_edge_dropping");
891 84 : var_params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
892 56 : var_params.set<Real>("penalty") = getParam<Real>("penalty");
893 :
894 : // AL parameters
895 56 : var_params.set<Real>("max_penalty_multiplier") = getParam<Real>("max_penalty_multiplier");
896 56 : var_params.set<MooseEnum>("adaptivity_penalty_normal") =
897 56 : getParam<MooseEnum>("adaptivity_penalty_normal");
898 56 : if (isParamValid("al_penetration_tolerance"))
899 0 : var_params.set<Real>("penetration_tolerance") = getParam<Real>("al_penetration_tolerance");
900 56 : if (isParamValid("penalty_multiplier"))
901 56 : 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 28 : var_params.set<bool>("use_physical_gap") = true;
905 :
906 28 : if (_use_dual)
907 30 : var_params.set<std::vector<VariableName>>("aux_lm") = {auxiliary_lagrange_multiplier_name};
908 :
909 28 : if (ndisp > 2)
910 0 : var_params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
911 28 : var_params.set<bool>("use_displaced_mesh") = true;
912 :
913 84 : _problem->addUserObject(
914 28 : "PenaltyWeightedGapUserObject", "penalty_weightedgap_object_" + name(), var_params);
915 28 : _problem->haveADObjects(true);
916 28 : }
917 445 : else if (_model == ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR_PENALTY)
918 : {
919 88 : auto var_params = _factory.getValidParams("PenaltyFrictionUserObject");
920 88 : var_params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
921 88 : var_params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
922 88 : var_params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
923 88 : var_params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
924 132 : var_params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
925 88 : var_params.set<bool>("correct_edge_dropping") = getParam<bool>("correct_edge_dropping");
926 132 : var_params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
927 44 : if (ndisp > 2)
928 0 : var_params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
929 :
930 44 : var_params.set<VariableName>("secondary_variable") = displacements[0];
931 44 : var_params.set<bool>("use_displaced_mesh") = true;
932 88 : var_params.set<Real>("friction_coefficient") = getParam<Real>("friction_coefficient");
933 88 : var_params.set<Real>("penalty") = getParam<Real>("penalty");
934 88 : var_params.set<Real>("penalty_friction") = getParam<Real>("penalty_friction");
935 :
936 : // AL parameters
937 88 : var_params.set<Real>("max_penalty_multiplier") = getParam<Real>("max_penalty_multiplier");
938 88 : var_params.set<MooseEnum>("adaptivity_penalty_normal") =
939 88 : getParam<MooseEnum>("adaptivity_penalty_normal");
940 88 : var_params.set<MooseEnum>("adaptivity_penalty_friction") =
941 88 : getParam<MooseEnum>("adaptivity_penalty_friction");
942 88 : if (isParamValid("al_penetration_tolerance"))
943 72 : var_params.set<Real>("penetration_tolerance") = getParam<Real>("al_penetration_tolerance");
944 88 : if (isParamValid("penalty_multiplier"))
945 88 : var_params.set<Real>("penalty_multiplier") = getParam<Real>("penalty_multiplier");
946 88 : if (isParamValid("penalty_multiplier_friction"))
947 44 : var_params.set<Real>("penalty_multiplier_friction") =
948 132 : getParam<Real>("penalty_multiplier_friction");
949 :
950 88 : if (isParamValid("al_incremental_slip_tolerance"))
951 72 : 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 44 : var_params.set<bool>("use_physical_gap") = true;
955 :
956 44 : if (_use_dual)
957 108 : var_params.set<std::vector<VariableName>>("aux_lm") = {auxiliary_lagrange_multiplier_name};
958 :
959 176 : var_params.applySpecificParameters(parameters(),
960 : {"friction_coefficient", "penalty", "penalty_friction"});
961 :
962 88 : _problem->addUserObject(
963 44 : "PenaltyFrictionUserObject", "penalty_friction_object_" + name(), var_params);
964 44 : _problem->haveADObjects(true);
965 44 : }
966 : }
967 :
968 4267 : if (_current_task == "add_constraint")
969 : {
970 : // Prepare problem for enforcement with Lagrange multipliers
971 473 : if (_model != ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR)
972 : {
973 : std::string mortar_constraint_name;
974 :
975 203 : if (!_mortar_dynamics)
976 : mortar_constraint_name = "ComputeWeightedGapLMMechanicalContact";
977 : else
978 : mortar_constraint_name = "ComputeDynamicWeightedGapLMMechanicalContact";
979 :
980 203 : InputParameters params = _factory.getValidParams(mortar_constraint_name);
981 203 : if (_mortar_dynamics)
982 40 : params.applySpecificParameters(
983 : parameters(), {"newmark_beta", "newmark_gamma", "capture_tolerance", "wear_depth"});
984 :
985 : else // We need user objects for quasistatic constraints
986 780 : params.set<UserObjectName>("weighted_gap_uo") = "lm_weightedgap_object_" + name();
987 :
988 406 : params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
989 406 : params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
990 406 : params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
991 406 : params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
992 406 : params.set<NonlinearVariableName>("variable") = normal_lagrange_multiplier_name;
993 609 : params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
994 406 : params.set<Real>("c") = getParam<Real>("c_normal");
995 :
996 203 : if (ndisp > 1)
997 609 : params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
998 203 : if (ndisp > 2)
999 48 : params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
1000 :
1001 203 : params.set<bool>("use_displaced_mesh") = true;
1002 :
1003 1015 : params.applySpecificParameters(parameters(),
1004 : {"correct_edge_dropping",
1005 : "normalize_c",
1006 : "extra_vector_tags",
1007 : "absolute_value_vector_tags"});
1008 :
1009 203 : _problem->addConstraint(
1010 203 : mortar_constraint_name, action_name + "_normal_lm_weighted_gap", params);
1011 203 : _problem->haveADObjects(true);
1012 406 : }
1013 : // Add the tangential and normal Lagrange's multiplier constraints on the secondary boundary.
1014 270 : else if (_model == ContactModel::COULOMB && _formulation == ContactFormulation::MORTAR)
1015 : {
1016 : std::string mortar_constraint_name;
1017 :
1018 198 : if (!_mortar_dynamics)
1019 : mortar_constraint_name = "ComputeFrictionalForceLMMechanicalContact";
1020 : else
1021 : mortar_constraint_name = "ComputeDynamicFrictionalForceLMMechanicalContact";
1022 :
1023 198 : InputParameters params = _factory.getValidParams(mortar_constraint_name);
1024 198 : if (_mortar_dynamics)
1025 80 : params.applySpecificParameters(
1026 : parameters(), {"newmark_beta", "newmark_gamma", "capture_tolerance", "wear_depth"});
1027 : else
1028 : { // We need user objects for quasistatic constraints
1029 546 : params.set<UserObjectName>("weighted_gap_uo") = "lm_weightedvelocities_object_" + name();
1030 364 : params.set<UserObjectName>("weighted_velocities_uo") =
1031 364 : "lm_weightedvelocities_object_" + name();
1032 : }
1033 :
1034 396 : params.set<bool>("correct_edge_dropping") = getParam<bool>("correct_edge_dropping");
1035 396 : params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
1036 396 : params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
1037 396 : params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
1038 396 : params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
1039 198 : params.set<bool>("use_displaced_mesh") = true;
1040 396 : params.set<Real>("c_t") = getParam<Real>("c_tangential");
1041 396 : params.set<Real>("c") = getParam<Real>("c_normal");
1042 396 : params.set<bool>("normalize_c") = getParam<bool>("normalize_c");
1043 198 : params.set<bool>("compute_primal_residuals") = false;
1044 :
1045 594 : params.set<std::vector<VariableName>>("disp_x") = {displacements[0]};
1046 :
1047 198 : if (ndisp > 1)
1048 594 : params.set<std::vector<VariableName>>("disp_y") = {displacements[1]};
1049 198 : if (ndisp > 2)
1050 48 : params.set<std::vector<VariableName>>("disp_z") = {displacements[2]};
1051 :
1052 396 : params.set<NonlinearVariableName>("variable") = normal_lagrange_multiplier_name;
1053 594 : params.set<std::vector<VariableName>>("friction_lm") = {tangential_lagrange_multiplier_name};
1054 :
1055 198 : if (ndisp > 2)
1056 32 : params.set<std::vector<VariableName>>("friction_lm_dir") = {
1057 64 : tangential_lagrange_multiplier_3d_name};
1058 :
1059 396 : params.set<Real>("mu") = getParam<Real>("friction_coefficient");
1060 594 : params.applySpecificParameters(parameters(),
1061 : {"extra_vector_tags", "absolute_value_vector_tags"});
1062 :
1063 198 : _problem->addConstraint(mortar_constraint_name, action_name + "_tangential_lm", params);
1064 198 : _problem->haveADObjects(true);
1065 198 : }
1066 :
1067 : const auto addMechanicalContactConstraints =
1068 731 : [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 731 : const bool is_normal_constraint)
1074 : {
1075 731 : InputParameters params = _factory.getValidParams(constraint_type);
1076 :
1077 1462 : params.set<bool>("correct_edge_dropping") = getParam<bool>("correct_edge_dropping");
1078 1462 : params.set<BoundaryName>("primary_boundary") = _boundary_pairs[0].first;
1079 1462 : params.set<BoundaryName>("secondary_boundary") = _boundary_pairs[0].second;
1080 1462 : params.set<SubdomainName>("primary_subdomain") = primary_subdomain_name;
1081 1462 : params.set<SubdomainName>("secondary_subdomain") = secondary_subdomain_name;
1082 :
1083 731 : if (_formulation == ContactFormulation::MORTAR)
1084 1230 : params.set<NonlinearVariableName>("variable") = variable_name;
1085 :
1086 731 : params.set<bool>("use_displaced_mesh") = true;
1087 731 : 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 731 : if (is_additional_frictional_constraint)
1092 32 : params.set<MooseEnum>("direction") = "direction_2";
1093 2193 : params.applySpecificParameters(parameters(),
1094 : {"extra_vector_tags", "absolute_value_vector_tags"});
1095 :
1096 2257 : for (unsigned int i = 0; i < displacements.size(); ++i)
1097 : {
1098 3052 : std::string constraint_name = constraint_prefix + Moose::stringify(i);
1099 :
1100 1526 : params.set<VariableName>("secondary_variable") = displacements[i];
1101 1526 : params.set<MooseEnum>("component") = i;
1102 :
1103 1526 : if (is_normal_constraint && _model != ContactModel::COULOMB &&
1104 478 : _formulation == ContactFormulation::MORTAR)
1105 1688 : params.set<UserObjectName>("weighted_gap_uo") = "lm_weightedgap_object_" + name();
1106 556 : else if (is_normal_constraint && _model == ContactModel::COULOMB &&
1107 : _formulation == ContactFormulation::MORTAR)
1108 1648 : params.set<UserObjectName>("weighted_gap_uo") = "lm_weightedvelocities_object_" + name();
1109 692 : else if (_formulation == ContactFormulation::MORTAR)
1110 920 : params.set<UserObjectName>("weighted_velocities_uo") =
1111 920 : "lm_weightedvelocities_object_" + name();
1112 232 : else if (is_normal_constraint && _model != ContactModel::COULOMB &&
1113 : _formulation == ContactFormulation::MORTAR_PENALTY)
1114 224 : params.set<UserObjectName>("weighted_gap_uo") = "penalty_weightedgap_object_" + name();
1115 88 : else if (is_normal_constraint && _model == ContactModel::COULOMB &&
1116 : _formulation == ContactFormulation::MORTAR_PENALTY)
1117 352 : params.set<UserObjectName>("weighted_gap_uo") = "penalty_friction_object_" + name();
1118 88 : else if (_formulation == ContactFormulation::MORTAR_PENALTY)
1119 176 : params.set<UserObjectName>("weighted_velocities_uo") =
1120 176 : "penalty_friction_object_" + name();
1121 :
1122 1526 : _problem->addConstraint(constraint_type, constraint_name, params);
1123 : }
1124 731 : _problem->haveADObjects(true);
1125 2193 : };
1126 :
1127 : // Add mortar mechanical contact constraint objects for primal variables
1128 946 : addMechanicalContactConstraints(normal_lagrange_multiplier_name,
1129 946 : action_name + "_normal_constraint_",
1130 : "NormalMortarMechanicalContact",
1131 : /* is_additional_frictional_constraint = */ false,
1132 : /* is_normal_constraint = */ true);
1133 :
1134 473 : if (_model == ContactModel::COULOMB)
1135 : {
1136 484 : addMechanicalContactConstraints(tangential_lagrange_multiplier_name,
1137 484 : action_name + "_tangential_constraint_",
1138 : "TangentialMortarMechanicalContact",
1139 : /* is_additional_frictional_constraint = */ false,
1140 : /* is_normal_constraint = */ false);
1141 242 : if (ndisp > 2)
1142 48 : addMechanicalContactConstraints(tangential_lagrange_multiplier_3d_name,
1143 32 : action_name + "_tangential_constraint_3d_",
1144 : "TangentialMortarMechanicalContact",
1145 : /* is_additional_frictional_constraint = */ true,
1146 : /* is_normal_constraint = */ false);
1147 : }
1148 : }
1149 9548 : }
1150 :
1151 : void
1152 17009 : ContactAction::addNodeFaceContact()
1153 : {
1154 17009 : if (_current_task == "post_mesh_prepared" && _automatic_pairing_boundaries.size() > 0)
1155 : {
1156 80 : if (getParam<MooseEnum>("automatic_pairing_method").getEnum<ProximityMethod>() ==
1157 : ProximityMethod::NODE)
1158 30 : createSidesetsFromNodeProximity();
1159 20 : else if (getParam<MooseEnum>("automatic_pairing_method").getEnum<ProximityMethod>() ==
1160 : ProximityMethod::CENTROID)
1161 10 : createSidesetPairsFromGeometry();
1162 : }
1163 :
1164 17009 : if (_current_task != "add_constraint")
1165 15120 : return;
1166 :
1167 1889 : std::string action_name = MooseUtils::shortName(name());
1168 5667 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
1169 1889 : const unsigned int ndisp = displacements.size();
1170 :
1171 : std::string constraint_type;
1172 :
1173 1889 : if (_formulation == ContactFormulation::RANFS)
1174 : constraint_type = "RANFSNormalMechanicalContact";
1175 : else
1176 : constraint_type = "MechanicalContactConstraint";
1177 :
1178 1889 : InputParameters params = _factory.getValidParams(constraint_type);
1179 :
1180 11334 : params.applyParameters(parameters(),
1181 : {"displacements",
1182 : "secondary_gap_offset",
1183 : "mapped_primary_gap_offset",
1184 : "primary",
1185 : "secondary"});
1186 :
1187 1889 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
1188 1889 : .system()
1189 1889 : .variable_type(displacements[0])
1190 : .order.get_order();
1191 :
1192 1889 : params.set<std::vector<VariableName>>("displacements") = displacements;
1193 1889 : params.set<bool>("use_displaced_mesh") = true;
1194 3778 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
1195 :
1196 3972 : for (const auto & contact_pair : _boundary_pairs)
1197 : {
1198 2083 : if (_formulation != ContactFormulation::RANFS)
1199 : {
1200 6195 : params.set<std::vector<VariableName>>("nodal_area") = {"nodal_area"};
1201 2065 : params.set<BoundaryName>("boundary") = contact_pair.first;
1202 4130 : if (isParamValid("secondary_gap_offset"))
1203 22 : params.set<std::vector<VariableName>>("secondary_gap_offset") = {
1204 44 : getParam<VariableName>("secondary_gap_offset")};
1205 4130 : if (isParamValid("mapped_primary_gap_offset"))
1206 22 : params.set<std::vector<VariableName>>("mapped_primary_gap_offset") = {
1207 44 : getParam<VariableName>("mapped_primary_gap_offset")};
1208 : }
1209 :
1210 6681 : for (unsigned int i = 0; i < ndisp; ++i)
1211 : {
1212 13794 : std::string name = action_name + "_constraint_" + Moose::stringify(contact_pair, "_") + "_" +
1213 4598 : Moose::stringify(i);
1214 :
1215 4598 : if (_formulation == ContactFormulation::RANFS)
1216 72 : params.set<MooseEnum>("component") = i;
1217 : else
1218 4562 : params.set<unsigned int>("component") = i;
1219 :
1220 4598 : params.set<BoundaryName>("primary") = contact_pair.first;
1221 4598 : params.set<BoundaryName>("secondary") = contact_pair.second;
1222 9196 : params.set<NonlinearVariableName>("variable") = displacements[i];
1223 13794 : params.set<std::vector<VariableName>>("primary_variable") = {displacements[i]};
1224 13794 : params.applySpecificParameters(parameters(),
1225 : {"extra_vector_tags", "absolute_value_vector_tags"});
1226 4598 : _problem->addConstraint(constraint_type, name, params);
1227 : }
1228 : }
1229 16796 : }
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 30 : ContactAction::createSidesetsFromNodeProximity()
1242 : {
1243 30 : 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 30 : if (!_mesh)
1247 0 : mooseError("Failed to obtain mesh for automatically generating contact pairs.");
1248 :
1249 30 : 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 120 : for (const auto & sideset_name : _automatic_pairing_boundaries)
1257 90 : _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 30 : const ConstBndNodeRange & bnd_nodes = *_mesh->getBoundaryNodeRange();
1264 :
1265 3000 : for (const auto & bnode : bnd_nodes)
1266 : {
1267 2970 : const BoundaryID boundary_id = bnode->_bnd_id;
1268 2970 : const Node * node_ptr = bnode->_node;
1269 :
1270 : // Make sure node is on a boundary chosen for contact mechanics
1271 2970 : auto it = std::find(_automatic_pairing_boundaries_id.begin(),
1272 : _automatic_pairing_boundaries_id.end(),
1273 : boundary_id);
1274 :
1275 2970 : if (it != _automatic_pairing_boundaries_id.end())
1276 1170 : node_boundary_id_vector.emplace_back(node_ptr, boundary_id);
1277 : }
1278 :
1279 : // sort by increasing boundary id
1280 30 : std::sort(node_boundary_id_vector.begin(),
1281 : node_boundary_id_vector.end(),
1282 : [](const NodeBoundaryIDInfo & first_pair, const NodeBoundaryIDInfo & second_pair)
1283 4530 : { 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 : node_boundary_id_vector.end());
1298 : auto kd_tree = std::make_unique<KDTreeType>(
1299 30 : LIBMESH_DIM, point_list, nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size));
1300 :
1301 30 : if (!kd_tree)
1302 0 : mooseError("Internal error. KDTree was not properly initialized in the contact action.");
1303 :
1304 30 : 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 60 : const auto radius_for_search = getParam<Real>("automatic_pairing_distance");
1311 :
1312 : // For all nodes
1313 1200 : for (const auto & pair : node_boundary_id_vector)
1314 : {
1315 : // clear result buffer
1316 : ret_matches.clear();
1317 :
1318 : // position where we expect a periodic partner for the current node and boundary
1319 1170 : const Point search_point = *pair.first;
1320 :
1321 : // search at the expected point
1322 1170 : kd_tree->radiusSearch(
1323 1170 : &(search_point)(0), radius_for_search * radius_for_search, ret_matches, search_params);
1324 :
1325 17400 : for (auto & match_pair : ret_matches)
1326 : {
1327 16230 : 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 16230 : auto it = std::find(_automatic_pairing_boundaries_id.begin(),
1335 : _automatic_pairing_boundaries_id.end(),
1336 16230 : match.second);
1337 :
1338 : // If nodes are on the same boundary, pass.
1339 16230 : 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 10380 : if (it != _automatic_pairing_boundaries_id.end())
1346 : {
1347 : const auto index_one = cast_int<int>(it - _automatic_pairing_boundaries_id.begin());
1348 10380 : auto it_other = std::find(_automatic_pairing_boundaries_id.begin(),
1349 : _automatic_pairing_boundaries_id.end(),
1350 10380 : 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 10380 : if (pair.second > match.second)
1359 5190 : _boundary_pairs.push_back(
1360 5190 : {_automatic_pairing_boundaries[index_two], _automatic_pairing_boundaries[index_one]});
1361 : else
1362 5190 : _boundary_pairs.push_back(
1363 5190 : {_automatic_pairing_boundaries[index_one], _automatic_pairing_boundaries[index_two]});
1364 : }
1365 : }
1366 : }
1367 :
1368 : // Let's remove likely repeated pairs
1369 30 : removeRepeatedPairs();
1370 :
1371 30 : mooseInfo(
1372 : "The following boundary pairs were created by the contact action using nodal proximity: ");
1373 90 : for (const auto & [primary, secondary] : _boundary_pairs)
1374 : mooseInfoRepeated(
1375 : "Primary boundary ID: ", primary, " and secondary boundary ID: ", secondary, ".");
1376 60 : }
1377 :
1378 : void
1379 10 : ContactAction::createSidesetPairsFromGeometry()
1380 : {
1381 10 : 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 10 : if (!_mesh)
1385 0 : mooseError("Failed to obtain mesh for automatically generating contact pairs.");
1386 :
1387 10 : 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 10 : const auto & sideset_ids = _mesh->meshSidesetIds();
1395 :
1396 10 : const auto & bnd_to_elem_map = _mesh->getBoundariesToActiveSemiLocalElemIds();
1397 :
1398 40 : 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 30 : const auto find_set = sideset_ids.find(_mesh->getBoundaryID(sideset_name));
1402 30 : 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 30 : 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 30 : std::unique_ptr<const Elem> side_ptr;
1415 : const std::unordered_set<dof_id_type> & bnd_elems = dofs_set->second;
1416 :
1417 390 : for (auto elem_id : bnd_elems)
1418 : {
1419 360 : const Elem * elem = _mesh->elemPtr(elem_id);
1420 360 : unsigned int side = _mesh->sideWithBoundaryID(elem, _mesh->getBoundaryID(sideset_name));
1421 :
1422 : // update side_ptr
1423 360 : elem->side_ptr(side_ptr, side);
1424 :
1425 : // area of the (linearized) side
1426 360 : const auto side_area = side_ptr->volume();
1427 :
1428 : // position of the side
1429 360 : const auto side_position = side_ptr->true_centroid();
1430 :
1431 : center_of_gravity += side_position * side_area;
1432 360 : 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 30 : automatic_pairing_boundaries_cog.emplace_back(sideset_name, center_of_gravity);
1440 30 : }
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 30 : for (std::size_t i = 0; i < automatic_pairing_boundaries_cog.size() - 1; i++)
1447 50 : 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 30 : if (automatic_pairing_boundaries_cog[i].first != automatic_pairing_boundaries_cog[j].first)
1453 : {
1454 30 : const Real distance = distance_vector.norm();
1455 30 : const std::pair pair = std::make_pair(automatic_pairing_boundaries_cog[i].first,
1456 30 : automatic_pairing_boundaries_cog[j].first);
1457 60 : pairs_distances.emplace_back(std::make_pair(pair, distance));
1458 : }
1459 : }
1460 :
1461 20 : 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 40 : for (const auto & pair_distance : pairs_distances)
1466 30 : if (pair_distance.second <= automatic_pairing_distance)
1467 : {
1468 30 : lean_pairs_distances.emplace_back(pair_distance);
1469 : mooseInfoRepeated("Generating contact pair primary--secondary ",
1470 30 : pair_distance.first.first,
1471 : "--",
1472 30 : pair_distance.first.second,
1473 : ", with a relative distance of ",
1474 30 : pair_distance.second);
1475 : }
1476 :
1477 : // Create the boundary pairs (possibly with repeated pairs depending on user input)
1478 40 : 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 60 : if (_mesh->getBoundaryID(lean_pairs_distance.first.first) >
1484 30 : _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 60 : _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 10 : removeRepeatedPairs();
1494 10 : }
1495 :
1496 : MooseEnum
1497 11177 : ContactAction::getModelEnum()
1498 : {
1499 22354 : return MooseEnum("frictionless glued coulomb", "frictionless");
1500 : }
1501 :
1502 : MooseEnum
1503 2372 : ContactAction::getProximityMethod()
1504 : {
1505 4744 : return MooseEnum("node centroid");
1506 : }
1507 :
1508 : MooseEnum
1509 8805 : ContactAction::getFormulationEnum()
1510 : {
1511 : auto formulations = MooseEnum(
1512 : "ranfs kinematic penalty augmented_lagrange tangential_penalty mortar mortar_penalty",
1513 17610 : "kinematic");
1514 :
1515 17610 : 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 17610 : 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 17610 : 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 17610 : 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 17610 : 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 17610 : 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 17610 : 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 8805 : return formulations;
1553 0 : }
1554 :
1555 : MooseEnum
1556 0 : ContactAction::getSystemEnum()
1557 : {
1558 0 : return MooseEnum("Constraint", "Constraint");
1559 : }
1560 :
1561 : MooseEnum
1562 8805 : ContactAction::getSmoothingEnum()
1563 : {
1564 17610 : return MooseEnum("edge_based nodal_normal_based", "");
1565 : }
1566 :
1567 : InputParameters
1568 8805 : ContactAction::commonParameters()
1569 : {
1570 8805 : InputParameters params = emptyInputParameters();
1571 :
1572 17610 : params.addParam<MooseEnum>("normal_smoothing_method",
1573 17610 : ContactAction::getSmoothingEnum(),
1574 : "Method to use to smooth normals");
1575 17610 : params.addParam<Real>(
1576 : "normal_smoothing_distance",
1577 : "Distance from edge in parametric coordinates over which to smooth contact normal");
1578 :
1579 17610 : params.addParam<MooseEnum>(
1580 17610 : "formulation", ContactAction::getFormulationEnum(), "The contact formulation");
1581 :
1582 17610 : params.addParam<MooseEnum>("model", ContactAction::getModelEnum(), "The contact model to use");
1583 :
1584 8805 : return params;
1585 0 : }
|
