Line data Source code
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 :
10 : #include "ComputeFrictionalForceCartesianLMMechanicalContact.h"
11 : #include "MortarContactUtils.h"
12 : #include "DisplacedProblem.h"
13 : #include "Assembly.h"
14 : #include "AutomaticMortarGeneration.h"
15 : #include "metaphysicl/metaphysicl_version.h"
16 : #include "metaphysicl/dualsemidynamicsparsenumberarray.h"
17 : #include "metaphysicl/parallel_dualnumber.h"
18 : #if METAPHYSICL_MAJOR_VERSION < 2
19 : #include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
20 : #else
21 : #include "metaphysicl/parallel_dynamic_array_wrapper.h"
22 : #endif
23 : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
24 : #include "timpi/parallel_sync.h"
25 :
26 : registerMooseObject("ContactApp", ComputeFrictionalForceCartesianLMMechanicalContact);
27 :
28 : namespace
29 : {
30 : const InputParameters &
31 78 : assignVarsInParamsFriction(const InputParameters & params_in)
32 : {
33 : InputParameters & ret = const_cast<InputParameters &>(params_in);
34 78 : const auto & disp_x_name = ret.get<std::vector<VariableName>>("disp_x");
35 78 : if (disp_x_name.size() != 1)
36 0 : mooseError("We require that the disp_x parameter have exactly one coupled name");
37 :
38 : // We do this so we don't get any variable errors during MortarConstraint(Base) construction
39 156 : ret.set<VariableName>("secondary_variable") = disp_x_name[0];
40 78 : ret.set<VariableName>("primary_variable") = disp_x_name[0];
41 :
42 78 : return ret;
43 : }
44 : }
45 :
46 : InputParameters
47 156 : ComputeFrictionalForceCartesianLMMechanicalContact::validParams()
48 : {
49 156 : InputParameters params = ComputeWeightedGapCartesianLMMechanicalContact::validParams();
50 156 : params.addClassDescription("Computes mortar frictional forces.");
51 312 : params.addParam<Real>("c_t", 1e0, "Numerical parameter for tangential constraints");
52 312 : params.addParam<Real>(
53 : "epsilon",
54 312 : 1.0e-7,
55 : "Minimum value of contact pressure that will trigger frictional enforcement");
56 312 : params.addRangeCheckedParam<Real>(
57 : "mu", "mu > 0", "The friction coefficient for the Coulomb friction law");
58 156 : return params;
59 0 : }
60 :
61 78 : ComputeFrictionalForceCartesianLMMechanicalContact::
62 78 : ComputeFrictionalForceCartesianLMMechanicalContact(const InputParameters & parameters)
63 : : ComputeWeightedGapCartesianLMMechanicalContact(assignVarsInParamsFriction(parameters)),
64 156 : _c_t(getParam<Real>("c_t")),
65 78 : _secondary_x_dot(adCoupledDot("disp_x")),
66 78 : _primary_x_dot(adCoupledNeighborValueDot("disp_x")),
67 78 : _secondary_y_dot(adCoupledDot("disp_y")),
68 78 : _primary_y_dot(adCoupledNeighborValueDot("disp_y")),
69 78 : _secondary_z_dot(_has_disp_z ? &adCoupledDot("disp_z") : nullptr),
70 78 : _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot("disp_z") : nullptr),
71 156 : _mu(getParam<Real>("mu")),
72 234 : _epsilon(getParam<Real>("epsilon"))
73 : {
74 78 : }
75 :
76 : void
77 564604 : ComputeFrictionalForceCartesianLMMechanicalContact::computeQpProperties()
78 : {
79 564604 : ComputeWeightedGapCartesianLMMechanicalContact::computeQpProperties();
80 :
81 : // Trim derivatives
82 : const auto & primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
83 564604 : *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
84 : const auto & secondary_ip_lowerd_map =
85 564604 : amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
86 :
87 564604 : std::array<const MooseVariable *, 3> var_array{{_disp_x_var, _disp_y_var, _disp_z_var}};
88 : std::array<ADReal, 3> primary_disp_dot{
89 564604 : {_primary_x_dot[_qp], _primary_y_dot[_qp], _has_disp_z ? (*_primary_z_dot)[_qp] : 0}};
90 : std::array<ADReal, 3> secondary_disp_dot{
91 564604 : {_secondary_x_dot[_qp], _secondary_y_dot[_qp], _has_disp_z ? (*_secondary_z_dot)[_qp] : 0}};
92 :
93 564604 : trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp_dot, false);
94 564604 : trimInteriorNodeDerivatives(secondary_ip_lowerd_map, var_array, secondary_disp_dot, true);
95 :
96 : const ADReal & prim_x_dot = primary_disp_dot[0];
97 : const ADReal & prim_y_dot = primary_disp_dot[1];
98 : const ADReal * prim_z_dot = nullptr;
99 564604 : if (_has_disp_z)
100 : prim_z_dot = &primary_disp_dot[2];
101 :
102 : const ADReal & sec_x_dot = secondary_disp_dot[0];
103 : const ADReal & sec_y_dot = secondary_disp_dot[1];
104 : const ADReal * sec_z_dot = nullptr;
105 564604 : if (_has_disp_z)
106 : sec_z_dot = &secondary_disp_dot[2];
107 :
108 : // Build relative velocity vector
109 : ADRealVectorValue relative_velocity;
110 :
111 564604 : if (_has_disp_z)
112 169984 : relative_velocity = {sec_x_dot - prim_x_dot, sec_y_dot - prim_y_dot, *sec_z_dot - *prim_z_dot};
113 : else
114 789240 : relative_velocity = {sec_x_dot - prim_x_dot, sec_y_dot - prim_y_dot, 0.0};
115 :
116 1129208 : _qp_tangential_velocity_nodal = relative_velocity * (_JxW_msm[_qp] * _coord[_qp]);
117 564604 : }
118 :
119 : void
120 1469176 : ComputeFrictionalForceCartesianLMMechanicalContact::computeQpIProperties()
121 : {
122 1469176 : ComputeWeightedGapCartesianLMMechanicalContact::computeQpIProperties();
123 :
124 1469176 : const auto & nodal_tangents = amg().getNodalTangents(*_lower_secondary_elem);
125 : // Get the _dof_to_weighted_tangential_velocity map
126 : const DofObject * const dof =
127 1469176 : _lm_vars[0]->isNodal() ? static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i))
128 0 : : static_cast<const DofObject *>(_lower_secondary_elem);
129 :
130 : _dof_to_weighted_tangential_velocity[dof][0] +=
131 2938352 : _test[_i][_qp] * _qp_tangential_velocity_nodal * nodal_tangents[0][_i];
132 :
133 : // Get the _dof_to_weighted_tangential_velocity map for a second direction
134 1469176 : if (_has_disp_z)
135 : _dof_to_weighted_tangential_velocity[dof][1] +=
136 1359872 : _test[_i][_qp] * _qp_tangential_velocity_nodal * nodal_tangents[1][_i];
137 1469176 : }
138 :
139 : void
140 6069 : ComputeFrictionalForceCartesianLMMechanicalContact::residualSetup()
141 : {
142 6069 : ComputeWeightedGapCartesianLMMechanicalContact::residualSetup();
143 : _dof_to_weighted_tangential_velocity.clear();
144 6069 : }
145 :
146 : void
147 4483 : ComputeFrictionalForceCartesianLMMechanicalContact::post()
148 : {
149 4483 : Moose::Mortar::Contact::communicateGaps(
150 4483 : _dof_to_weighted_gap, _mesh, _nodal, _normalize_c, _communicator, false);
151 4483 : Moose::Mortar::Contact::communicateVelocities(
152 4483 : _dof_to_weighted_tangential_velocity, _mesh, _nodal, _communicator, false);
153 :
154 : // Enforce frictional complementarity constraints
155 75326 : for (const auto & pr : _dof_to_weighted_tangential_velocity)
156 : {
157 70843 : const DofObject * const dof = pr.first;
158 :
159 70843 : if (dof->processor_id() != this->processor_id())
160 1064 : continue;
161 :
162 : auto & weighted_gap_pr = _dof_to_weighted_gap[dof];
163 69779 : _weighted_gap_ptr = &weighted_gap_pr.first;
164 69779 : _normalization_ptr = &weighted_gap_pr.second;
165 69779 : _tangential_vel_ptr[0] = &(pr.second[0]);
166 :
167 69779 : if (_has_disp_z)
168 2656 : _tangential_vel_ptr[1] = &(pr.second[1]);
169 :
170 69779 : enforceConstraintOnDof(dof);
171 : }
172 4483 : }
173 :
174 : void
175 1586 : ComputeFrictionalForceCartesianLMMechanicalContact::incorrectEdgeDroppingPost(
176 : const std::unordered_set<const Node *> & inactive_lm_nodes)
177 : {
178 1586 : Moose::Mortar::Contact::communicateGaps(
179 1586 : _dof_to_weighted_gap, _mesh, _nodal, _normalize_c, _communicator, false);
180 1586 : Moose::Mortar::Contact::communicateVelocities(
181 1586 : _dof_to_weighted_tangential_velocity, _mesh, _nodal, _communicator, false);
182 :
183 : // Enforce frictional complementarity constraints
184 114074 : for (const auto & pr : _dof_to_weighted_tangential_velocity)
185 : {
186 112488 : const DofObject * const dof = pr.first;
187 :
188 : // If node inactive, skip
189 112488 : if ((inactive_lm_nodes.find(static_cast<const Node *>(dof)) != inactive_lm_nodes.end()) ||
190 110268 : (dof->processor_id() != this->processor_id()))
191 2696 : continue;
192 :
193 109792 : _weighted_gap_ptr = &_dof_to_weighted_gap[dof].first;
194 109792 : _normalization_ptr = &_dof_to_weighted_gap[dof].second;
195 109792 : _tangential_vel_ptr[0] = &pr.second[0];
196 :
197 109792 : if (_has_disp_z)
198 0 : _tangential_vel_ptr[1] = &pr.second[1];
199 :
200 109792 : enforceConstraintOnDof(dof);
201 : }
202 1586 : }
203 :
204 : void
205 179571 : ComputeFrictionalForceCartesianLMMechanicalContact::enforceConstraintOnDof(
206 : const DofObject * const dof)
207 : {
208 : using std::abs, std::sqrt, std::max, std::min;
209 :
210 179571 : const auto & weighted_gap = *_weighted_gap_ptr;
211 179571 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
212 179571 : const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
213 :
214 179571 : const auto dof_index_x = dof->dof_number(_sys.number(), _lm_vars[0]->number(), 0);
215 179571 : const auto dof_index_y = dof->dof_number(_sys.number(), _lm_vars[1]->number(), 0);
216 179571 : const Real scaling_factor_x = _lm_vars[0]->scalingFactor();
217 179571 : const Real scaling_factor_y = _lm_vars[1]->scalingFactor();
218 : Real scaling_factor_z = 1;
219 :
220 179571 : ADReal lm_x = (*_sys.currentSolution())(dof_index_x);
221 179571 : ADReal lm_y = (*_sys.currentSolution())(dof_index_y);
222 :
223 : Moose::derivInsert(lm_x.derivatives(), dof_index_x, 1.);
224 : Moose::derivInsert(lm_y.derivatives(), dof_index_y, 1.);
225 :
226 : dof_id_type dof_index_z(-1);
227 : ADReal lm_z;
228 179571 : if (_has_disp_z)
229 : {
230 2656 : dof_index_z = dof->dof_number(_sys.number(), _lm_vars[2]->number(), 0);
231 2656 : lm_z = (*_sys.currentSolution())(dof_index_z);
232 : Moose::derivInsert(lm_z.derivatives(), dof_index_z, 1.);
233 2656 : scaling_factor_z = _lm_vars[2]->scalingFactor();
234 : }
235 :
236 : ADReal normal_pressure_value =
237 179571 : lm_x * _dof_to_normal_vector[dof](0) + lm_y * _dof_to_normal_vector[dof](1);
238 : ADReal tangential_pressure_value =
239 179571 : lm_x * _dof_to_tangent_vectors[dof][0](0) + lm_y * _dof_to_tangent_vectors[dof][0](1);
240 :
241 : ADReal tangential_pressure_value_dir;
242 :
243 179571 : if (_has_disp_z)
244 : {
245 2656 : normal_pressure_value += lm_z * _dof_to_normal_vector[dof](2);
246 2656 : tangential_pressure_value += lm_z * _dof_to_tangent_vectors[dof][0](2);
247 2656 : tangential_pressure_value_dir = lm_x * _dof_to_tangent_vectors[dof][1](0) +
248 2656 : lm_y * _dof_to_tangent_vectors[dof][1](1) +
249 2656 : lm_z * _dof_to_tangent_vectors[dof][1](2);
250 : }
251 :
252 179571 : ADReal normal_dof_residual = min(normal_pressure_value, weighted_gap * c);
253 : ADReal tangential_dof_residual;
254 : ADReal tangential_dof_residual_dir;
255 :
256 179571 : const ADReal & tangential_vel = *_tangential_vel_ptr[0];
257 :
258 : // Primal-dual active set strategy (PDASS)
259 179571 : if (!_has_disp_z)
260 : {
261 176915 : if (normal_pressure_value.value() < _epsilon)
262 120699 : tangential_dof_residual = tangential_pressure_value;
263 : else
264 : {
265 112432 : const auto term_1 = max(_mu * (normal_pressure_value + c * weighted_gap),
266 112432 : abs(tangential_pressure_value + c_t * tangential_vel * _dt)) *
267 : tangential_pressure_value;
268 56216 : const auto term_2 = _mu * max(0.0, normal_pressure_value + c * weighted_gap) *
269 112432 : (tangential_pressure_value + c_t * tangential_vel * _dt);
270 :
271 56216 : tangential_dof_residual = term_1 - term_2;
272 : }
273 : }
274 : else
275 : {
276 2656 : if (normal_pressure_value.value() < _epsilon)
277 : {
278 500 : tangential_dof_residual = tangential_pressure_value;
279 500 : tangential_dof_residual_dir = tangential_pressure_value_dir;
280 : }
281 : else
282 : {
283 : const Real epsilon_sqrt = 1.0e-48;
284 :
285 2156 : const auto lamdba_plus_cg = normal_pressure_value + c * weighted_gap;
286 :
287 : std::array<ADReal, 2> lambda_t_plus_ctu;
288 4312 : lambda_t_plus_ctu[0] = tangential_pressure_value + c_t * (*_tangential_vel_ptr[0]) * _dt;
289 4312 : lambda_t_plus_ctu[1] = tangential_pressure_value_dir + c_t * (*_tangential_vel_ptr[1]) * _dt;
290 :
291 4312 : const auto term_1_x = max(_mu * lamdba_plus_cg,
292 2156 : sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
293 2156 : lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
294 : tangential_pressure_value;
295 :
296 2156 : const auto term_1_y = max(_mu * lamdba_plus_cg,
297 2156 : sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
298 2156 : lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
299 : tangential_pressure_value_dir;
300 :
301 2156 : const auto term_2_x = _mu * max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[0];
302 :
303 4312 : const auto term_2_y = _mu * max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[1];
304 :
305 2156 : tangential_dof_residual = term_1_x - term_2_x;
306 2156 : tangential_dof_residual_dir = term_1_y - term_2_y;
307 : }
308 : }
309 :
310 : // Compute the friction coefficient (constant or function)
311 :
312 : // Get index for normal constraint.
313 : // We do this to get a decent Jacobian structure, which is key for the use of iterative solvers.
314 : // Using old normal vector to avoid changes in the Jacobian structure within one time step
315 :
316 : Real ny, nz;
317 : // Intially, use the current normal vector
318 179571 : if (_dof_to_old_normal_vector[dof].norm() < TOLERANCE)
319 : {
320 21790 : ny = _dof_to_normal_vector[dof](1);
321 21790 : nz = _dof_to_normal_vector[dof](2);
322 : }
323 : else
324 : {
325 157781 : ny = _dof_to_old_normal_vector[dof](1);
326 157781 : nz = _dof_to_old_normal_vector[dof](2);
327 : }
328 :
329 : unsigned int component_normal = 0;
330 :
331 : // Consider constraint orientation to improve Jacobian structure
332 179571 : const Real threshold_for_Jacobian = _has_disp_z ? 1.0 / sqrt(3.0) : 1.0 / sqrt(2.0);
333 :
334 179571 : if (abs(ny) > threshold_for_Jacobian)
335 : component_normal = 1;
336 71901 : else if (abs(nz) > threshold_for_Jacobian)
337 : component_normal = 2;
338 :
339 359142 : addResidualsAndJacobian(
340 : _assembly,
341 359142 : std::array<ADReal, 1>{{normal_dof_residual}},
342 179571 : std::array<dof_id_type, 1>{{component_normal == 0
343 110326 : ? dof_index_x
344 : : (component_normal == 1 ? dof_index_y : dof_index_z)}},
345 : component_normal == 0 ? scaling_factor_x
346 : : (component_normal == 1 ? scaling_factor_y : scaling_factor_z));
347 :
348 179571 : addResidualsAndJacobian(
349 : _assembly,
350 359142 : std::array<ADReal, 1>{{tangential_dof_residual}},
351 179571 : std::array<dof_id_type, 1>{
352 179571 : {(component_normal == 0 || component_normal == 2) ? dof_index_y : dof_index_x}},
353 179571 : (component_normal == 0 || component_normal == 2) ? scaling_factor_y : scaling_factor_x);
354 :
355 179571 : if (_has_disp_z)
356 5312 : addResidualsAndJacobian(
357 : _assembly,
358 5312 : std::array<ADReal, 1>{{tangential_dof_residual_dir}},
359 5312 : std::array<dof_id_type, 1>{
360 : {(component_normal == 0 || component_normal == 1) ? dof_index_z : dof_index_x}},
361 : (component_normal == 0 || component_normal == 1) ? scaling_factor_z : scaling_factor_x);
362 179571 : }
|
