https://mooseframework.inl.gov
ComputeWeightedGapLMMechanicalContact.C
Go to the documentation of this file.
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 
11 #include "DisplacedProblem.h"
12 #include "Assembly.h"
13 #include "MortarContactUtils.h"
14 #include "WeightedGapUserObject.h"
15 #include "metaphysicl/dualsemidynamicsparsenumberarray.h"
16 #include "metaphysicl/parallel_dualnumber.h"
17 #include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
18 #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
19 #include "timpi/parallel_sync.h"
20 
22 
23 namespace
24 {
25 const InputParameters &
26 assignVarsInParams(const InputParameters & params_in)
27 {
28  InputParameters & ret = const_cast<InputParameters &>(params_in);
29  const auto & disp_x_name = ret.get<std::vector<VariableName>>("disp_x");
30  if (disp_x_name.size() != 1)
31  mooseError("We require that the disp_x parameter have exactly one coupled name");
32 
33  // We do this so we don't get any variable errors during MortarConstraint(Base) construction
34  ret.set<VariableName>("secondary_variable") = disp_x_name[0];
35  ret.set<VariableName>("primary_variable") = disp_x_name[0];
36 
37  return ret;
38 }
39 }
40 
43 {
45  params.addClassDescription("Computes the weighted gap that will later be used to enforce the "
46  "zero-penetration mechanical contact conditions");
47  params.suppressParameter<VariableName>("secondary_variable");
48  params.suppressParameter<VariableName>("primary_variable");
49  params.addRequiredCoupledVar("disp_x", "The x displacement variable");
50  params.addRequiredCoupledVar("disp_y", "The y displacement variable");
51  params.addCoupledVar("disp_z", "The z displacement variable");
52  params.addParam<Real>(
53  "c", 1e6, "Parameter for balancing the size of the gap and contact pressure");
54  params.addParam<bool>(
55  "normalize_c",
56  false,
57  "Whether to normalize c by weighting function norm. When unnormalized "
58  "the value of c effectively depends on element size since in the constraint we compare nodal "
59  "Lagrange Multiplier values to integrated gap values (LM nodal value is independent of "
60  "element size, where integrated values are dependent on element size).");
61  params.set<bool>("use_displaced_mesh") = true;
62  params.set<bool>("interpolate_normals") = false;
63  params.addRequiredParam<UserObjectName>("weighted_gap_uo", "The weighted gap user object");
64  return params;
65 }
66 
68  const InputParameters & parameters)
69  : ADMortarConstraint(assignVarsInParams(parameters)),
70  _secondary_disp_x(adCoupledValue("disp_x")),
71  _primary_disp_x(adCoupledNeighborValue("disp_x")),
72  _secondary_disp_y(adCoupledValue("disp_y")),
73  _primary_disp_y(adCoupledNeighborValue("disp_y")),
74  _has_disp_z(isCoupled("disp_z")),
75  _secondary_disp_z(_has_disp_z ? &adCoupledValue("disp_z") : nullptr),
76  _primary_disp_z(_has_disp_z ? &adCoupledNeighborValue("disp_z") : nullptr),
77  _c(getParam<Real>("c")),
78  _normalize_c(getParam<bool>("normalize_c")),
79  _nodal(getVar("disp_x", 0)->feType().family == LAGRANGE),
80  _disp_x_var(getVar("disp_x", 0)),
81  _disp_y_var(getVar("disp_y", 0)),
82  _disp_z_var(_has_disp_z ? getVar("disp_z", 0) : nullptr),
83  _weighted_gap_uo(getUserObject<WeightedGapUserObject>("weighted_gap_uo"))
84 {
85  if (!getParam<bool>("use_displaced_mesh"))
86  paramError(
87  "use_displaced_mesh",
88  "'use_displaced_mesh' must be true for the ComputeWeightedGapLMMechanicalContact object");
89 
90  if (!_var->isNodal())
91  if (_var->feType().order != static_cast<Order>(0))
92  mooseError("Normal contact constraints only support elemental variables of CONSTANT order");
93 }
94 
96 {
97  mooseError("We should never call computeQpResidual for ComputeWeightedGapLMMechanicalContact");
98 }
99 
100 void
102 {
103 }
104 
105 void
107 {
108 }
109 
110 void
112 {
113 }
114 
115 void
117 {
118 }
119 
120 void
122 {
123 }
124 
125 void
127 {
128  // During "computeResidual" and "computeJacobian" we are actually just computing properties on the
129  // mortar segment element mesh. We are *not* actually assembling into the residual/Jacobian. For
130  // the zero-penetration constraint, the property of interest is the map from node to weighted gap.
131  // Computation of the properties proceeds identically for residual and Jacobian evaluation hence
132  // why we simply call computeResidual here. We will assemble into the residual/Jacobian later from
133  // the post() method
134  computeResidual(mortar_type);
135 }
136 
137 void
139 {
140  parallel_object_only();
141 
142  const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
143 
144  for (const auto & [dof_object, weighted_gap_pr] : dof_to_weighted_gap)
145  {
146  if (dof_object->processor_id() != this->processor_id())
147  continue;
148 
149  const auto & [weighted_gap, normalization] = weighted_gap_pr;
150  _weighted_gap_ptr = &weighted_gap;
151  _normalization_ptr = &normalization;
152 
153  enforceConstraintOnDof(dof_object);
154  }
155 }
156 
157 void
159  const std::unordered_set<const Node *> & inactive_lm_nodes)
160 {
161  const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
162 
163  for (const auto & [dof_object, weighted_gap_pr] : dof_to_weighted_gap)
164  {
165  if ((inactive_lm_nodes.find(static_cast<const Node *>(dof_object)) !=
166  inactive_lm_nodes.end()) ||
167  (dof_object->processor_id() != this->processor_id()))
168  continue;
169 
170  const auto & [weighted_gap, normalization] = weighted_gap_pr;
171  _weighted_gap_ptr = &weighted_gap;
172  _normalization_ptr = &normalization;
173 
174  enforceConstraintOnDof(dof_object);
175  }
176 }
177 
178 void
180 {
181  const auto & weighted_gap = *_weighted_gap_ptr;
182  const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
183 
184  const auto dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
185  ADReal lm_value = (*_sys.currentSolution())(dof_index);
186  Moose::derivInsert(lm_value.derivatives(), dof_index, 1.);
187 
188  const ADReal dof_residual = std::min(lm_value, weighted_gap * c);
189 
191  std::array<ADReal, 1>{{dof_residual}},
192  std::array<dof_id_type, 1>{{dof_index}},
193  _var->scalingFactor());
194 }
LAGRANGE
virtual const NumericVector< Number > *const & currentSolution() const=0
const libMesh::FEType & feType() const
Order
void addResidualsAndJacobian(Assembly &assembly, const Residuals &residuals, const Indices &dof_indices, Real scaling_factor)
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
void mooseError(Args &&... args)
unsigned int number() const
std::vector< std::pair< R1, R2 > > get(const std::string &param1, const std::string &param2) const
T & set(const std::string &name, bool quiet_mode=false)
OrderWrapper order
ComputeWeightedGapLMMechanicalContact(const InputParameters &parameters)
const Real _c
This factor multiplies the weighted gap.
const WeightedGapUserObject & _weighted_gap_uo
The weighted gap user object.
DualNumber< Real, DNDerivativeType, true > ADReal
virtual void computeResidual() override
void addRequiredParam(const std::string &name, const std::string &doc_string)
virtual void computeQpProperties()
Computes properties that are functions only of the current quadrature point (_qp), e.g.
void suppressParameter(const std::string &name)
SystemBase & _sys
virtual void computeQpIProperties()
Computes properties that are functions both of _qp and _i, for example the weighted gap...
registerMooseObject("ContactApp", ComputeWeightedGapLMMechanicalContact)
ADReal computeQpResidual(Moose::MortarType mortar_type) final
const ADReal * _weighted_gap_ptr
A pointer members that can be used to help avoid copying ADReals.
void paramError(const std::string &param, Args... args) const
unsigned int number() const
const std::unordered_map< const DofObject *, std::pair< ADReal, Real > > & dofToWeightedGap() const
Get the degree of freedom to weighted gap information.
Assembly & _assembly
void incorrectEdgeDroppingPost(const std::unordered_set< const Node *> &inactive_lm_nodes) override
Copy of the post routine but that skips assembling inactive nodes.
void addCoupledVar(const std::string &name, const std::string &doc_string)
bool isNodal() const override
void addRequiredCoupledVar(const std::string &name, const std::string &doc_string)
MooseVariable *const _var
Creates dof object to weighted gap map.
const bool _normalize_c
Whether to normalize weighted gap by weighting function norm.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void computeJacobian() override
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
static InputParameters validParams()
virtual void enforceConstraintOnDof(const DofObject *const dof)
Method called from post().
void derivInsert(SemiDynamicSparseNumberArray< Real, libMesh::dof_id_type, NWrapper< N >> &derivs, libMesh::dof_id_type index, Real value)
processor_id_type processor_id() const
Computes the weighted gap that will later be used to enforce the zero-penetration mechanical contact ...
void scalingFactor(const std::vector< Real > &factor)