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 "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 :
20 : registerMooseObject("ContactApp", ComputeFrictionalForceCartesianLMMechanicalContact);
21 :
22 : namespace
23 : {
24 : const InputParameters &
25 71 : assignVarsInParamsFriction(const InputParameters & params_in)
26 : {
27 : InputParameters & ret = const_cast<InputParameters &>(params_in);
28 71 : const auto & disp_x_name = ret.get<std::vector<VariableName>>("disp_x");
29 71 : if (disp_x_name.size() != 1)
30 0 : 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 142 : ret.set<VariableName>("secondary_variable") = disp_x_name[0];
34 71 : ret.set<VariableName>("primary_variable") = disp_x_name[0];
35 :
36 71 : return ret;
37 : }
38 : }
39 :
40 : InputParameters
41 142 : ComputeFrictionalForceCartesianLMMechanicalContact::validParams()
42 : {
43 142 : InputParameters params = ComputeWeightedGapCartesianLMMechanicalContact::validParams();
44 142 : params.addClassDescription("Computes mortar frictional forces.");
45 284 : params.addParam<Real>("c_t", 1e0, "Numerical parameter for tangential constraints");
46 284 : params.addParam<Real>(
47 : "epsilon",
48 284 : 1.0e-7,
49 : "Minimum value of contact pressure that will trigger frictional enforcement");
50 284 : params.addRangeCheckedParam<Real>(
51 : "mu", "mu > 0", "The friction coefficient for the Coulomb friction law");
52 142 : return params;
53 0 : }
54 :
55 71 : ComputeFrictionalForceCartesianLMMechanicalContact::
56 71 : ComputeFrictionalForceCartesianLMMechanicalContact(const InputParameters & parameters)
57 : : ComputeWeightedGapCartesianLMMechanicalContact(assignVarsInParamsFriction(parameters)),
58 142 : _c_t(getParam<Real>("c_t")),
59 71 : _secondary_x_dot(adCoupledDot("disp_x")),
60 71 : _primary_x_dot(adCoupledNeighborValueDot("disp_x")),
61 71 : _secondary_y_dot(adCoupledDot("disp_y")),
62 71 : _primary_y_dot(adCoupledNeighborValueDot("disp_y")),
63 71 : _secondary_z_dot(_has_disp_z ? &adCoupledDot("disp_z") : nullptr),
64 71 : _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot("disp_z") : nullptr),
65 142 : _mu(getParam<Real>("mu")),
66 213 : _epsilon(getParam<Real>("epsilon"))
67 : {
68 71 : }
69 :
70 : void
71 440324 : ComputeFrictionalForceCartesianLMMechanicalContact::computeQpProperties()
72 : {
73 440324 : ComputeWeightedGapCartesianLMMechanicalContact::computeQpProperties();
74 :
75 : // Trim derivatives
76 : const auto & primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
77 440324 : *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
78 : const auto & secondary_ip_lowerd_map =
79 440324 : amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
80 :
81 440324 : 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 440324 : {_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 440324 : {_secondary_x_dot[_qp], _secondary_y_dot[_qp], _has_disp_z ? (*_secondary_z_dot)[_qp] : 0}};
86 :
87 440324 : trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp_dot, false);
88 440324 : 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 440324 : 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 440324 : if (_has_disp_z)
100 : sec_z_dot = &secondary_disp_dot[2];
101 :
102 : // Build relative velocity vector
103 : ADRealVectorValue relative_velocity;
104 :
105 440324 : if (_has_disp_z)
106 110592 : relative_velocity = {sec_x_dot - prim_x_dot, sec_y_dot - prim_y_dot, *sec_z_dot - *prim_z_dot};
107 : else
108 659464 : relative_velocity = {sec_x_dot - prim_x_dot, sec_y_dot - prim_y_dot, 0.0};
109 :
110 880648 : _qp_tangential_velocity_nodal = relative_velocity * (_JxW_msm[_qp] * _coord[_qp]);
111 440324 : }
112 :
113 : void
114 1101832 : ComputeFrictionalForceCartesianLMMechanicalContact::computeQpIProperties()
115 : {
116 1101832 : ComputeWeightedGapCartesianLMMechanicalContact::computeQpIProperties();
117 :
118 1101832 : const auto & nodal_tangents = amg().getNodalTangents(*_lower_secondary_elem);
119 : // Get the _dof_to_weighted_tangential_velocity map
120 : const DofObject * const dof =
121 1101832 : _lm_vars[0]->isNodal() ? static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i))
122 0 : : static_cast<const DofObject *>(_lower_secondary_elem);
123 :
124 : _dof_to_weighted_tangential_velocity[dof][0] +=
125 2203664 : _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 1101832 : if (_has_disp_z)
129 : _dof_to_weighted_tangential_velocity[dof][1] +=
130 884736 : _test[_i][_qp] * _qp_tangential_velocity_nodal * nodal_tangents[1][_i];
131 1101832 : }
132 :
133 : void
134 5334 : ComputeFrictionalForceCartesianLMMechanicalContact::residualSetup()
135 : {
136 5334 : ComputeWeightedGapCartesianLMMechanicalContact::residualSetup();
137 : _dof_to_weighted_tangential_velocity.clear();
138 5334 : }
139 :
140 : void
141 3950 : ComputeFrictionalForceCartesianLMMechanicalContact::post()
142 : {
143 3950 : Moose::Mortar::Contact::communicateGaps(
144 3950 : _dof_to_weighted_gap, _mesh, _nodal, _normalize_c, _communicator, false);
145 3950 : Moose::Mortar::Contact::communicateVelocities(
146 3950 : _dof_to_weighted_tangential_velocity, _mesh, _nodal, _communicator, false);
147 :
148 : // Enforce frictional complementarity constraints
149 64312 : for (const auto & pr : _dof_to_weighted_tangential_velocity)
150 : {
151 60362 : const DofObject * const dof = pr.first;
152 :
153 60362 : if (dof->processor_id() != this->processor_id())
154 1064 : continue;
155 :
156 : auto & weighted_gap_pr = _dof_to_weighted_gap[dof];
157 59298 : _weighted_gap_ptr = &weighted_gap_pr.first;
158 59298 : _normalization_ptr = &weighted_gap_pr.second;
159 59298 : _tangential_vel_ptr[0] = &(pr.second[0]);
160 :
161 59298 : if (_has_disp_z)
162 1728 : _tangential_vel_ptr[1] = &(pr.second[1]);
163 :
164 59298 : enforceConstraintOnDof(dof);
165 : }
166 3950 : }
167 :
168 : void
169 1384 : ComputeFrictionalForceCartesianLMMechanicalContact::incorrectEdgeDroppingPost(
170 : const std::unordered_set<const Node *> & inactive_lm_nodes)
171 : {
172 1384 : Moose::Mortar::Contact::communicateGaps(
173 1384 : _dof_to_weighted_gap, _mesh, _nodal, _normalize_c, _communicator, false);
174 1384 : Moose::Mortar::Contact::communicateVelocities(
175 1384 : _dof_to_weighted_tangential_velocity, _mesh, _nodal, _communicator, false);
176 :
177 : // Enforce frictional complementarity constraints
178 93490 : for (const auto & pr : _dof_to_weighted_tangential_velocity)
179 : {
180 92106 : const DofObject * const dof = pr.first;
181 :
182 : // If node inactive, skip
183 92106 : if ((inactive_lm_nodes.find(static_cast<const Node *>(dof)) != inactive_lm_nodes.end()) ||
184 : (dof->processor_id() != this->processor_id()))
185 2292 : continue;
186 :
187 89814 : _weighted_gap_ptr = &_dof_to_weighted_gap[dof].first;
188 89814 : _normalization_ptr = &_dof_to_weighted_gap[dof].second;
189 89814 : _tangential_vel_ptr[0] = &pr.second[0];
190 :
191 89814 : if (_has_disp_z)
192 0 : _tangential_vel_ptr[1] = &pr.second[1];
193 :
194 89814 : enforceConstraintOnDof(dof);
195 : }
196 1384 : }
197 :
198 : void
199 149112 : ComputeFrictionalForceCartesianLMMechanicalContact::enforceConstraintOnDof(
200 : const DofObject * const dof)
201 : {
202 149112 : const auto & weighted_gap = *_weighted_gap_ptr;
203 149112 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
204 149112 : const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
205 :
206 149112 : const auto dof_index_x = dof->dof_number(_sys.number(), _lm_vars[0]->number(), 0);
207 149112 : const auto dof_index_y = dof->dof_number(_sys.number(), _lm_vars[1]->number(), 0);
208 149112 : const Real scaling_factor_x = _lm_vars[0]->scalingFactor();
209 149112 : const Real scaling_factor_y = _lm_vars[1]->scalingFactor();
210 : Real scaling_factor_z = 1;
211 :
212 149112 : ADReal lm_x = (*_sys.currentSolution())(dof_index_x);
213 149112 : 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 149112 : if (_has_disp_z)
221 : {
222 1728 : dof_index_z = dof->dof_number(_sys.number(), _lm_vars[2]->number(), 0);
223 1728 : lm_z = (*_sys.currentSolution())(dof_index_z);
224 : Moose::derivInsert(lm_z.derivatives(), dof_index_z, 1.);
225 1728 : scaling_factor_z = _lm_vars[2]->scalingFactor();
226 : }
227 :
228 : ADReal normal_pressure_value =
229 149112 : lm_x * _dof_to_normal_vector[dof](0) + lm_y * _dof_to_normal_vector[dof](1);
230 : ADReal tangential_pressure_value =
231 149112 : 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 149112 : if (_has_disp_z)
236 : {
237 1728 : normal_pressure_value += lm_z * _dof_to_normal_vector[dof](2);
238 1728 : tangential_pressure_value += lm_z * _dof_to_tangent_vectors[dof][0](2);
239 1728 : tangential_pressure_value_dir = lm_x * _dof_to_tangent_vectors[dof][1](0) +
240 1728 : lm_y * _dof_to_tangent_vectors[dof][1](1) +
241 1728 : lm_z * _dof_to_tangent_vectors[dof][1](2);
242 : }
243 :
244 149112 : 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 149112 : const ADReal & tangential_vel = *_tangential_vel_ptr[0];
249 :
250 : // Primal-dual active set strategy (PDASS)
251 149112 : if (!_has_disp_z)
252 : {
253 147384 : if (normal_pressure_value.value() < _epsilon)
254 99147 : tangential_dof_residual = tangential_pressure_value;
255 : else
256 : {
257 : const auto term_1 =
258 96474 : std::max(_mu * (normal_pressure_value + c * weighted_gap),
259 96474 : std::abs(tangential_pressure_value + c_t * tangential_vel * _dt)) *
260 : tangential_pressure_value;
261 48237 : const auto term_2 = _mu * std::max(0.0, normal_pressure_value + c * weighted_gap) *
262 96474 : (tangential_pressure_value + c_t * tangential_vel * _dt);
263 :
264 48237 : tangential_dof_residual = term_1 - term_2;
265 : }
266 : }
267 : else
268 : {
269 1728 : if (normal_pressure_value.value() < _epsilon)
270 : {
271 300 : tangential_dof_residual = tangential_pressure_value;
272 300 : tangential_dof_residual_dir = tangential_pressure_value_dir;
273 : }
274 : else
275 : {
276 : const Real epsilon_sqrt = 1.0e-48;
277 :
278 1428 : const auto lamdba_plus_cg = normal_pressure_value + c * weighted_gap;
279 :
280 : std::array<ADReal, 2> lambda_t_plus_ctu;
281 2856 : lambda_t_plus_ctu[0] = tangential_pressure_value + c_t * (*_tangential_vel_ptr[0]) * _dt;
282 2856 : lambda_t_plus_ctu[1] = tangential_pressure_value_dir + c_t * (*_tangential_vel_ptr[1]) * _dt;
283 :
284 : const auto term_1_x =
285 2856 : std::max(_mu * lamdba_plus_cg,
286 1428 : std::sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
287 1428 : lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
288 : tangential_pressure_value;
289 :
290 : const auto term_1_y =
291 1428 : std::max(_mu * lamdba_plus_cg,
292 1428 : std::sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
293 1428 : lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
294 : tangential_pressure_value_dir;
295 :
296 1428 : const auto term_2_x = _mu * std::max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[0];
297 :
298 2856 : const auto term_2_y = _mu * std::max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[1];
299 :
300 1428 : tangential_dof_residual = term_1_x - term_2_x;
301 1428 : 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 149112 : if (_dof_to_old_normal_vector[dof].norm() < TOLERANCE)
314 : {
315 20745 : ny = _dof_to_normal_vector[dof](1);
316 20745 : nz = _dof_to_normal_vector[dof](2);
317 : }
318 : else
319 : {
320 128367 : ny = _dof_to_old_normal_vector[dof](1);
321 128367 : 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 149112 : const Real threshold_for_Jacobian = _has_disp_z ? 1.0 / std::sqrt(3.0) : 1.0 / std::sqrt(2.0);
328 :
329 149112 : if (std::abs(ny) > threshold_for_Jacobian)
330 : component_normal = 1;
331 61036 : else if (std::abs(nz) > threshold_for_Jacobian)
332 : component_normal = 2;
333 :
334 298224 : addResidualsAndJacobian(
335 : _assembly,
336 149112 : std::array<ADReal, 1>{{normal_dof_residual}},
337 298224 : std::array<dof_id_type, 1>{{component_normal == 0
338 89804 : ? 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 :
343 149112 : addResidualsAndJacobian(
344 : _assembly,
345 149112 : std::array<ADReal, 1>{{tangential_dof_residual}},
346 298224 : std::array<dof_id_type, 1>{
347 149112 : {(component_normal == 0 || component_normal == 2) ? dof_index_y : dof_index_x}},
348 149112 : (component_normal == 0 || component_normal == 2) ? scaling_factor_y : scaling_factor_x);
349 :
350 149112 : if (_has_disp_z)
351 3456 : addResidualsAndJacobian(
352 : _assembly,
353 1728 : std::array<ADReal, 1>{{tangential_dof_residual_dir}},
354 1728 : 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 449064 : }
|
