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