26 #include <unordered_map> 29 #include "libmesh/petsc_nonlinear_solver.h" 30 #include "libmesh/string_to_enum.h" 33 #if NANOFLANN_VERSION < 0x150 38 template <
typename T,
typename U>
77 params.
addParam<std::vector<BoundaryName>>(
78 "primary",
"The list of boundary IDs referring to primary sidesets");
79 params.
addParam<std::vector<BoundaryName>>(
80 "secondary",
"The list of boundary IDs referring to secondary sidesets");
81 params.
addParam<std::vector<BoundaryName>>(
82 "automatic_pairing_boundaries",
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. ");
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.");
99 "The mesh generator for mortar method",
100 "This parameter is not used anymore and can simply be removed");
101 params.
addParam<VariableName>(
"secondary_gap_offset",
102 "Offset to gap distance from secondary side");
103 params.
addParam<VariableName>(
"mapped_primary_gap_offset",
104 "Offset to gap distance mapped from primary side");
105 params.
addParam<std::vector<VariableName>>(
108 "The displacements appropriate for the simulation geometry and coordinate system");
112 "The penalty to apply. This can vary depending on the stiffness of your materials");
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");
119 "penalty_multiplier",
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 " 126 "penalty_multiplier_friction",
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 params.
addParam<
Real>(
"friction_coefficient", 0,
"The friction coefficient");
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.");
140 "Tangential distance to extend edges of contact surfaces");
143 "Normal distance from surface within which nodes are captured. This " 144 "parameter is used for node-face and mortar formulations.");
146 "normal_smoothing_distance",
147 "Distance from edge in parametric coordinates over which to smooth contact normal");
149 params.
addParam<
bool>(
"normalize_penalty",
151 "Whether to normalize the penalty parameter with the nodal area.");
153 "primary_secondary_jacobian",
155 "Whether to include Jacobian entries coupling primary and secondary nodes.");
157 "The tolerance of the penetration for augmented Lagrangian method.");
159 "The tolerance of the incremental slip for augmented Lagrangian method.");
161 "max_penalty_multiplier",
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 MooseEnum adaptivity_penalty_normal(
"SIMPLE BUSSETTA",
"SIMPLE");
169 adaptivity_penalty_normal.addDocumentation(
170 "SIMPLE",
"Keep multiplying by the penalty multiplier between AL iterations");
171 adaptivity_penalty_normal.addDocumentation(
173 "Modify the penalty using an algorithm from Bussetta et al, 2012, Comput Mech 49:259-275 " 174 "between AL iterations.");
176 "adaptivity_penalty_normal",
177 adaptivity_penalty_normal,
178 "The augmented Lagrange update strategy used on the normal penalty coefficient.");
179 MooseEnum adaptivity_penalty_friction(
"SIMPLE FRICTION_LIMIT",
"FRICTION_LIMIT");
180 adaptivity_penalty_friction.addDocumentation(
181 "SIMPLE",
"Keep multiplying by the frictional penalty multiplier between AL iterations");
182 adaptivity_penalty_friction.addDocumentation(
184 "This strategy will be guided by the Coulomb limit and be less reliant on the initial " 185 "penalty factor provided by the user.");
187 "adaptivity_penalty_friction",
188 adaptivity_penalty_friction,
189 "The augmented Lagrange update strategy used on the frictional penalty coefficient.");
191 "The tolerance of the frictional force for augmented Lagrangian method.");
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.");
201 "c_tangential", 1,
"Numerical parameter for nonlinear mortar frictional constraints");
202 params.
addParam<
bool>(
"ping_pong_protection",
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");
211 "Scaling factor to apply to the normal LM variable for a mortar formulation");
213 "tangential_lm_scaling",
215 "Scaling factor to apply to the tangential LM variable for a mortar formulation");
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 params.
addClassDescription(
"Sets up all objects needed for mechanical contact enforcement");
226 "Whether to use the dual mortar approach within a mortar formulation. It is defaulted to " 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.");
233 "correct_edge_dropping",
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.");
239 "generate_mortar_mesh",
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.");
245 "The proximity method used for automatic pairing of boundaries.");
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");
255 "Newmark-beta beta parameter for its inclusion in the weighted gap update formula");
259 "Newmark-beta gamma parameter for its inclusion in the weighted gap update formula");
261 "The name of the mortar auxiliary variable that is used to modify the " 262 "weighted gap definition");
263 params.
addParam<std::vector<TagName>>(
265 "The tag names for extra vectors that residual data should be saved into");
266 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");
271 "use_petrov_galerkin",
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.");
281 _boundary_pairs(getParam<BoundaryName, BoundaryName>(
"primary",
"secondary")),
284 _generate_mortar_mesh(getParam<bool>(
"generate_mortar_mesh")),
285 _mortar_dynamics(getParam<bool>(
"mortar_dynamics"))
288 if (
getParam<std::vector<BoundaryName>>(
"automatic_pairing_boundaries").size() > 1)
290 getParam<std::vector<BoundaryName>>(
"automatic_pairing_boundaries");
294 "For automatic selection of contact pairs (for particular geometries) in contact " 295 "action, 'automatic_pairing_distance' needs to be provided.");
299 "For automatic selection of contact pairs (for particular geometries) in contact " 300 "action, 'automatic_pairing_method' needs to be provided.");
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.");
309 "'primary' and 'secondary' surfaces or a list of boundaries for automatic pair " 310 "generation need to be provided.");
315 paramError(
"formulation",
"When using mortar, a vector of contact pairs cannot be used");
319 "The 'tangential_penalty' formulation can only be used with the 'coulomb' model");
330 paramError(
"model",
"The 'mortar_penalty' formulation does not support glued contact");
332 if (getParam<bool>(
"mortar_dynamics"))
334 "The 'mortar_penalty' formulation does not support implicit dynamic simulations");
336 if (getParam<bool>(
"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.");
347 paramError(
"model",
"The 'mortar' formulation does not support glued contact (yet)");
355 if (!getParam<bool>(
"mortar_dynamics"))
358 paramError(
"newmark_beta",
"newmark_beta can only be used with the mortar_dynamics option");
362 "newmark_gamma can only be used with the mortar_dynamics option");
369 "correct_edge_dropping",
370 "The 'correct_edge_dropping' option can only be used with the 'mortar' formulation " 375 "The 'use_dual' option can only be used with the 'mortar' formulation");
378 "The 'c_normal' option can only be used with the 'mortar' formulation");
381 "The 'c_tangential' option can only be used with the 'mortar' formulation");
384 "The 'mortar_dynamics' constraint option can only be used with the 'mortar' " 385 "formulation and in dynamic simulations using Newmark-beta");
392 "The 'secondary_gap_offset' option can only be used with the " 393 "'MechanicalContactConstraint'");
396 "The 'mapped_primary_gap_offset' option can only be used with the " 397 "'MechanicalContactConstraint'");
399 else if (getParam<bool>(
"ping_pong_protection"))
401 "The 'ping_pong_protection' option can only be used with the 'ranfs' formulation");
413 "Number of contact pairs in the contact action is zero. Please revise your input file.");
416 std::vector<std::pair<BoundaryName, BoundaryName>> lean_boundary_pairs;
421 auto it = std::find_if(lean_boundary_pairs.begin(),
422 lean_boundary_pairs.end(),
423 [&, primary_copy = primary, secondary_copy = secondary](
424 const std::pair<BoundaryName, BoundaryName> & lean_pair)
426 const bool match_one = lean_pair.second == secondary_copy &&
427 lean_pair.first == primary_copy;
428 const bool match_two = lean_pair.second == primary_copy &&
429 lean_pair.first == secondary_copy;
430 const bool exist = match_one || match_two;
434 if (it == lean_boundary_pairs.end())
435 lean_boundary_pairs.emplace_back(primary, secondary);
441 " has been removed from the contact interaction list due to " 442 "duplicates in the input file.");
454 if (!
_problem->getDisplacedProblem())
456 "Contact requires updated coordinates. Use the 'displacements = ...' parameter in the " 463 if (!
_problem->isSNESMFReuseBaseSetbyUser())
464 _problem->setSNESMFReuseBase(
false,
false);
475 if (!
_problem->getDisplacedProblem())
476 mooseError(
"Contact requires updated coordinates. Use the 'displacements = ...' line in the " 484 {
"secondary_gap_offset",
"mapped_primary_gap_offset",
"order"});
486 std::vector<VariableName> displacements =
487 getParam<std::vector<VariableName>>(
"displacements");
488 const auto order =
_problem->systemBaseNonlinear(0)
490 .variable_type(displacements[0])
493 params.
set<
MooseEnum>(
"order") = Utility::enum_to_string<Order>(OrderWrapper{order});
495 params.
set<std::vector<BoundaryName>>(
"boundary") = {contact_pair.second};
496 params.
set<BoundaryName>(
"paired_boundary") = contact_pair.first;
497 params.
set<AuxVariableName>(
"variable") =
"penetration";
499 params.
set<std::vector<VariableName>>(
"secondary_gap_offset") = {
500 getParam<VariableName>(
"secondary_gap_offset")};
502 params.
set<std::vector<VariableName>>(
"mapped_primary_gap_offset") = {
503 getParam<VariableName>(
"mapped_primary_gap_offset")};
504 params.
set<
bool>(
"use_displaced_mesh") =
true;
513 const unsigned int ndisp = getParam<std::vector<VariableName>>(
"displacements").size();
527 const std::string tangential_lagrange_multiplier_name = action_name +
"_tangential_lm";
528 const std::string tangential_lagrange_multiplier_3d_name =
529 action_name +
"_tangential_3d_lm";
531 params.
set<std::vector<VariableName>>(
"tangent_one") = {
532 tangential_lagrange_multiplier_name};
533 params.
set<std::vector<VariableName>>(
"tangent_two") = {
534 tangential_lagrange_multiplier_3d_name};
536 std::vector<std::string> disp_components({
"x",
"y",
"z"});
537 unsigned component_index = 0;
540 for (
const auto & disp_component : disp_components)
542 params.
set<AuxVariableName>(
"variable") =
_name +
"_tangent_" + disp_component;
543 params.
set<
unsigned int>(
"component") = component_index;
545 std::string
name =
_name +
"_mortar_frictional_pressure_" + disp_component +
"_" +
548 _problem->addAuxKernel(
"MortarFrictionalPressureVectorAux",
name, params);
557 std::vector<VariableName> displacements = getParam<std::vector<VariableName>>(
"displacements");
558 const auto order =
_problem->systemBaseNonlinear(0)
560 .variable_type(displacements[0])
565 var_params.
set<
MooseEnum>(
"order") = Utility::enum_to_string<Order>(OrderWrapper{order});
566 var_params.set<
MooseEnum>(
"family") =
"LAGRANGE";
568 _problem->addAuxVariable(
"MooseVariable",
"penetration", var_params);
573 var_params.
set<
MooseEnum>(
"order") = Utility::enum_to_string<Order>(OrderWrapper{order});
574 var_params.set<
MooseEnum>(
"family") =
"LAGRANGE";
576 _problem->addAuxVariable(
"MooseVariable",
"contact_pressure", var_params);
581 var_params.
set<
MooseEnum>(
"order") = Utility::enum_to_string<Order>(OrderWrapper{order});
582 var_params.set<
MooseEnum>(
"family") =
"LAGRANGE";
584 _problem->addAuxVariable(
"MooseVariable",
"nodal_area", var_params);
587 const unsigned int ndisp = getParam<std::vector<VariableName>>(
"displacements").size();
593 std::vector<std::string> disp_components({
"x",
"y",
"z"});
595 for (
const auto & disp_component : disp_components)
598 var_params.
set<
MooseEnum>(
"order") = Utility::enum_to_string<Order>(OrderWrapper{order});
599 var_params.set<
MooseEnum>(
"family") =
"LAGRANGE";
602 "MooseVariable",
_name +
"_tangent_" + disp_component, var_params);
616 std::vector<BoundaryName> secondary_boundary_vector;
617 for (
const auto *
const action : actions)
618 for (
const auto j :
index_range(action->_boundary_pairs))
619 secondary_boundary_vector.push_back(action->_boundary_pairs[
j].second);
621 var_params.set<std::vector<BoundaryName>>(
"boundary") = secondary_boundary_vector;
622 var_params.set<std::vector<VariableName>>(
"variable") = {
"nodal_area"};
624 mooseAssert(
_problem,
"Problem pointer is NULL");
626 var_params.set<
bool>(
"use_displaced_mesh") =
true;
628 _problem->addUserObject(
"NodalArea",
646 std::vector<BoundaryName> boundary_vector;
647 std::vector<BoundaryName> pair_boundary_vector;
649 for (
const auto *
const action : actions)
650 for (
const auto j :
index_range(action->_boundary_pairs))
652 boundary_vector.push_back(action->_boundary_pairs[
j].second);
653 pair_boundary_vector.push_back(action->_boundary_pairs[
j].first);
659 std::vector<VariableName> displacements = getParam<std::vector<VariableName>>(
"displacements");
660 const auto order =
_problem->systemBaseNonlinear(0)
662 .variable_type(displacements[0])
665 params.
set<
MooseEnum>(
"order") = Utility::enum_to_string<Order>(OrderWrapper{order});
666 params.
set<std::vector<BoundaryName>>(
"boundary") = boundary_vector;
667 params.
set<std::vector<BoundaryName>>(
"paired_boundary") = pair_boundary_vector;
668 params.
set<AuxVariableName>(
"variable") =
"contact_pressure";
670 params.
set<std::vector<VariableName>>(
"nodal_area") = {
"nodal_area"};
671 params.
set<
bool>(
"use_displaced_mesh") =
true;
673 std::string
name =
_name +
"_contact_pressure";
676 _problem->addAuxKernel(
"ContactPressureAux",
name, params);
687 params.set<
bool>(
"use_displaced_mesh") =
true;
689 const std::string primary_subdomain_name = action_name +
"_primary_subdomain";
690 const std::string secondary_subdomain_name = action_name +
"_secondary_subdomain";
691 params.set<BoundaryName>(
"primary_boundary") =
_boundary_pairs[0].first;
692 params.set<BoundaryName>(
"secondary_boundary") =
_boundary_pairs[0].second;
693 params.set<SubdomainName>(
"primary_subdomain") = primary_subdomain_name;
694 params.set<SubdomainName>(
"secondary_subdomain") = secondary_subdomain_name;
695 params.set<
bool>(
"use_petrov_galerkin") = getParam<bool>(
"use_petrov_galerkin");
705 std::vector<VariableName> displacements = getParam<std::vector<VariableName>>(
"displacements");
706 const unsigned int ndisp = displacements.size();
709 const std::string primary_subdomain_name = action_name +
"_primary_subdomain";
710 const std::string secondary_subdomain_name = action_name +
"_secondary_subdomain";
711 const std::string normal_lagrange_multiplier_name = action_name +
"_normal_lm";
712 const std::string tangential_lagrange_multiplier_name = action_name +
"_tangential_lm";
713 const std::string tangential_lagrange_multiplier_3d_name = action_name +
"_tangential_3d_lm";
714 const std::string auxiliary_lagrange_multiplier_name = action_name +
"_aux_lm";
724 const MeshGeneratorName primary_name = primary_subdomain_name +
"_generator";
725 const MeshGeneratorName secondary_name = secondary_subdomain_name +
"_generator";
730 primary_params.
set<SubdomainName>(
"new_block_name") = primary_subdomain_name;
731 secondary_params.set<SubdomainName>(
"new_block_name") = secondary_subdomain_name;
733 primary_params.set<std::vector<BoundaryName>>(
"sidesets") = {
_boundary_pairs[0].first};
734 secondary_params.set<std::vector<BoundaryName>>(
"sidesets") = {
_boundary_pairs[0].second};
742 const auto addLagrangeMultiplier =
743 [
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)
753 if (!add_aux_lm || penalty_traction)
756 mooseAssert(
_problem->systemBaseNonlinear(0).hasVariable(displacements[0]),
757 "Displacement variable is missing");
758 const auto primal_type =
759 _problem->systemBaseNonlinear(0).system().variable_type(displacements[0]);
761 const int lm_order = primal_type.order.get_order();
765 params.
set<
MooseEnum>(
"family") = Utility::enum_to_string<FEFamily>(primal_type.family);
766 params.
set<
MooseEnum>(
"order") = Utility::enum_to_string<Order>(OrderWrapper{lm_order});
769 mooseError(
"Invalid bases for mortar contact.");
771 params.
set<std::vector<SubdomainName>>(
"block") = {secondary_subdomain_name};
772 if (!(add_aux_lm || penalty_traction))
773 params.
set<std::vector<Real>>(
"scaling") = {scaling_factor};
777 if (add_aux_lm || penalty_traction)
778 _problem->addAuxVariable(var_type, variable_name, params);
780 _problem->addVariable(var_type, variable_name, params);
785 addLagrangeMultiplier(
786 normal_lagrange_multiplier_name, getParam<Real>(
"normal_lm_scaling"),
false,
false);
790 addLagrangeMultiplier(tangential_lagrange_multiplier_name,
791 getParam<Real>(
"tangential_lm_scaling"),
795 addLagrangeMultiplier(tangential_lagrange_multiplier_3d_name,
796 getParam<Real>(
"tangential_lm_scaling"),
801 if (getParam<bool>(
"use_petrov_galerkin"))
802 addLagrangeMultiplier(auxiliary_lagrange_multiplier_name, 1.0,
true,
false);
808 addLagrangeMultiplier(auxiliary_lagrange_multiplier_name, 1.0,
false,
true);
815 !dynamic_cast<AugmentedLagrangianContactProblemInterface *>(
_problem.get()))
817 const std::vector<std::string> params = {
"penalty_multiplier",
818 "penalty_multiplier_friction",
819 "al_penetration_tolerance",
820 "al_incremental_slip_tolerance",
821 "al_frictional_force_tolerance"};
822 for (
const auto & param : params)
825 "Augmented Lagrange parameter was specified, but the selected problem type " 826 "does not support Augmented Lagrange iterations.");
834 var_params.set<BoundaryName>(
"secondary_boundary") =
_boundary_pairs[0].second;
835 var_params.set<SubdomainName>(
"primary_subdomain") = primary_subdomain_name;
836 var_params.set<SubdomainName>(
"secondary_subdomain") = secondary_subdomain_name;
837 var_params.set<std::vector<VariableName>>(
"disp_x") = {displacements[0]};
838 var_params.set<
bool>(
"correct_edge_dropping") = getParam<bool>(
"correct_edge_dropping");
839 var_params.set<std::vector<VariableName>>(
"disp_y") = {displacements[1]};
841 var_params.set<std::vector<VariableName>>(
"disp_z") = {displacements[2]};
842 var_params.set<
bool>(
"use_displaced_mesh") =
true;
843 var_params.set<std::vector<VariableName>>(
"lm_variable") = {normal_lagrange_multiplier_name};
844 var_params.set<
bool>(
"use_petrov_galerkin") = getParam<bool>(
"use_petrov_galerkin");
845 if (getParam<bool>(
"use_petrov_galerkin"))
846 var_params.set<std::vector<VariableName>>(
"aux_lm") = {auxiliary_lagrange_multiplier_name};
849 "LMWeightedGapUserObject",
"lm_weightedgap_object_" +
name(), var_params);
855 var_params.set<BoundaryName>(
"secondary_boundary") =
_boundary_pairs[0].second;
856 var_params.set<SubdomainName>(
"primary_subdomain") = primary_subdomain_name;
857 var_params.set<SubdomainName>(
"secondary_subdomain") = secondary_subdomain_name;
858 var_params.set<std::vector<VariableName>>(
"disp_x") = {displacements[0]};
859 var_params.set<
bool>(
"correct_edge_dropping") = getParam<bool>(
"correct_edge_dropping");
860 var_params.set<std::vector<VariableName>>(
"disp_y") = {displacements[1]};
862 var_params.set<std::vector<VariableName>>(
"disp_z") = {displacements[2]};
864 var_params.set<VariableName>(
"secondary_variable") = displacements[0];
865 var_params.set<
bool>(
"use_displaced_mesh") =
true;
866 var_params.set<std::vector<VariableName>>(
"lm_variable_normal") = {
867 normal_lagrange_multiplier_name};
868 var_params.set<std::vector<VariableName>>(
"lm_variable_tangential_one") = {
869 tangential_lagrange_multiplier_name};
871 var_params.set<std::vector<VariableName>>(
"lm_variable_tangential_two") = {
872 tangential_lagrange_multiplier_3d_name};
873 var_params.set<
bool>(
"use_petrov_galerkin") = getParam<bool>(
"use_petrov_galerkin");
874 if (getParam<bool>(
"use_petrov_galerkin"))
875 var_params.set<std::vector<VariableName>>(
"aux_lm") = {auxiliary_lagrange_multiplier_name};
878 "LMWeightedVelocitiesUserObject",
"lm_weightedvelocities_object_" +
name(), var_params);
886 var_params.set<BoundaryName>(
"secondary_boundary") =
_boundary_pairs[0].second;
887 var_params.set<SubdomainName>(
"primary_subdomain") = primary_subdomain_name;
888 var_params.set<SubdomainName>(
"secondary_subdomain") = secondary_subdomain_name;
889 var_params.set<std::vector<VariableName>>(
"disp_x") = {displacements[0]};
890 var_params.set<
bool>(
"correct_edge_dropping") = getParam<bool>(
"correct_edge_dropping");
891 var_params.set<std::vector<VariableName>>(
"disp_y") = {displacements[1]};
892 var_params.set<
Real>(
"penalty") = getParam<Real>(
"penalty");
895 var_params.set<
Real>(
"max_penalty_multiplier") = getParam<Real>(
"max_penalty_multiplier");
896 var_params.set<
MooseEnum>(
"adaptivity_penalty_normal") =
897 getParam<MooseEnum>(
"adaptivity_penalty_normal");
899 var_params.set<
Real>(
"penetration_tolerance") = getParam<Real>(
"al_penetration_tolerance");
901 var_params.set<
Real>(
"penalty_multiplier") = getParam<Real>(
"penalty_multiplier");
904 var_params.set<
bool>(
"use_physical_gap") =
true;
907 var_params.set<std::vector<VariableName>>(
"aux_lm") = {auxiliary_lagrange_multiplier_name};
910 var_params.set<std::vector<VariableName>>(
"disp_z") = {displacements[2]};
911 var_params.set<
bool>(
"use_displaced_mesh") =
true;
914 "PenaltyWeightedGapUserObject",
"penalty_weightedgap_object_" +
name(), var_params);
921 var_params.set<BoundaryName>(
"secondary_boundary") =
_boundary_pairs[0].second;
922 var_params.set<SubdomainName>(
"primary_subdomain") = primary_subdomain_name;
923 var_params.set<SubdomainName>(
"secondary_subdomain") = secondary_subdomain_name;
924 var_params.set<std::vector<VariableName>>(
"disp_x") = {displacements[0]};
925 var_params.set<
bool>(
"correct_edge_dropping") = getParam<bool>(
"correct_edge_dropping");
926 var_params.set<std::vector<VariableName>>(
"disp_y") = {displacements[1]};
928 var_params.set<std::vector<VariableName>>(
"disp_z") = {displacements[2]};
930 var_params.set<VariableName>(
"secondary_variable") = displacements[0];
931 var_params.set<
bool>(
"use_displaced_mesh") =
true;
932 var_params.set<
Real>(
"friction_coefficient") = getParam<Real>(
"friction_coefficient");
933 var_params.set<
Real>(
"penalty") = getParam<Real>(
"penalty");
934 var_params.set<
Real>(
"penalty_friction") = getParam<Real>(
"penalty_friction");
937 var_params.set<
Real>(
"max_penalty_multiplier") = getParam<Real>(
"max_penalty_multiplier");
938 var_params.set<
MooseEnum>(
"adaptivity_penalty_normal") =
939 getParam<MooseEnum>(
"adaptivity_penalty_normal");
940 var_params.set<
MooseEnum>(
"adaptivity_penalty_friction") =
941 getParam<MooseEnum>(
"adaptivity_penalty_friction");
943 var_params.set<
Real>(
"penetration_tolerance") = getParam<Real>(
"al_penetration_tolerance");
945 var_params.set<
Real>(
"penalty_multiplier") = getParam<Real>(
"penalty_multiplier");
947 var_params.set<
Real>(
"penalty_multiplier_friction") =
948 getParam<Real>(
"penalty_multiplier_friction");
951 var_params.set<
Real>(
"slip_tolerance") = getParam<Real>(
"al_incremental_slip_tolerance");
954 var_params.set<
bool>(
"use_physical_gap") =
true;
957 var_params.set<std::vector<VariableName>>(
"aux_lm") = {auxiliary_lagrange_multiplier_name};
959 var_params.applySpecificParameters(
parameters(),
960 {
"friction_coefficient",
"penalty",
"penalty_friction"});
963 "PenaltyFrictionUserObject",
"penalty_friction_object_" +
name(), var_params);
973 std::string mortar_constraint_name;
976 mortar_constraint_name =
"ComputeWeightedGapLMMechanicalContact";
978 mortar_constraint_name =
"ComputeDynamicWeightedGapLMMechanicalContact";
983 parameters(), {
"newmark_beta",
"newmark_gamma",
"capture_tolerance",
"wear_depth"});
986 params.
set<UserObjectName>(
"weighted_gap_uo") =
"lm_weightedgap_object_" +
name();
990 params.
set<SubdomainName>(
"primary_subdomain") = primary_subdomain_name;
991 params.
set<SubdomainName>(
"secondary_subdomain") = secondary_subdomain_name;
992 params.
set<NonlinearVariableName>(
"variable") = normal_lagrange_multiplier_name;
993 params.
set<std::vector<VariableName>>(
"disp_x") = {displacements[0]};
994 params.
set<
Real>(
"c") = getParam<Real>(
"c_normal");
997 params.
set<std::vector<VariableName>>(
"disp_y") = {displacements[1]};
999 params.
set<std::vector<VariableName>>(
"disp_z") = {displacements[2]};
1001 params.
set<
bool>(
"use_displaced_mesh") =
true;
1004 {
"correct_edge_dropping",
1006 "extra_vector_tags",
1007 "absolute_value_vector_tags"});
1010 mortar_constraint_name, action_name +
"_normal_lm_weighted_gap", params);
1016 std::string mortar_constraint_name;
1019 mortar_constraint_name =
"ComputeFrictionalForceLMMechanicalContact";
1021 mortar_constraint_name =
"ComputeDynamicFrictionalForceLMMechanicalContact";
1026 parameters(), {
"newmark_beta",
"newmark_gamma",
"capture_tolerance",
"wear_depth"});
1029 params.
set<UserObjectName>(
"weighted_gap_uo") =
"lm_weightedvelocities_object_" +
name();
1030 params.
set<UserObjectName>(
"weighted_velocities_uo") =
1031 "lm_weightedvelocities_object_" +
name();
1034 params.
set<
bool>(
"correct_edge_dropping") = getParam<bool>(
"correct_edge_dropping");
1037 params.
set<SubdomainName>(
"primary_subdomain") = primary_subdomain_name;
1038 params.
set<SubdomainName>(
"secondary_subdomain") = secondary_subdomain_name;
1039 params.
set<
bool>(
"use_displaced_mesh") =
true;
1040 params.
set<
Real>(
"c_t") = getParam<Real>(
"c_tangential");
1041 params.
set<
Real>(
"c") = getParam<Real>(
"c_normal");
1042 params.
set<
bool>(
"normalize_c") = getParam<bool>(
"normalize_c");
1043 params.
set<
bool>(
"compute_primal_residuals") =
false;
1045 params.
set<std::vector<VariableName>>(
"disp_x") = {displacements[0]};
1048 params.
set<std::vector<VariableName>>(
"disp_y") = {displacements[1]};
1050 params.
set<std::vector<VariableName>>(
"disp_z") = {displacements[2]};
1052 params.
set<NonlinearVariableName>(
"variable") = normal_lagrange_multiplier_name;
1053 params.
set<std::vector<VariableName>>(
"friction_lm") = {tangential_lagrange_multiplier_name};
1056 params.
set<std::vector<VariableName>>(
"friction_lm_dir") = {
1057 tangential_lagrange_multiplier_3d_name};
1059 params.
set<
Real>(
"mu") = getParam<Real>(
"friction_coefficient");
1061 {
"extra_vector_tags",
"absolute_value_vector_tags"});
1063 _problem->addConstraint(mortar_constraint_name, action_name +
"_tangential_lm", params);
1067 const auto addMechanicalContactConstraints =
1068 [
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)
1077 params.
set<
bool>(
"correct_edge_dropping") = getParam<bool>(
"correct_edge_dropping");
1080 params.
set<SubdomainName>(
"primary_subdomain") = primary_subdomain_name;
1081 params.
set<SubdomainName>(
"secondary_subdomain") = secondary_subdomain_name;
1084 params.
set<NonlinearVariableName>(
"variable") = variable_name;
1086 params.
set<
bool>(
"use_displaced_mesh") =
true;
1087 params.
set<
bool>(
"compute_lm_residuals") =
false;
1091 if (is_additional_frictional_constraint)
1094 {
"extra_vector_tags",
"absolute_value_vector_tags"});
1096 for (
unsigned int i = 0; i < displacements.size(); ++i)
1100 params.
set<VariableName>(
"secondary_variable") = displacements[i];
1105 params.
set<UserObjectName>(
"weighted_gap_uo") =
"lm_weightedgap_object_" +
name();
1108 params.
set<UserObjectName>(
"weighted_gap_uo") =
"lm_weightedvelocities_object_" +
name();
1110 params.
set<UserObjectName>(
"weighted_velocities_uo") =
1111 "lm_weightedvelocities_object_" +
name();
1114 params.
set<UserObjectName>(
"weighted_gap_uo") =
"penalty_weightedgap_object_" +
name();
1117 params.
set<UserObjectName>(
"weighted_gap_uo") =
"penalty_friction_object_" +
name();
1119 params.
set<UserObjectName>(
"weighted_velocities_uo") =
1120 "penalty_friction_object_" +
name();
1122 _problem->addConstraint(constraint_type, constraint_name, params);
1128 addMechanicalContactConstraints(normal_lagrange_multiplier_name,
1129 action_name +
"_normal_constraint_",
1130 "NormalMortarMechanicalContact",
1136 addMechanicalContactConstraints(tangential_lagrange_multiplier_name,
1137 action_name +
"_tangential_constraint_",
1138 "TangentialMortarMechanicalContact",
1142 addMechanicalContactConstraints(tangential_lagrange_multiplier_3d_name,
1143 action_name +
"_tangential_constraint_3d_",
1144 "TangentialMortarMechanicalContact",
1156 if (getParam<MooseEnum>(
"automatic_pairing_method").getEnum<
ProximityMethod>() ==
1159 else if (getParam<MooseEnum>(
"automatic_pairing_method").getEnum<ProximityMethod>() ==
1168 std::vector<VariableName> displacements = getParam<std::vector<VariableName>>(
"displacements");
1169 const unsigned int ndisp = displacements.size();
1171 std::string constraint_type;
1174 constraint_type =
"RANFSNormalMechanicalContact";
1176 constraint_type =
"MechanicalContactConstraint";
1182 "secondary_gap_offset",
1183 "mapped_primary_gap_offset",
1187 const auto order =
_problem->systemBaseNonlinear(0)
1189 .variable_type(displacements[0])
1192 params.
set<std::vector<VariableName>>(
"displacements") = displacements;
1193 params.
set<
bool>(
"use_displaced_mesh") =
true;
1194 params.
set<
MooseEnum>(
"order") = Utility::enum_to_string<Order>(OrderWrapper{order});
1200 params.
set<std::vector<VariableName>>(
"nodal_area") = {
"nodal_area"};
1201 params.
set<BoundaryName>(
"boundary") = contact_pair.first;
1203 params.
set<std::vector<VariableName>>(
"secondary_gap_offset") = {
1204 getParam<VariableName>(
"secondary_gap_offset")};
1206 params.
set<std::vector<VariableName>>(
"mapped_primary_gap_offset") = {
1207 getParam<VariableName>(
"mapped_primary_gap_offset")};
1210 for (
unsigned int i = 0; i < ndisp; ++i)
1218 params.
set<
unsigned int>(
"component") = i;
1220 params.
set<BoundaryName>(
"primary") = contact_pair.first;
1221 params.
set<BoundaryName>(
"secondary") = contact_pair.second;
1222 params.
set<NonlinearVariableName>(
"variable") = displacements[i];
1223 params.
set<std::vector<VariableName>>(
"primary_variable") = {displacements[i]};
1225 {
"extra_vector_tags",
"absolute_value_vector_tags"});
1226 _problem->addConstraint(constraint_type,
name, params);
1234 inline const Point &
1237 return *(item.first);
1243 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.");
1247 mooseError(
"Failed to obtain mesh for automatically generating contact pairs.");
1249 if (!
_mesh->getMesh().is_serial())
1251 "automatic_pairing_boundaries",
1252 "The generation of automatic contact pairs in the contact action requires a serial mesh.");
1255 std::vector<BoundaryID> _automatic_pairing_boundaries_id;
1257 _automatic_pairing_boundaries_id.emplace_back(
_mesh->getBoundaryID(sideset_name));
1260 std::vector<NodeBoundaryIDInfo> node_boundary_id_vector;
1265 for (
const auto & bnode : bnd_nodes)
1267 const BoundaryID boundary_id = bnode->_bnd_id;
1268 const Node * node_ptr = bnode->_node;
1271 auto it = std::find(_automatic_pairing_boundaries_id.begin(),
1272 _automatic_pairing_boundaries_id.end(),
1275 if (it != _automatic_pairing_boundaries_id.end())
1276 node_boundary_id_vector.emplace_back(node_ptr, boundary_id);
1280 std::sort(node_boundary_id_vector.begin(),
1281 node_boundary_id_vector.end(),
1283 {
return first_pair.second < second_pair.second; });
1286 using KDTreeType = nanoflann::KDTreeSingleIndexAdaptor<
1287 nanoflann::L2_Simple_Adaptor<Real, PointListAdaptor<NodeBoundaryIDInfo>,
Real, std::size_t>,
1293 const unsigned int max_leaf_size = 20;
1297 node_boundary_id_vector.end());
1298 auto kd_tree = std::make_unique<KDTreeType>(
1299 LIBMESH_DIM, point_list, nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size));
1302 mooseError(
"Internal error. KDTree was not properly initialized in the contact action.");
1304 kd_tree->buildIndex();
1308 std::vector<nanoflann::ResultItem<std::size_t, Real>> ret_matches;
1310 const auto radius_for_search = getParam<Real>(
"automatic_pairing_distance");
1313 for (
const auto & pair : node_boundary_id_vector)
1316 ret_matches.clear();
1319 const Point search_point = *pair.first;
1322 kd_tree->radiusSearch(
1323 &(search_point)(0), radius_for_search * radius_for_search, ret_matches, search_params);
1325 for (
auto & match_pair : ret_matches)
1327 const auto & match = node_boundary_id_vector[match_pair.first];
1334 auto it = std::find(_automatic_pairing_boundaries_id.begin(),
1335 _automatic_pairing_boundaries_id.end(),
1339 if (match.second == pair.second)
1345 if (it != _automatic_pairing_boundaries_id.end())
1347 const auto index_one = cast_int<int>(it - _automatic_pairing_boundaries_id.begin());
1348 auto it_other = std::find(_automatic_pairing_boundaries_id.begin(),
1349 _automatic_pairing_boundaries_id.end(),
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.");
1356 const auto index_two = cast_int<int>(it_other - _automatic_pairing_boundaries_id.begin());
1358 if (pair.second > match.second)
1372 "The following boundary pairs were created by the contact action using nodal proximity: ");
1375 "Primary boundary ID: ", primary,
" and secondary boundary ID: ", secondary,
".");
1381 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.");
1385 mooseError(
"Failed to obtain mesh for automatically generating contact pairs.");
1387 if (!
_mesh->getMesh().is_serial())
1389 "automatic_pairing_boundaries",
1390 "The generation of automatic contact pairs in the contact action requires a serial mesh.");
1393 std::vector<std::pair<BoundaryName, Point>> automatic_pairing_boundaries_cog;
1394 const auto & sideset_ids =
_mesh->meshSidesetIds();
1396 const auto & bnd_to_elem_map =
_mesh->getBoundariesToActiveSemiLocalElemIds();
1401 const auto find_set = sideset_ids.find(
_mesh->getBoundaryID(sideset_name));
1402 if (find_set == sideset_ids.end())
1405 " is not defined as a sideset in the mesh.");
1407 auto dofs_set = bnd_to_elem_map.find(
_mesh->getBoundaryID(sideset_name));
1410 Point center_of_gravity(0, 0, 0);
1411 Real accumulated_sideset_area(0);
1414 std::unique_ptr<const Elem> side_ptr;
1415 const std::unordered_set<dof_id_type> & bnd_elems = dofs_set->second;
1417 for (
auto elem_id : bnd_elems)
1419 const Elem * elem =
_mesh->elemPtr(elem_id);
1420 unsigned int side =
_mesh->sideWithBoundaryID(elem,
_mesh->getBoundaryID(sideset_name));
1423 elem->side_ptr(side_ptr, side);
1426 const auto side_area = side_ptr->volume();
1429 const auto side_position = side_ptr->true_centroid();
1431 center_of_gravity += side_position * side_area;
1432 accumulated_sideset_area += side_area;
1436 center_of_gravity /= accumulated_sideset_area;
1439 automatic_pairing_boundaries_cog.emplace_back(sideset_name, center_of_gravity);
1443 std::vector<std::pair<std::pair<BoundaryName, BoundaryName>,
Real>> pairs_distances;
1446 for (std::size_t i = 0; i < automatic_pairing_boundaries_cog.size() - 1; i++)
1447 for (std::size_t
j = i + 1;
j < automatic_pairing_boundaries_cog.size();
j++)
1449 const Point & distance_vector =
1450 automatic_pairing_boundaries_cog[i].second - automatic_pairing_boundaries_cog[
j].second;
1452 if (automatic_pairing_boundaries_cog[i].first != automatic_pairing_boundaries_cog[
j].first)
1455 const std::pair pair = std::make_pair(automatic_pairing_boundaries_cog[i].first,
1456 automatic_pairing_boundaries_cog[
j].first);
1457 pairs_distances.emplace_back(std::make_pair(pair,
distance));
1461 const auto automatic_pairing_distance = getParam<Real>(
"automatic_pairing_distance");
1464 std::vector<std::pair<std::pair<BoundaryName, BoundaryName>,
Real>> lean_pairs_distances;
1465 for (
const auto & pair_distance : pairs_distances)
1466 if (pair_distance.second <= automatic_pairing_distance)
1468 lean_pairs_distances.emplace_back(pair_distance);
1470 pair_distance.first.first,
1472 pair_distance.first.second,
1473 ", with a relative distance of ",
1474 pair_distance.second);
1478 for (
const auto & lean_pairs_distance : lean_pairs_distances)
1483 if (
_mesh->getBoundaryID(lean_pairs_distance.first.first) >
1484 _mesh->getBoundaryID(lean_pairs_distance.first.second))
1486 {lean_pairs_distance.first.first, lean_pairs_distance.first.second});
1489 {lean_pairs_distance.first.second, lean_pairs_distance.first.first});
1499 return MooseEnum(
"frictionless glued coulomb",
"frictionless");
1512 "ranfs kinematic penalty augmented_lagrange tangential_penalty mortar mortar_penalty",
1515 formulations.addDocumentation(
1517 "Reduced Active Nonlinear Function Set scheme for node-on-face contact. Provides exact " 1518 "enforcement without Lagrange multipliers or penalty terms.");
1519 formulations.addDocumentation(
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 formulations.addDocumentation(
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 " 1529 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 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 formulations.addDocumentation(
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 formulations.addDocumentation(
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.");
1552 return formulations;
1558 return MooseEnum(
"Constraint",
"Constraint");
1564 return MooseEnum(
"edge_based nodal_normal_based",
"");
1574 "Method to use to smooth normals");
1576 "normal_smoothing_distance",
1577 "Distance from edge in parametric coordinates over which to smooth contact normal");
bool isUltimateMaster() const
static std::string variableType(const FEType &fe_type, const bool is_fv=false, const bool is_array=false)
void mooseInfo(Args &&... args) const
const Point & getPoint(const PointObject &item) const
InputParameters getValidParams(const std::string &name) const
void mooseInfoRepeated(Args &&... args)
const ExecFlagType EXEC_TIMESTEP_END
Real distance(const Point &p)
virtual const std::string & name() const
std::string shortName(const std::string &name)
static InputParameters validParams()
bool isParamValid(const std::string &name) const
static InputParameters validParams()
const ExecFlagType EXEC_TIMESTEP_BEGIN
boundary_id_type BoundaryID
const T & getParam(const std::string &name) const
const std::string & _current_task
const MooseMesh * masterMesh() const
void paramError(const std::string ¶m, Args... args) const
const ExecFlagType EXEC_LINEAR
std::string stringify(const T &t)
const MeshGenerator & appendMeshGenerator(const std::string &type, const std::string &name, InputParameters params)
std::pair< T, U > ResultItem
const ExecFlagType EXEC_NONLINEAR
std::shared_ptr< MooseMesh > & _mesh
static FEType feType(const InputParameters ¶ms)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void mooseError(Args &&... args) const
std::shared_ptr< FEProblemBase > & _problem
const InputParameters & parameters() const
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
StoredRange< MooseMesh::const_bnd_node_iterator, const BndNode *> ConstBndNodeRange
std::vector< const T *> getActions()
bool isRecovering() const
SearchParams SearchParameters
auto index_range(const T &sizable)
const ExecFlagType EXEC_INITIAL