https://mooseframework.inl.gov
ComputeFrictionalForceLMMechanicalContact.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"
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 <cmath>
23 
25 
28 {
30  params.addClassDescription("Computes the tangential frictional forces");
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");
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.addRangeCheckedParam<Real>(
44  "mu", "mu > 0", "The friction coefficient for the Coulomb friction law");
45  params.addRequiredParam<UserObjectName>("weighted_velocities_uo",
46  "The weighted tangential velocities user object.");
47 
48  return params;
49 }
50 
52  const InputParameters & parameters)
54  _weighted_velocities_uo(getUserObject<WeightedVelocitiesUserObject>("weighted_velocities_uo")),
55  _c_t(getParam<Real>("c_t")),
56  _secondary_x_dot(_secondary_var.adUDot()),
57  _primary_x_dot(_primary_var.adUDotNeighbor()),
58  _secondary_y_dot(adCoupledDot("disp_y")),
59  _primary_y_dot(adCoupledNeighborValueDot("disp_y")),
60  _secondary_z_dot(_has_disp_z ? &adCoupledDot("disp_z") : nullptr),
61  _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot("disp_z") : nullptr),
62  _epsilon(getParam<Real>("epsilon")),
63  _mu(isParamValid("function_friction") ? std::numeric_limits<double>::quiet_NaN()
64  : getParam<Real>("mu")),
65  _function_friction(isParamValid("function_friction") ? &getFunction("function_friction")
66  : nullptr),
67  _has_friction_function(isParamValid("function_friction")),
68  _3d(_has_disp_z)
69 
70 {
72  paramError(
73  "mu",
74  "Please only provide friction either as a function or as a constant value, but not both.");
76  paramError("mu", "Please provide a value or a function for the coefficient of friction.");
77 
78  if (!getParam<bool>("use_displaced_mesh"))
79  paramError("use_displaced_mesh",
80  "'use_displaced_mesh' must be true for the "
81  "ComputeFrictionalForceLMMechanicalContact object");
82 
83  if (_3d && !isParamValid("friction_lm_dir"))
84  paramError("friction_lm_dir",
85  "Three-dimensional mortar frictional contact simulations require an additional "
86  "frictional Lagrange's multiplier to enforce a second tangential pressure");
87 
88  _friction_vars.push_back(getVar("friction_lm", 0));
89 
90  if (_3d)
91  _friction_vars.push_back(getVar("friction_lm_dir", 0));
92 
93  if (!_friction_vars[0]->isNodal())
94  if (_friction_vars[0]->feType().order != static_cast<Order>(0))
95  paramError(
96  "friction_lm",
97  "Frictional contact constraints only support elemental variables of CONSTANT order");
98 }
99 
100 void
102 {
103 }
104 
105 void
107 {
108 }
109 
110 void
112 {
113 }
114 
115 void
117 {
118  const auto & dof_to_weighted_tangential_velocity =
120 
121  const std::unordered_map<const DofObject *, std::pair<ADReal, Real>> & dof_to_weighted_gap =
123 
124  // Enforce frictional constraints
125 
126  for (const auto & [dof_object, weighted_velocities_pr] : dof_to_weighted_tangential_velocity)
127  {
128  if (dof_object->processor_id() != this->processor_id())
129  continue;
130 
131  const auto & [weighted_gap_pr, normalization] =
132  libmesh_map_find(dof_to_weighted_gap, dof_object);
133  _weighted_gap_ptr = &weighted_gap_pr;
134  _normalization_ptr = &normalization;
135  _tangential_vel_ptr[0] = &(weighted_velocities_pr[0]);
136 
137  if (_3d)
138  {
139  _tangential_vel_ptr[1] = &(weighted_velocities_pr[1]);
140  enforceConstraintOnDof3d(dof_object);
141  }
142  else
143  enforceConstraintOnDof(dof_object);
144  }
145 }
146 
147 void
149  const std::unordered_set<const Node *> & inactive_lm_nodes)
150 {
151  const auto & dof_to_weighted_tangential_velocity =
153  const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
154  // Enforce frictional complementarity constraints
155  for (const auto & [dof_object, weighted_velocities_pr] : dof_to_weighted_tangential_velocity)
156  {
157  // If node inactive, skip
158  if ((inactive_lm_nodes.find(static_cast<const Node *>(dof_object)) !=
159  inactive_lm_nodes.end()) ||
160  (dof_object->processor_id() != this->processor_id()))
161  continue;
162 
163  _weighted_gap_ptr = &dof_to_weighted_gap.at(dof_object).first;
164  _normalization_ptr = &dof_to_weighted_gap.at(dof_object).second;
165  _tangential_vel_ptr[0] = &weighted_velocities_pr[0];
166 
167  if (_3d)
168  {
169  _tangential_vel_ptr[1] = &weighted_velocities_pr[1];
170  enforceConstraintOnDof3d(dof_object);
171  }
172  else
173  enforceConstraintOnDof(dof_object);
174  }
175 }
176 
177 void
179 {
181 
182  // Get normal LM
183  const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
184  const ADReal & weighted_gap = *_weighted_gap_ptr;
185  ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
186  Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
187 
188  // Get friction LMs
189  std::array<const ADReal *, 2> & tangential_vel = _tangential_vel_ptr;
190  std::array<dof_id_type, 2> friction_dof_indices;
191  std::array<ADReal, 2> friction_lm_values;
192 
193  const unsigned int num_tangents = 2;
194  for (const auto i : make_range(num_tangents))
195  {
196  friction_dof_indices[i] = dof->dof_number(_sys.number(), _friction_vars[i]->number(), 0);
197  friction_lm_values[i] = (*_sys.currentSolution())(friction_dof_indices[i]);
198  Moose::derivInsert(friction_lm_values[i].derivatives(), friction_dof_indices[i], 1.);
199  }
200 
201  // Get normalized c and c_t values (if normalization specified
202  const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
203  const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
204 
205  // Compute the friction coefficient (constant or function)
206  ADReal mu_ad = computeFrictionValue(contact_pressure,
209 
210  ADReal dof_residual;
211  ADReal dof_residual_dir;
212 
213  // Primal-dual active set strategy (PDASS)
214  if (contact_pressure < _epsilon)
215  {
216  dof_residual = friction_lm_values[0];
217  dof_residual_dir = friction_lm_values[1];
218  }
219  else
220  {
221  // Espilon to avoid automatic differentiation singularity
222  const Real epsilon_sqrt = 1.0e-48;
223 
224  const auto lamdba_plus_cg = contact_pressure + c * weighted_gap;
225  std::array<ADReal, 2> lambda_t_plus_ctu;
226  lambda_t_plus_ctu[0] = friction_lm_values[0] + c_t * *tangential_vel[0] * _dt;
227  lambda_t_plus_ctu[1] = friction_lm_values[1] + c_t * *tangential_vel[1] * _dt;
228 
229  const auto term_1_x =
230  std::max(mu_ad * lamdba_plus_cg,
231  std::sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
232  lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
233  friction_lm_values[0];
234 
235  const auto term_1_y =
236  std::max(mu_ad * lamdba_plus_cg,
237  std::sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
238  lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
239  friction_lm_values[1];
240 
241  const auto term_2_x = mu_ad * std::max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[0];
242 
243  const auto term_2_y = mu_ad * std::max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[1];
244 
245  dof_residual = term_1_x - term_2_x;
246  dof_residual_dir = term_1_y - term_2_y;
247  }
248 
250  std::array<ADReal, 1>{{dof_residual}},
251  std::array<dof_id_type, 1>{{friction_dof_indices[0]}},
252  _friction_vars[0]->scalingFactor());
254  std::array<ADReal, 1>{{dof_residual_dir}},
255  std::array<dof_id_type, 1>{{friction_dof_indices[1]}},
256  _friction_vars[1]->scalingFactor());
257 }
258 
259 void
261 {
263 
264  // Get friction LM
265  const auto friction_dof_index = dof->dof_number(_sys.number(), _friction_vars[0]->number(), 0);
266  const ADReal & tangential_vel = *_tangential_vel_ptr[0];
267  ADReal friction_lm_value = (*_sys.currentSolution())(friction_dof_index);
268  Moose::derivInsert(friction_lm_value.derivatives(), friction_dof_index, 1.);
269 
270  // Get normal LM
271  const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
272  const ADReal & weighted_gap = *_weighted_gap_ptr;
273  ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
274  Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
275 
276  // Get normalized c and c_t values (if normalization specified
277  const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
278  const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
279 
280  // Compute the friction coefficient (constant or function)
281  ADReal mu_ad =
282  computeFrictionValue(contact_pressure, _dof_to_real_tangential_velocity[dof][0], 0.0);
283 
284  ADReal dof_residual;
285  // Primal-dual active set strategy (PDASS)
286  if (contact_pressure < _epsilon)
287  dof_residual = friction_lm_value;
288  else
289  {
290  const auto term_1 = std::max(mu_ad * (contact_pressure + c * weighted_gap),
291  std::abs(friction_lm_value + c_t * tangential_vel * _dt)) *
292  friction_lm_value;
293  const auto term_2 = mu_ad * std::max(0.0, contact_pressure + c * weighted_gap) *
294  (friction_lm_value + c_t * tangential_vel * _dt);
295 
296  dof_residual = term_1 - term_2;
297  }
298 
300  std::array<ADReal, 1>{{dof_residual}},
301  std::array<dof_id_type, 1>{{friction_dof_index}},
302  _friction_vars[0]->scalingFactor());
303 }
304 
305 ADReal
307  const ADReal & tangential_vel,
308  const ADReal & tangential_vel_dir)
309 {
310  // TODO: Introduce temperature dependence in the function. Do this when we have an example.
311  ADReal mu_ad;
312 
314  mu_ad = _mu;
315  else
316  {
317  ADReal tangential_vel_magnitude = std::sqrt(tangential_vel * tangential_vel +
318  tangential_vel_dir * tangential_vel_dir + 1.0e-24);
319  mu_ad = _function_friction->value<ADReal>(0.0, contact_pressure, tangential_vel_magnitude, 0.0);
320  }
321 
322  return mu_ad;
323 }
ComputeFrictionalForceLMMechanicalContact(const InputParameters &parameters)
virtual void enforceConstraintOnDof(const DofObject *const dof) override
Method called from post().
virtual const NumericVector< Number > *const & currentSolution() const=0
void addResidualsAndJacobian(Assembly &assembly, const Residuals &residuals, const Indices &dof_indices, Real scaling_factor)
std::unordered_map< const DofObject *, std::array< ADReal, 2 > > _dof_to_real_tangential_velocity
A map from node to two interpolated, physical tangential velocities.
const WeightedVelocitiesUserObject & _weighted_velocities_uo
The weighted gap user object.
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
unsigned int number() const
const bool _has_friction_function
Boolean to determine whether the friction coefficient is taken from a function.
const Real _c_t
Numerical factor used in the tangential constraints for convergence purposes.
virtual void computeQpIProperties() override
Computes properties that are functions both of _qp and _i, for example the weighted gap...
MooseVariable * getVar(const std::string &var_name, unsigned int comp)
Creates dof object to weighted tangential velocities map.
Computes frictional constraints (and normal contact constraints by calling its parent object) ...
bool _3d
Automatic flag to determine whether we are doing three-dimensional work.
const Real _c
This factor multiplies the weighted gap.
const WeightedGapUserObject & _weighted_gap_uo
The weighted gap user object.
DualNumber< Real, DNDerivativeType, true > ADReal
void addRequiredParam(const std::string &name, const std::string &doc_string)
bool isParamValid(const std::string &name) const
SystemBase & _sys
const std::unordered_map< const DofObject *, std::array< ADReal, 2 > > & dofToWeightedVelocities() const
Get the degree of freedom to weighted velocities information.
std::vector< MooseVariable * > _friction_vars
Frictional Lagrange&#39;s multiplier variable pointers.
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.
const Real _epsilon
Small contact pressure value to trigger computation of frictional forces.
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)
bool isParamSetByUser(const std::string &name) const
MooseVariable *const _var
const bool _normalize_c
Whether to normalize weighted gap by weighting function norm.
ADReal computeFrictionValue(const ADReal &contact_pressure, const ADReal &tangential_vel, const ADReal &tangential_vel_dir)
Compute coefficient of friction.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)
registerMooseObject("ContactApp", ComputeFrictionalForceLMMechanicalContact)
void addClassDescription(const std::string &doc_string)
const InputParameters & parameters() const
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)
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
virtual void computeQpProperties() override
Computes properties that are functions only of the current quadrature point (_qp), e.g.
virtual Real value(Real t, const Point &p) const
processor_id_type processor_id() const
void incorrectEdgeDroppingPost(const std::unordered_set< const Node *> &inactive_lm_nodes) override
Copy of the post routine but that skips assembling inactive nodes.
Computes the weighted gap that will later be used to enforce the zero-penetration mechanical contact ...
std::array< const ADReal *, 2 > _tangential_vel_ptr
An array of two pointers to avoid copies.
virtual void enforceConstraintOnDof3d(const DofObject *const dof)
Method called from post().