https://mooseframework.inl.gov
ComputeDynamicWeightedGapLMMechanicalContact.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 "MortarContactUtils.h"
12 #include "DisplacedProblem.h"
13 #include "Assembly.h"
14 
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 assignVarsInParamsDyn(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 }
42 {
44  params.addClassDescription(
45  "Computes the normal contact mortar constraints for dynamic simulations");
46  params.addRangeCheckedParam<Real>("capture_tolerance",
47  1.0e-5,
48  "capture_tolerance>=0",
49  "Parameter describing a gap threshold for the application of "
50  "the persistency constraint in dynamic simulations.");
51  params.addCoupledVar("wear_depth",
52  "The name of the mortar auxiliary variable that is used to modify the "
53  "weighted gap definition");
55  "newmark_beta", "newmark_beta > 0", "Beta parameter for the Newmark time integrator");
57  "newmark_gamma", "newmark_gamma >= 0.0", "Gamma parameter for the Newmark time integrator");
58  params.suppressParameter<VariableName>("secondary_variable");
59  params.suppressParameter<VariableName>("primary_variable");
60  params.addRequiredCoupledVar("disp_x", "The x displacement variable");
61  params.addRequiredCoupledVar("disp_y", "The y displacement variable");
62  params.addCoupledVar("disp_z", "The z displacement variable");
63  params.addParam<Real>(
64  "c", 1e6, "Parameter for balancing the size of the gap and contact pressure");
65  params.addParam<bool>(
66  "normalize_c",
67  false,
68  "Whether to normalize c by weighting function norm. When unnormalized "
69  "the value of c effectively depends on element size since in the constraint we compare nodal "
70  "Lagrange Multiplier values to integrated gap values (LM nodal value is independent of "
71  "element size, where integrated values are dependent on element size).");
72  params.set<bool>("use_displaced_mesh") = true;
73  params.set<bool>("interpolate_normals") = false;
74  return params;
75 }
76 
78  const InputParameters & parameters)
79  : ADMortarConstraint(assignVarsInParamsDyn(parameters)),
80  _secondary_disp_x(adCoupledValue("disp_x")),
81  _primary_disp_x(adCoupledNeighborValue("disp_x")),
82  _secondary_disp_y(adCoupledValue("disp_y")),
83  _primary_disp_y(adCoupledNeighborValue("disp_y")),
84  _has_disp_z(isCoupled("disp_z")),
85  _secondary_disp_z(_has_disp_z ? &adCoupledValue("disp_z") : nullptr),
86  _primary_disp_z(_has_disp_z ? &adCoupledNeighborValue("disp_z") : nullptr),
87  _c(getParam<Real>("c")),
88  _normalize_c(getParam<bool>("normalize_c")),
89  _nodal(getVar("disp_x", 0)->feType().family == LAGRANGE),
90  _disp_x_var(getVar("disp_x", 0)),
91  _disp_y_var(getVar("disp_y", 0)),
92  _disp_z_var(_has_disp_z ? getVar("disp_z", 0) : nullptr),
93  _capture_tolerance(getParam<Real>("capture_tolerance")),
94  _secondary_x_dot(adCoupledDot("disp_x")),
95  _primary_x_dot(adCoupledNeighborValueDot("disp_x")),
96  _secondary_y_dot(adCoupledDot("disp_y")),
97  _primary_y_dot(adCoupledNeighborValueDot("disp_y")),
98  _secondary_z_dot(_has_disp_z ? &adCoupledDot("disp_z") : nullptr),
99  _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot("disp_z") : nullptr),
100  _has_wear(isParamValid("wear_depth")),
101  _wear_depth(_has_wear ? coupledValueLower("wear_depth") : _zero),
102  _newmark_beta(getParam<Real>("newmark_beta")),
103  _newmark_gamma(getParam<Real>("newmark_gamma"))
104 {
105  if (!useDual())
106  mooseError("Dynamic mortar contact constraints requires the use of Lagrange multipliers dual "
107  "interpolation");
108 }
109 
110 void
112 {
113  // Trim interior node variable derivatives
114  const auto & primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
116  const auto & secondary_ip_lowerd_map =
118 
119  std::array<const MooseVariable *, 3> var_array{{_disp_x_var, _disp_y_var, _disp_z_var}};
120  std::array<ADReal, 3> primary_disp{
121  {_primary_disp_x[_qp], _primary_disp_y[_qp], _has_disp_z ? (*_primary_disp_z)[_qp] : 0}};
122  std::array<ADReal, 3> secondary_disp{{_secondary_disp_x[_qp],
124  _has_disp_z ? (*_secondary_disp_z)[_qp] : 0}};
125 
126  trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp, false);
127  trimInteriorNodeDerivatives(secondary_ip_lowerd_map, var_array, secondary_disp, true);
128 
129  const ADReal & prim_x = primary_disp[0];
130  const ADReal & prim_y = primary_disp[1];
131  const ADReal * prim_z = nullptr;
132  if (_has_disp_z)
133  prim_z = &primary_disp[2];
134 
135  const ADReal & sec_x = secondary_disp[0];
136  const ADReal & sec_y = secondary_disp[1];
137  const ADReal * sec_z = nullptr;
138  if (_has_disp_z)
139  sec_z = &secondary_disp[2];
140 
141  std::array<ADReal, 3> primary_disp_dot{
142  {_primary_x_dot[_qp], _primary_y_dot[_qp], _has_disp_z ? (*_primary_z_dot)[_qp] : 0}};
143  std::array<ADReal, 3> secondary_disp_dot{
144  {_secondary_x_dot[_qp], _secondary_y_dot[_qp], _has_disp_z ? (*_secondary_z_dot)[_qp] : 0}};
145 
146  trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp_dot, false);
147  trimInteriorNodeDerivatives(secondary_ip_lowerd_map, var_array, secondary_disp_dot, true);
148 
149  const ADReal & prim_x_dot = primary_disp_dot[0];
150  const ADReal & prim_y_dot = primary_disp_dot[1];
151  const ADReal * prim_z_dot = nullptr;
152  if (_has_disp_z)
153  prim_z_dot = &primary_disp_dot[2];
154 
155  const ADReal & sec_x_dot = secondary_disp_dot[0];
156  const ADReal & sec_y_dot = secondary_disp_dot[1];
157  const ADReal * sec_z_dot = nullptr;
158  if (_has_disp_z)
159  sec_z_dot = &secondary_disp_dot[2];
160 
161  // Compute dynamic constraint-related quantities
163  gap_vec(0).derivatives() = prim_x.derivatives() - sec_x.derivatives();
164  gap_vec(1).derivatives() = prim_y.derivatives() - sec_y.derivatives();
165 
166  _relative_velocity = ADRealVectorValue(prim_x_dot - sec_x_dot, prim_y_dot - sec_y_dot, 0.0);
167 
168  if (_has_disp_z)
169  {
170  gap_vec(2).derivatives() = prim_z->derivatives() - sec_z->derivatives();
171  _relative_velocity(2) = *prim_z_dot - *sec_z_dot;
172  }
173 
174  _qp_gap_nodal = gap_vec * (_JxW_msm[_qp] * _coord[_qp]);
176 
177  // Current part of the gap velocity Newmark-beta time discretization
179  (_newmark_gamma / _newmark_beta * gap_vec / _dt) * (_JxW_msm[_qp] * _coord[_qp]);
180 
181  // To do normalization of constraint coefficient (c_n)
183 }
184 
185 ADReal
187 {
188  mooseError(
189  "We should never call computeQpResidual for ComputeDynamicWeightedGapLMMechanicalContact");
190 }
191 
192 void
194 {
195  mooseAssert(_normals.size() == _lower_secondary_elem->n_nodes(),
196  "Making sure that _normals is the expected size");
197 
198  // Get the _dof_to_weighted_gap map
199  const DofObject * dof = _var->isNodal()
200  ? static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i))
201  : static_cast<const DofObject *>(_lower_secondary_elem);
202 
203  // Regular normal contact constraint: Use before contact is established for contact detection
204  _dof_to_weighted_gap[dof].first += _test[_i][_qp] * (_qp_gap_nodal * _normals[_i]);
205 
206  // Integrated part of the "persistency" constraint
209 
211 
212  if (_normalize_c)
213  _dof_to_weighted_gap[dof].second += _test[_i][_qp] * _qp_factor;
214 }
215 
216 void
218 {
219  // These dof maps are not recoverable as they are maps of pointers, and recovering old pointers
220  // would be wrong. We would need to create a custom dataStore() and dataLoad()
221  if (_app.isRecovering())
222  mooseError("This object does not support recovering");
223 
224  _dof_to_old_weighted_gap.clear();
225  _dof_to_old_velocity.clear();
227 
228  for (auto & map_pr : _dof_to_weighted_gap)
229  _dof_to_old_weighted_gap.emplace(map_pr.first, std::move(map_pr.second.first));
230 
231  for (auto & map_pr : _dof_to_velocity)
232  _dof_to_old_velocity.emplace(map_pr);
233 
234  for (auto & map_pr : _dof_to_nodal_wear_depth)
235  _dof_to_nodal_old_wear_depth.emplace(map_pr);
236 }
237 
238 void
240 {
241  _dof_to_weighted_gap.clear();
243  _dof_to_velocity.clear();
244 
245  // Wear
246  _dof_to_nodal_wear_depth.clear();
247 }
248 
249 void
251 {
254 
255  if (_has_wear)
256  communicateWear();
257 
258  // There is a need for the dynamic constraint to uncouple the computation of the weighted gap from
259  // the computation of the constraint itself since we are switching from gap constraint to
260  // persistency constraint.
261  for (const auto & pr : _dof_to_weighted_gap)
262  {
263  if (pr.first->processor_id() != this->processor_id())
264  continue;
265 
266  //
267  _dof_to_weighted_gap[pr.first].first += _dof_to_nodal_wear_depth[pr.first];
268  _dof_to_weighted_gap_dynamics[pr.first] +=
270  //
271 
272  const auto is_dof_on_map = _dof_to_old_weighted_gap.find(pr.first);
273 
274  // If is_dof_on_map isn't on map, it means it's an initial step
275  if (is_dof_on_map == _dof_to_old_weighted_gap.end() ||
277  _weighted_gap_ptr = &pr.second.first;
278  else
279  {
281  term += _dof_to_old_velocity[pr.first];
282  _dof_to_weighted_gap_dynamics[pr.first] += term;
284  }
285 
286  _normalization_ptr = &pr.second.second;
287 
289  }
290 }
291 
292 void
294  const std::unordered_set<const Node *> & inactive_lm_nodes)
295 {
298 
299  if (_has_wear)
300  communicateWear();
301 
302  for (const auto & pr : _dof_to_weighted_gap)
303  {
304  if ((inactive_lm_nodes.find(static_cast<const Node *>(pr.first)) != inactive_lm_nodes.end()) ||
305  (pr.first->processor_id() != this->processor_id()))
306  continue;
307 
308  //
309  _dof_to_weighted_gap[pr.first].first += _dof_to_nodal_wear_depth[pr.first];
310  _dof_to_weighted_gap_dynamics[pr.first] +=
312  const auto is_dof_on_map = _dof_to_old_weighted_gap.find(pr.first);
313 
314  // If is_dof_on_map isn't on map, it means it's an initial step
315  if (is_dof_on_map == _dof_to_old_weighted_gap.end() ||
317  {
318  // If this is the first step or the previous step gap is not identified as in contact, apply
319  // regular conditions
320  _weighted_gap_ptr = &pr.second.first;
321  }
322  else
323  {
324  ADReal term = _dof_to_weighted_gap[pr.first].first * _newmark_gamma / (_newmark_beta * _dt);
325  term -= _dof_to_old_weighted_gap[pr.first] * _newmark_gamma / (_newmark_beta * _dt);
326  term -= _dof_to_old_velocity[pr.first];
327  _dof_to_weighted_gap_dynamics[pr.first] = term;
328  // Enable the application of persistency condition
330  }
331 
332  _normalization_ptr = &pr.second.second;
333 
335  }
336 }
337 
338 void
340 {
341  // We may have wear depth information that should go to other processes that own the dofs
342  using Datum = std::pair<dof_id_type, ADReal>;
343  std::unordered_map<processor_id_type, std::vector<Datum>> push_data;
344 
345  for (auto & pr : _dof_to_nodal_wear_depth)
346  {
347  const auto * const dof_object = pr.first;
348  const auto proc_id = dof_object->processor_id();
349  if (proc_id == this->processor_id())
350  continue;
351 
352  push_data[proc_id].push_back(std::make_pair(dof_object->id(), std::move(pr.second)));
353  }
354 
355  const auto & lm_mesh = _mesh.getMesh();
356 
357  auto action_functor = [this, &lm_mesh](const processor_id_type libmesh_dbg_var(pid),
358  const std::vector<Datum> & sent_data)
359  {
360  mooseAssert(pid != this->processor_id(), "We do not send messages to ourself here");
361  for (auto & pr : sent_data)
362  {
363  const auto dof_id = pr.first;
364  const auto * const dof_object =
365  _nodal ? static_cast<const DofObject *>(lm_mesh.node_ptr(dof_id))
366  : static_cast<const DofObject *>(lm_mesh.elem_ptr(dof_id));
367  mooseAssert(dof_object, "This should be non-null");
368  _dof_to_nodal_wear_depth[dof_object] += std::move(pr.second);
369  }
370  };
371 
372  TIMPI::push_parallel_vector_data(_communicator, push_data, action_functor);
373 }
374 
375 void
377 {
378  const auto & weighted_gap = *_weighted_gap_ptr;
379  const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
380 
381  const auto dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
382  ADReal lm_value = (*_sys.currentSolution())(dof_index);
383  Moose::derivInsert(lm_value.derivatives(), dof_index, 1.);
384 
385  const ADReal dof_residual = std::min(lm_value, weighted_gap * c);
386 
388  std::array<ADReal, 1>{{dof_residual}},
389  std::array<dof_id_type, 1>{{dof_index}},
390  _var->scalingFactor());
391 }
392 
393 void
395 {
396  if (mortar_type != Moose::MortarType::Lower)
397  return;
398 
399  mooseAssert(_var, "LM variable is null");
400 
401  for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
402  {
404  for (_i = 0; _i < _test.size(); ++_i)
406  }
407 }
408 
409 void
411 {
412  residualSetup();
413 }
414 
415 void
417 {
418  // During "computeResidual" and "computeJacobian" we are actually just computing properties on the
419  // mortar segment element mesh. We are *not* actually assembling into the residual/Jacobian. For
420  // the zero-penetration constraint, the property of interest is the map from node to weighted gap.
421  // Computation of the properties proceeds identically for residual and Jacobian evaluation hence
422  // why we simply call computeResidual here. We will assemble into the residual/Jacobian later from
423  // the post() method
424  computeResidual(mortar_type);
425 }
MooseMesh & _mesh
virtual void computeQpProperties()
Computes properties that are functions only of the current quadrature point (_qp), e.g.
const VariableValue & _wear_depth
Wear depth to include contact.
LAGRANGE
virtual const NumericVector< Number > *const & currentSolution() const=0
void addRequiredRangeCheckedParam(const std::string &name, const std::string &parsed_function, const std::string &doc_string)
void addResidualsAndJacobian(Assembly &assembly, const Residuals &residuals, const Indices &dof_indices, Real scaling_factor)
ADRealVectorValue _qp_velocity
Vector for computation of weighted gap velocity to fulfill "persistency" condition.
const ADVariableValue & _secondary_disp_x
x-displacement on the secondary face
const ADVariableValue & _primary_disp_y
y-displacement on the primary face
virtual void enforceConstraintOnDof(const DofObject *const dof)
Method called from post().
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
const ADVariableValue & _primary_disp_x
x-displacement on the primary face
void mooseError(Args &&... args)
unsigned int number() const
std::map< unsigned int, unsigned int > getPrimaryIpToLowerElementMap(const Elem &primary_elem, const Elem &primary_elem_ip, const Elem &lower_secondary_elem) const
std::vector< std::pair< R1, R2 > > get(const std::string &param1, const std::string &param2) const
std::unordered_map< const DofObject *, std::pair< ADReal, Real > > _dof_to_weighted_gap
A map from node to weighted gap and normalization (if requested)
std::unordered_map< const DofObject *, ADReal > _dof_to_old_weighted_gap
A map from dof-object to the old weighted gap.
T & set(const std::string &name, bool quiet_mode=false)
const bool _has_disp_z
For 2D mortar contact no displacement will be specified, so const pointers used.
const MooseArray< Point > & _phys_points_primary
Elem const *const & _lower_primary_elem
const std::vector< Real > & _JxW_msm
Computes the normal contact mortar constraints for dynamic simulations.
const Parallel::Communicator & _communicator
unsigned int _i
DualNumber< Real, DNDerivativeType, true > ADReal
virtual void computeResidual() override
const libMesh::QBase *const & _qrule_msm
const VariableTestValue & _test
void suppressParameter(const std::string &name)
std::unordered_map< const DofObject *, ADReal > _dof_to_old_velocity
A map from node to weighted gap velocity times _dt.
std::unordered_map< const DofObject *, ADReal > _dof_to_nodal_wear_depth
A map from node to wear in this step.
SystemBase & _sys
void communicateGaps(std::unordered_map< const DofObject *, std::pair< ADReal, Real >> &dof_to_weighted_gap, const MooseMesh &mesh, bool nodal, bool normalize_c, const Parallel::Communicator &communicator, bool send_data_back)
This function is used to communicate gaps across processes.
void push_parallel_vector_data(const Communicator &comm, MapToVectors &&data, const ActionFunctor &act_on_data)
uint8_t processor_id_type
Elem const *const & _lower_secondary_elem
const AutomaticMortarGeneration & amg() const
MeshBase & getMesh()
static void trimInteriorNodeDerivatives(const std::map< unsigned int, unsigned int > &primary_ip_lowerd_map, const Variables &moose_var, DualNumbers &ad_vars, const bool is_secondary)
std::unordered_map< const DofObject *, ADReal > _dof_to_weighted_gap_dynamics
A map from node to weighted gap velocity times _dt.
std::vector< Point > _normals
virtual void incorrectEdgeDroppingPost(const std::unordered_set< const Node *> &inactive_lm_nodes) override
bool useDual() const
const MooseVariable *const _disp_z_var
The z displacement variable.
const bool _nodal
Whether the dof objects are nodal; if they&#39;re not, then they&#39;re elemental.
registerMooseObject("ContactApp", ComputeDynamicWeightedGapLMMechanicalContact)
const ADVariableValue & _secondary_disp_y
y-displacement on the secondary face
unsigned int n_points() const
std::map< unsigned int, unsigned int > getSecondaryIpToLowerElementMap(const Elem &lower_secondary_elem) const
unsigned int number() const
const MooseArray< Point > & _phys_points_secondary
const ADReal * _weighted_gap_ptr
A pointer members that can be used to help avoid copying ADReals.
Assembly & _assembly
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)
const MooseVariable *const _disp_x_var
The x displacement variable.
ADRealVectorValue _qp_gap_nodal_dynamics
Vector for computation of weighted gap velocity to fulfill "persistency" condition.
MooseVariable *const _var
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool _normalize_c
Whether to normalize weighted gap by weighting function norm.
const Real _capture_tolerance
A small threshold gap value to consider that a node needs a "persistency" constraint.
std::unordered_map< const DofObject *, ADReal > _dof_to_nodal_old_wear_depth
A map from node to wear in old step.
virtual void computeJacobian() override
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
static InputParameters validParams()
void derivInsert(SemiDynamicSparseNumberArray< Real, libMesh::dof_id_type, NWrapper< N >> &derivs, libMesh::dof_id_type index, Real value)
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
std::unordered_map< const DofObject *, ADReal > _dof_to_velocity
A map from node to weighted gap velocity times _dt.
const bool _has_wear
Flag to determine whether wear needs to be included in the contact constraints.
const MooseArray< Real > & _coord
processor_id_type processor_id() const
Real _qp_factor
The value of the LM at the current quadrature point.
bool isRecovering() const
ADRealVectorValue _qp_gap_nodal
Vector for computation of weighted gap with nodal normals.
unsigned int _qp
const Real _c
This factor multiplies the weighted gap.
void scalingFactor(const std::vector< Real > &factor)
const MooseVariable *const _disp_y_var
The y displacement variable.