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