https://mooseframework.inl.gov
ComputeDynamicFrictionalForceLMMechanicalContact.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 "Function.h"
14 #include "MortarContactUtils.h"
15 
16 #include "metaphysicl/dualsemidynamicsparsenumberarray.h"
17 #include "metaphysicl/parallel_dualnumber.h"
18 #include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
19 #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
20 #include "timpi/parallel_sync.h"
21 
22 #include <limits>
23 
25 
28 {
30  params.addClassDescription("Computes the tangential frictional forces for dynamic simulations");
31  params.addRequiredCoupledVar("friction_lm", "The frictional Lagrange's multiplier");
32  params.addCoupledVar("friction_lm_dir",
33  "The frictional Lagrange's multiplier for an addtional direction.");
34  params.addParam<FunctionName>(
35  "function_friction",
36  "Coupled function to evaluate friction with values from contact pressure and relative "
37  "tangential velocities (from the previous step).");
38  params.addParam<Real>("c_t", 1e0, "Numerical parameter for tangential constraints");
39  params.addParam<Real>(
40  "epsilon",
41  1.0e-7,
42  "Minimum value of contact pressure that will trigger frictional enforcement");
43  params.addParam<Real>("mu", "The friction coefficient for the Coulomb friction law");
44  return params;
45 }
46 
48  const InputParameters & parameters)
50  _c_t(getParam<Real>("c_t")),
51  _secondary_x_dot(_secondary_var.adUDot()),
52  _primary_x_dot(_primary_var.adUDotNeighbor()),
53  _secondary_y_dot(adCoupledDot("disp_y")),
54  _primary_y_dot(adCoupledNeighborValueDot("disp_y")),
55  _secondary_z_dot(_has_disp_z ? &adCoupledDot("disp_z") : nullptr),
56  _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot("disp_z") : nullptr),
57  _epsilon(getParam<Real>("epsilon")),
58  _mu(isParamValid("mu") ? getParam<Real>("mu") : std::numeric_limits<double>::quiet_NaN()),
59  _function_friction(isParamValid("function_friction") ? &getFunction("function_friction")
60  : nullptr),
61  _has_friction_function(isParamValid("function_friction")),
62  _3d(_has_disp_z)
63 {
64  if (!_has_friction_function && !isParamValid("mu"))
65  mooseError(
66  "A coefficient of friction needs to be provided as a constant value of via a function.");
67 
69  paramError("mu",
70  "Either provide a constant coefficient of friction or a function defining the "
71  "coefficient of friction. Both inputs cannot be provided simultaneously.");
72 
73  if (!getParam<bool>("use_displaced_mesh"))
74  paramError("use_displaced_mesh",
75  "'use_displaced_mesh' must be true for the "
76  "ComputeFrictionalForceLMMechanicalContact object");
77 
78  if (_3d && !isParamValid("friction_lm_dir"))
79  paramError("friction_lm_dir",
80  "Three-dimensional mortar frictional contact simulations require an additional "
81  "frictional Lagrange's multiplier to enforce a second tangential pressure");
82 
83  _friction_vars.push_back(getVar("friction_lm", 0));
84 
85  if (_3d)
86  _friction_vars.push_back(getVar("friction_lm_dir", 0));
87 
88  if (!_friction_vars[0]->isNodal())
89  if (_friction_vars[0]->feType().order != static_cast<Order>(0))
90  paramError(
91  "friction_lm",
92  "Frictional contact constraints only support elemental variables of CONSTANT order");
93 
94  // Request the old solution state in unison
95  _sys.solutionOld();
96 }
97 
98 void
100 {
101  // Compute the value of _qp_gap
103 
104  // It appears that the relative velocity between weighted gap and this class have a sign
105  // difference
108 }
109 
110 void
112 {
113  // Get the _dof_to_weighted_gap map
115 
116  const auto & nodal_tangents = amg().getNodalTangents(*_lower_secondary_elem);
117 
118  // Get the _dof_to_weighted_tangential_velocity map
119  const DofObject * const dof =
120  _friction_vars[0]->isNodal()
121  ? static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i))
122  : static_cast<const DofObject *>(_lower_secondary_elem);
123 
125  _test[_i][_qp] * _qp_tangential_velocity_nodal * nodal_tangents[0][_i];
126 
129  nodal_tangents[0][_i];
130 
131  // Get the _dof_to_weighted_tangential_velocity map for a second direction
132  if (_3d)
133  {
135  _test[_i][_qp] * _qp_tangential_velocity_nodal * nodal_tangents[1][_i];
136 
139  nodal_tangents[1][_i];
140  }
141 }
142 
143 void
145 {
146  // Clear both maps
150 }
151 
152 void
154 {
155 
157 
159 
160  for (auto & map_pr : _dof_to_real_tangential_velocity)
161  _dof_to_old_real_tangential_velocity.emplace(map_pr);
162 }
163 
164 void
166 {
168 
171 
175 
176  // Enforce frictional complementarity constraints
177  for (const auto & pr : _dof_to_weighted_tangential_velocity)
178  {
179  const DofObject * const dof = pr.first;
180 
181  if (dof->processor_id() != this->processor_id())
182  continue;
183 
184  // Use always weighted gap for dynamic PDASS. Omit the dynamic weighted gap approach that is
185  // used in normal contact where the discretized gap velocity is enforced if a node has
186  // identified to be into contact.
187 
190  _tangential_vel_ptr[0] = &(pr.second[0]);
191 
192  if (_3d)
193  {
194  _tangential_vel_ptr[1] = &(pr.second[1]);
196  }
197  else
199  }
200 }
201 
202 void
204  const std::unordered_set<const Node *> & inactive_lm_nodes)
205 {
207 
210 
214 
215  // Enforce frictional complementarity constraints
216  for (const auto & pr : _dof_to_weighted_tangential_velocity)
217  {
218  const DofObject * const dof = pr.first;
219 
220  // If node inactive, skip
221  if ((inactive_lm_nodes.find(static_cast<const Node *>(dof)) != inactive_lm_nodes.end()) ||
222  (dof->processor_id() != this->processor_id()))
223  continue;
224 
225  // Use always weighted gap for dynamic PDASS
228  _tangential_vel_ptr[0] = &pr.second[0];
229 
230  if (_3d)
231  {
232  _tangential_vel_ptr[1] = &pr.second[1];
234  }
235  else
237  }
238 }
239 
240 void
242  const DofObject * const dof)
243 {
244  // Get normal LM
245  const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
246  const ADReal & weighted_gap = *_weighted_gap_ptr;
247  ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
248  Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
249 
250  // Get friction LMs
251  std::array<const ADReal *, 2> & tangential_vel = _tangential_vel_ptr;
252  std::array<dof_id_type, 2> friction_dof_indices;
253  std::array<ADReal, 2> friction_lm_values;
254 
255  const unsigned int num_tangents = 2;
256  for (const auto i : make_range(num_tangents))
257  {
258  friction_dof_indices[i] = dof->dof_number(_sys.number(), _friction_vars[i]->number(), 0);
259  friction_lm_values[i] = (*_sys.currentSolution())(friction_dof_indices[i]);
260  Moose::derivInsert(friction_lm_values[i].derivatives(), friction_dof_indices[i], 1.);
261  }
262 
263  // Get normalized c and c_t values (if normalization specified
264  const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
265  const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
266 
267  const Real contact_pressure_old = _sys.solutionOld()(normal_dof_index);
268 
269  // Compute the friction coefficient (constant or function)
270  ADReal mu_ad = computeFrictionValue(contact_pressure_old,
273 
274  ADReal dof_residual;
275  ADReal dof_residual_dir;
276 
277  // Primal-dual active set strategy (PDASS)
278  if (contact_pressure < _epsilon)
279  {
280  dof_residual = friction_lm_values[0];
281  dof_residual_dir = friction_lm_values[1];
282  }
283  else
284  {
285  const Real epsilon_sqrt = 1.0e-48;
286 
287  const auto lamdba_plus_cg = contact_pressure + c * weighted_gap;
288  std::array<ADReal, 2> lambda_t_plus_ctu;
289  lambda_t_plus_ctu[0] = friction_lm_values[0] + c_t * *tangential_vel[0] * _dt;
290  lambda_t_plus_ctu[1] = friction_lm_values[1] + c_t * *tangential_vel[1] * _dt;
291 
292  const auto term_1_x =
293  std::max(mu_ad * lamdba_plus_cg,
294  std::sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
295  lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
296  friction_lm_values[0];
297 
298  const auto term_1_y =
299  std::max(mu_ad * lamdba_plus_cg,
300  std::sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
301  lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
302  friction_lm_values[1];
303 
304  const auto term_2_x = mu_ad * std::max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[0];
305 
306  const auto term_2_y = mu_ad * std::max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[1];
307 
308  dof_residual = term_1_x - term_2_x;
309  dof_residual_dir = term_1_y - term_2_y;
310  }
311 
313  std::array<ADReal, 1>{{dof_residual}},
314  std::array<dof_id_type, 1>{{friction_dof_indices[0]}},
315  _friction_vars[0]->scalingFactor());
317  std::array<ADReal, 1>{{dof_residual_dir}},
318  std::array<dof_id_type, 1>{{friction_dof_indices[1]}},
319  _friction_vars[1]->scalingFactor());
320 }
321 
322 void
324  const DofObject * const dof)
325 {
326  // Get friction LM
327  const auto friction_dof_index = dof->dof_number(_sys.number(), _friction_vars[0]->number(), 0);
328  const ADReal & tangential_vel = *_tangential_vel_ptr[0];
329  ADReal friction_lm_value = (*_sys.currentSolution())(friction_dof_index);
330  Moose::derivInsert(friction_lm_value.derivatives(), friction_dof_index, 1.);
331 
332  // Get normal LM
333  const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
334  const ADReal & weighted_gap = *_weighted_gap_ptr;
335  ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
336  Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
337 
338  const Real contact_pressure_old = _sys.solutionOld()(normal_dof_index);
339 
340  // Get normalized c and c_t values (if normalization specified
341  const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
342  const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
343 
344  // Compute the friction coefficient (constant or function)
345  ADReal mu_ad =
346  computeFrictionValue(contact_pressure_old, _dof_to_old_real_tangential_velocity[dof][0], 0.0);
347 
348  ADReal dof_residual;
349  // Primal-dual active set strategy (PDASS)
350  if (contact_pressure < _epsilon)
351  dof_residual = friction_lm_value;
352  else
353  {
354  const auto term_1 = std::max(mu_ad * (contact_pressure + c * weighted_gap),
355  std::abs(friction_lm_value + c_t * tangential_vel * _dt)) *
356  friction_lm_value;
357  const auto term_2 = mu_ad * std::max(0.0, contact_pressure + c * weighted_gap) *
358  (friction_lm_value + c_t * tangential_vel * _dt);
359 
360  dof_residual = term_1 - term_2;
361  }
362 
364  std::array<ADReal, 1>{{dof_residual}},
365  std::array<dof_id_type, 1>{{friction_dof_index}},
366  _friction_vars[0]->scalingFactor());
367 }
368 
369 ADReal
371  const ADReal & contact_pressure, const Real & tangential_vel, const Real & tangential_vel_dir)
372 {
373  // TODO: Introduce temperature dependence in the function. Do this when we have an example.
374  ADReal mu_ad;
375 
377  mu_ad = _mu;
378  else
379  {
380  ADReal tangential_vel_magnitude = std::sqrt(tangential_vel * tangential_vel +
381  tangential_vel_dir * tangential_vel_dir + 1.0e-24);
382 
383  mu_ad = _function_friction->value<ADReal>(0.0, contact_pressure, tangential_vel_magnitude, 0.0);
384  }
385 
386  return mu_ad;
387 }
MooseMesh & _mesh
virtual void computeQpProperties() override
Computes properties that are functions only of the current quadrature point (_qp), e.g.
virtual void computeQpProperties()
Computes properties that are functions only of the current quadrature point (_qp), e.g.
virtual const NumericVector< Number > *const & currentSolution() const=0
ADRealVectorValue _qp_tangential_velocity_nodal
The value of the tangential velocity vectors at the current node.
ADRealVectorValue _qp_real_tangential_velocity_nodal
The value of the tangential velocity vectors at the current node.
void addResidualsAndJacobian(Assembly &assembly, const Residuals &residuals, const Indices &dof_indices, Real scaling_factor)
std::unordered_map< const DofObject *, std::array< Real, 2 > > _dof_to_real_tangential_velocity
A map from node to two tangential velocities. Required to have direct connection to physics...
void communicateVelocities(std::unordered_map< const DofObject *, T > &dof_map, const MooseMesh &mesh, const bool nodal, const Parallel::Communicator &communicator, const bool send_data_back)
This function is used to communicate velocities across processes.
const Real _epsilon
Small contact pressure value to trigger computation of frictional forces.
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
unsigned int number() const
registerMooseObject("ContactApp", ComputeDynamicFrictionalForceLMMechanicalContact)
std::unordered_map< const DofObject *, std::pair< ADReal, Real > > _dof_to_weighted_gap
A map from node to weighted gap and normalization (if requested)
ADReal computeFrictionValue(const ADReal &contact_pressure, const Real &tangential_vel, const Real &tangential_vel_dir)
Apply constant or function-based friction coefficient.
auto raw_value(const Eigen::Map< T > &in)
MooseVariable * getVar(const std::string &var_name, unsigned int comp)
std::vector< MooseVariable * > _friction_vars
Frictional Lagrange&#39;s multiplier variable pointers.
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
void incorrectEdgeDroppingPost(const std::unordered_set< const Node *> &inactive_lm_nodes) override
Copy of the post routine but that skips assembling inactive nodes.
const VariableTestValue & _test
bool isParamValid(const std::string &name) const
SystemBase & _sys
std::unordered_map< const DofObject *, std::array< ADReal, 2 > > _dof_to_weighted_tangential_velocity
A map from node to two weighted tangential velocities.
Elem const *const & _lower_secondary_elem
const AutomaticMortarGeneration & amg() const
virtual void incorrectEdgeDroppingPost(const std::unordered_set< const Node *> &inactive_lm_nodes) override
virtual void enforceConstraintOnDof3d(const DofObject *const dof)
Method called from post().
const bool _nodal
Whether the dof objects are nodal; if they&#39;re not, then they&#39;re elemental.
std::array< MooseUtils::SemidynamicVector< Point, 9 >, 2 > getNodalTangents(const Elem &secondary_elem) const
void paramError(const std::string &param, Args... args) const
unsigned int number() const
virtual void enforceConstraintOnDof(const DofObject *const dof) override
Method called from post().
const ADReal * _weighted_gap_ptr
A pointer members that can be used to help avoid copying ADReals.
Assembly & _assembly
bool isNodal() const
void addCoupledVar(const std::string &name, const std::string &doc_string)
void addRequiredCoupledVar(const std::string &name, const std::string &doc_string)
MooseVariable *const _var
std::unordered_map< const DofObject *, std::array< Real, 2 > > _dof_to_old_real_tangential_velocity
A map from node to two old tangential velocities. Required to have direct connection to physics...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool _normalize_c
Whether to normalize weighted gap by weighting function norm.
IntRange< T > make_range(T beg, T end)
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
void derivInsert(SemiDynamicSparseNumberArray< Real, libMesh::dof_id_type, NWrapper< N >> &derivs, libMesh::dof_id_type index, Real value)
std::array< const ADReal *, 2 > _tangential_vel_ptr
An array of two pointers to avoid copies.
const bool _has_friction_function
Boolean to determine whether the friction coefficient is taken from a function.
virtual Real value(Real t, const Point &p) const
const MooseArray< Real > & _coord
NumericVector< Number > & solutionOld()
processor_id_type processor_id() const
Computes the mortar tangential frictional forces for dynamic simulations.
unsigned int _qp
bool _3d
Automatic flag to determine whether we are doing three-dimensional work.
virtual void computeQpIProperties() override
Computes properties that are functions both of _qp and _i, for example the weighted gap...
const Real _c
This factor multiplies the weighted gap.
const Real _c_t
Numerical factor used in the tangential constraints for convergence purposes.