https://mooseframework.inl.gov
ComputeFrictionalForceCartesianLMMechanicalContact.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 #include "metaphysicl/dualsemidynamicsparsenumberarray.h"
15 #include "metaphysicl/parallel_dualnumber.h"
16 #include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
17 #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
18 #include "timpi/parallel_sync.h"
19 
21 
22 namespace
23 {
24 const InputParameters &
25 assignVarsInParamsFriction(const InputParameters & params_in)
26 {
27  InputParameters & ret = const_cast<InputParameters &>(params_in);
28  const auto & disp_x_name = ret.get<std::vector<VariableName>>("disp_x");
29  if (disp_x_name.size() != 1)
30  mooseError("We require that the disp_x parameter have exactly one coupled name");
31 
32  // We do this so we don't get any variable errors during MortarConstraint(Base) construction
33  ret.set<VariableName>("secondary_variable") = disp_x_name[0];
34  ret.set<VariableName>("primary_variable") = disp_x_name[0];
35 
36  return ret;
37 }
38 }
39 
42 {
44  params.addClassDescription("Computes mortar frictional forces.");
45  params.addParam<Real>("c_t", 1e0, "Numerical parameter for tangential constraints");
46  params.addParam<Real>(
47  "epsilon",
48  1.0e-7,
49  "Minimum value of contact pressure that will trigger frictional enforcement");
50  params.addRangeCheckedParam<Real>(
51  "mu", "mu > 0", "The friction coefficient for the Coulomb friction law");
52  return params;
53 }
54 
57  : ComputeWeightedGapCartesianLMMechanicalContact(assignVarsInParamsFriction(parameters)),
58  _c_t(getParam<Real>("c_t")),
59  _secondary_x_dot(adCoupledDot("disp_x")),
60  _primary_x_dot(adCoupledNeighborValueDot("disp_x")),
61  _secondary_y_dot(adCoupledDot("disp_y")),
62  _primary_y_dot(adCoupledNeighborValueDot("disp_y")),
63  _secondary_z_dot(_has_disp_z ? &adCoupledDot("disp_z") : nullptr),
64  _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot("disp_z") : nullptr),
65  _mu(getParam<Real>("mu")),
66  _epsilon(getParam<Real>("epsilon"))
67 {
68 }
69 
70 void
72 {
74 
75  // Trim derivatives
76  const auto & primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
78  const auto & secondary_ip_lowerd_map =
80 
81  std::array<const MooseVariable *, 3> var_array{{_disp_x_var, _disp_y_var, _disp_z_var}};
82  std::array<ADReal, 3> primary_disp_dot{
83  {_primary_x_dot[_qp], _primary_y_dot[_qp], _has_disp_z ? (*_primary_z_dot)[_qp] : 0}};
84  std::array<ADReal, 3> secondary_disp_dot{
85  {_secondary_x_dot[_qp], _secondary_y_dot[_qp], _has_disp_z ? (*_secondary_z_dot)[_qp] : 0}};
86 
87  trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp_dot, false);
88  trimInteriorNodeDerivatives(secondary_ip_lowerd_map, var_array, secondary_disp_dot, true);
89 
90  const ADReal & prim_x_dot = primary_disp_dot[0];
91  const ADReal & prim_y_dot = primary_disp_dot[1];
92  const ADReal * prim_z_dot = nullptr;
93  if (_has_disp_z)
94  prim_z_dot = &primary_disp_dot[2];
95 
96  const ADReal & sec_x_dot = secondary_disp_dot[0];
97  const ADReal & sec_y_dot = secondary_disp_dot[1];
98  const ADReal * sec_z_dot = nullptr;
99  if (_has_disp_z)
100  sec_z_dot = &secondary_disp_dot[2];
101 
102  // Build relative velocity vector
103  ADRealVectorValue relative_velocity;
104 
105  if (_has_disp_z)
106  relative_velocity = {sec_x_dot - prim_x_dot, sec_y_dot - prim_y_dot, *sec_z_dot - *prim_z_dot};
107  else
108  relative_velocity = {sec_x_dot - prim_x_dot, sec_y_dot - prim_y_dot, 0.0};
109 
110  _qp_tangential_velocity_nodal = relative_velocity * (_JxW_msm[_qp] * _coord[_qp]);
111 }
112 
113 void
115 {
117 
118  const auto & nodal_tangents = amg().getNodalTangents(*_lower_secondary_elem);
119  // Get the _dof_to_weighted_tangential_velocity map
120  const DofObject * const dof =
121  _lm_vars[0]->isNodal() ? 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 
127  // Get the _dof_to_weighted_tangential_velocity map for a second direction
128  if (_has_disp_z)
130  _test[_i][_qp] * _qp_tangential_velocity_nodal * nodal_tangents[1][_i];
131 }
132 
133 void
135 {
138 }
139 
140 void
142 {
147 
148  // Enforce frictional complementarity constraints
149  for (const auto & pr : _dof_to_weighted_tangential_velocity)
150  {
151  const DofObject * const dof = pr.first;
152 
153  if (dof->processor_id() != this->processor_id())
154  continue;
155 
156  auto & weighted_gap_pr = _dof_to_weighted_gap[dof];
157  _weighted_gap_ptr = &weighted_gap_pr.first;
158  _normalization_ptr = &weighted_gap_pr.second;
159  _tangential_vel_ptr[0] = &(pr.second[0]);
160 
161  if (_has_disp_z)
162  _tangential_vel_ptr[1] = &(pr.second[1]);
163 
165  }
166 }
167 
168 void
170  const std::unordered_set<const Node *> & inactive_lm_nodes)
171 {
176 
177  // Enforce frictional complementarity constraints
178  for (const auto & pr : _dof_to_weighted_tangential_velocity)
179  {
180  const DofObject * const dof = pr.first;
181 
182  // If node inactive, skip
183  if ((inactive_lm_nodes.find(static_cast<const Node *>(dof)) != inactive_lm_nodes.end()) ||
184  (dof->processor_id() != this->processor_id()))
185  continue;
186 
189  _tangential_vel_ptr[0] = &pr.second[0];
190 
191  if (_has_disp_z)
192  _tangential_vel_ptr[1] = &pr.second[1];
193 
195  }
196 }
197 
198 void
200  const DofObject * const dof)
201 {
202  const auto & weighted_gap = *_weighted_gap_ptr;
203  const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
204  const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
205 
206  const auto dof_index_x = dof->dof_number(_sys.number(), _lm_vars[0]->number(), 0);
207  const auto dof_index_y = dof->dof_number(_sys.number(), _lm_vars[1]->number(), 0);
208  const Real scaling_factor_x = _lm_vars[0]->scalingFactor();
209  const Real scaling_factor_y = _lm_vars[1]->scalingFactor();
210  Real scaling_factor_z = 1;
211 
212  ADReal lm_x = (*_sys.currentSolution())(dof_index_x);
213  ADReal lm_y = (*_sys.currentSolution())(dof_index_y);
214 
215  Moose::derivInsert(lm_x.derivatives(), dof_index_x, 1.);
216  Moose::derivInsert(lm_y.derivatives(), dof_index_y, 1.);
217 
218  dof_id_type dof_index_z(-1);
219  ADReal lm_z;
220  if (_has_disp_z)
221  {
222  dof_index_z = dof->dof_number(_sys.number(), _lm_vars[2]->number(), 0);
223  lm_z = (*_sys.currentSolution())(dof_index_z);
224  Moose::derivInsert(lm_z.derivatives(), dof_index_z, 1.);
225  scaling_factor_z = _lm_vars[2]->scalingFactor();
226  }
227 
228  ADReal normal_pressure_value =
229  lm_x * _dof_to_normal_vector[dof](0) + lm_y * _dof_to_normal_vector[dof](1);
230  ADReal tangential_pressure_value =
231  lm_x * _dof_to_tangent_vectors[dof][0](0) + lm_y * _dof_to_tangent_vectors[dof][0](1);
232 
233  ADReal tangential_pressure_value_dir;
234 
235  if (_has_disp_z)
236  {
237  normal_pressure_value += lm_z * _dof_to_normal_vector[dof](2);
238  tangential_pressure_value += lm_z * _dof_to_tangent_vectors[dof][0](2);
239  tangential_pressure_value_dir = lm_x * _dof_to_tangent_vectors[dof][1](0) +
240  lm_y * _dof_to_tangent_vectors[dof][1](1) +
241  lm_z * _dof_to_tangent_vectors[dof][1](2);
242  }
243 
244  ADReal normal_dof_residual = std::min(normal_pressure_value, weighted_gap * c);
245  ADReal tangential_dof_residual;
246  ADReal tangential_dof_residual_dir;
247 
248  const ADReal & tangential_vel = *_tangential_vel_ptr[0];
249 
250  // Primal-dual active set strategy (PDASS)
251  if (!_has_disp_z)
252  {
253  if (normal_pressure_value.value() < _epsilon)
254  tangential_dof_residual = tangential_pressure_value;
255  else
256  {
257  const auto term_1 =
258  std::max(_mu * (normal_pressure_value + c * weighted_gap),
259  std::abs(tangential_pressure_value + c_t * tangential_vel * _dt)) *
260  tangential_pressure_value;
261  const auto term_2 = _mu * std::max(0.0, normal_pressure_value + c * weighted_gap) *
262  (tangential_pressure_value + c_t * tangential_vel * _dt);
263 
264  tangential_dof_residual = term_1 - term_2;
265  }
266  }
267  else
268  {
269  if (normal_pressure_value.value() < _epsilon)
270  {
271  tangential_dof_residual = tangential_pressure_value;
272  tangential_dof_residual_dir = tangential_pressure_value_dir;
273  }
274  else
275  {
276  const Real epsilon_sqrt = 1.0e-48;
277 
278  const auto lamdba_plus_cg = normal_pressure_value + c * weighted_gap;
279 
280  std::array<ADReal, 2> lambda_t_plus_ctu;
281  lambda_t_plus_ctu[0] = tangential_pressure_value + c_t * (*_tangential_vel_ptr[0]) * _dt;
282  lambda_t_plus_ctu[1] = tangential_pressure_value_dir + c_t * (*_tangential_vel_ptr[1]) * _dt;
283 
284  const auto term_1_x =
285  std::max(_mu * lamdba_plus_cg,
286  std::sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
287  lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
288  tangential_pressure_value;
289 
290  const auto term_1_y =
291  std::max(_mu * lamdba_plus_cg,
292  std::sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
293  lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
294  tangential_pressure_value_dir;
295 
296  const auto term_2_x = _mu * std::max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[0];
297 
298  const auto term_2_y = _mu * std::max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[1];
299 
300  tangential_dof_residual = term_1_x - term_2_x;
301  tangential_dof_residual_dir = term_1_y - term_2_y;
302  }
303  }
304 
305  // Compute the friction coefficient (constant or function)
306 
307  // Get index for normal constraint.
308  // We do this to get a decent Jacobian structure, which is key for the use of iterative solvers.
309  // Using old normal vector to avoid changes in the Jacobian structure within one time step
310 
311  Real ny, nz;
312  // Intially, use the current normal vector
313  if (_dof_to_old_normal_vector[dof].norm() < TOLERANCE)
314  {
315  ny = _dof_to_normal_vector[dof](1);
316  nz = _dof_to_normal_vector[dof](2);
317  }
318  else
319  {
320  ny = _dof_to_old_normal_vector[dof](1);
321  nz = _dof_to_old_normal_vector[dof](2);
322  }
323 
324  unsigned int component_normal = 0;
325 
326  // Consider constraint orientation to improve Jacobian structure
327  const Real threshold_for_Jacobian = _has_disp_z ? 1.0 / std::sqrt(3.0) : 1.0 / std::sqrt(2.0);
328 
329  if (std::abs(ny) > threshold_for_Jacobian)
330  component_normal = 1;
331  else if (std::abs(nz) > threshold_for_Jacobian)
332  component_normal = 2;
333 
335  _assembly,
336  std::array<ADReal, 1>{{normal_dof_residual}},
337  std::array<dof_id_type, 1>{{component_normal == 0
338  ? dof_index_x
339  : (component_normal == 1 ? dof_index_y : dof_index_z)}},
340  component_normal == 0 ? scaling_factor_x
341  : (component_normal == 1 ? scaling_factor_y : scaling_factor_z));
342 
344  _assembly,
345  std::array<ADReal, 1>{{tangential_dof_residual}},
346  std::array<dof_id_type, 1>{
347  {(component_normal == 0 || component_normal == 2) ? dof_index_y : dof_index_x}},
348  (component_normal == 0 || component_normal == 2) ? scaling_factor_y : scaling_factor_x);
349 
350  if (_has_disp_z)
352  _assembly,
353  std::array<ADReal, 1>{{tangential_dof_residual_dir}},
354  std::array<dof_id_type, 1>{
355  {(component_normal == 0 || component_normal == 1) ? dof_index_z : dof_index_x}},
356  (component_normal == 0 || component_normal == 1) ? scaling_factor_z : scaling_factor_x);
357 }
MooseMesh & _mesh
const ADVariableValue & _secondary_y_dot
y-velocity on the secondary face
virtual const NumericVector< Number > *const & currentSolution() const=0
virtual void computeQpProperties()
Computes properties that are functions only of the current quadrature point (_qp), e.g.
const MooseVariable *const _disp_z_var
The z displacement variable.
void addResidualsAndJacobian(Assembly &assembly, const Residuals &residuals, const Indices &dof_indices, Real scaling_factor)
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.
Computes the weighted gap that will later be used to enforce the zero-penetration mechanical contact ...
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
const bool _has_disp_z
For 2D mortar contact no displacement will be specified, so const pointers used.
std::unordered_map< const DofObject *, std::pair< ADReal, Real > > _dof_to_weighted_gap
A map from node to weighted gap and normalization (if requested)
void mooseError(Args &&... args)
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
const Real _c_t
Numerical factor used in the tangential constraints for convergence purposes.
const MooseVariable *const _disp_y_var
The y displacement variable.
T & set(const std::string &name, bool quiet_mode=false)
virtual void computeQpIProperties()
Computes properties that are functions both of _qp and _i, for example the weighted gap...
Elem const *const & _lower_primary_elem
const ADVariableValue & _primary_x_dot
x-velocity on the primary face
const std::vector< Real > & _JxW_msm
const Parallel::Communicator & _communicator
unsigned int _i
const bool _nodal
Whether the dof objects are nodal; if they&#39;re not, then they&#39;re elemental.
virtual void computeQpProperties() override
Computes properties that are functions only of the current quadrature point (_qp), e.g.
bool _normalize_c
Whether to normalize weighted gap by weighting function norm.
DualNumber< Real, DNDerivativeType, true > ADReal
const VariableTestValue & _test
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.
Elem const *const & _lower_secondary_elem
const AutomaticMortarGeneration & amg() const
std::unordered_map< const DofObject *, std::array< RealVectorValue, 2 > > _dof_to_tangent_vectors
A map from node to tangent vector (2D for now)
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)
const Real _epsilon
Small contact pressure value to trigger computation of frictional forces.
std::array< const ADReal *, 2 > _tangential_vel_ptr
An array of two pointers to avoid copies.
virtual void computeQpIProperties() override
Computes properties that are functions both of _qp and _i, for example the weighted gap...
std::array< MooseUtils::SemidynamicVector< Point, 9 >, 2 > getNodalTangents(const Elem &secondary_elem) const
std::map< unsigned int, unsigned int > getSecondaryIpToLowerElementMap(const Elem &lower_secondary_elem) const
std::unordered_map< const DofObject *, RealVectorValue > _dof_to_normal_vector
A map from node to normal vector (2D)
unsigned int number() const
auto norm(const T &a) -> decltype(std::abs(a))
Computes the weighted gap that will later be used to enforce the zero-penetration mechanical contact ...
const ADVariableValue & _primary_y_dot
y-velocity on the primary face
std::unordered_map< const DofObject *, std::array< ADReal, 2 > > _dof_to_weighted_tangential_velocity
A map from node to two weighted tangential velocities.
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.
std::vector< MooseVariable * > _lm_vars
Cartesian Lagrange multipliers for mechanical contact.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
registerMooseObject("ContactApp", ComputeFrictionalForceCartesianLMMechanicalContact)
const MooseVariable *const _disp_x_var
The x displacement variable.
ADRealVectorValue _qp_tangential_velocity_nodal
The value of the tangential velocity vectors at the current node.
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)
const ADVariableValue & _secondary_x_dot
x-velocity on the secondary face
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
const MooseArray< Real > & _coord
std::unordered_map< const DofObject *, RealVectorValue > _dof_to_old_normal_vector
A map from node to normal vector (2D) - old.
processor_id_type processor_id() const
unsigned int _qp
virtual void enforceConstraintOnDof(const DofObject *const dof) override
Method called from post().
uint8_t dof_id_type
const ADReal * _weighted_gap_ptr
A pointer members that can be used to help avoid copying ADReals.