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 "ComputeWeightedGapCartesianLMMechanicalContact.h"
11 : #include "DisplacedProblem.h"
12 : #include "MortarContactUtils.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", ComputeWeightedGapCartesianLMMechanicalContact);
21 :
22 : namespace
23 : {
24 : const InputParameters &
25 133 : assignVarsInParamsCartesianWeightedGap(const InputParameters & params_in)
26 : {
27 : InputParameters & ret = const_cast<InputParameters &>(params_in);
28 133 : const auto & disp_x_name = ret.get<std::vector<VariableName>>("disp_x");
29 133 : 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 266 : ret.set<VariableName>("secondary_variable") = disp_x_name[0];
34 133 : ret.set<VariableName>("primary_variable") = disp_x_name[0];
35 :
36 133 : return ret;
37 : }
38 : }
39 :
40 : InputParameters
41 266 : ComputeWeightedGapCartesianLMMechanicalContact::validParams()
42 : {
43 266 : InputParameters params = ADMortarConstraint::validParams();
44 266 : params.addClassDescription("Computes the weighted gap that will later be used to enforce the "
45 : "zero-penetration mechanical contact conditions");
46 266 : params.suppressParameter<VariableName>("secondary_variable");
47 266 : params.suppressParameter<VariableName>("primary_variable");
48 532 : params.addRequiredCoupledVar("disp_x", "The x displacement variable");
49 532 : params.addRequiredCoupledVar("disp_y", "The y displacement variable");
50 532 : params.addCoupledVar("disp_z", "The z displacement variable");
51 532 : params.addParam<Real>(
52 532 : "c", 1e6, "Parameter for balancing the size of the gap and contact pressure");
53 532 : params.addRequiredCoupledVar("lm_x",
54 : "Mechanical contact Lagrange multiplier along the x Cartesian axis");
55 532 : params.addRequiredCoupledVar(
56 : "lm_y", "Mechanical contact Lagrange multiplier along the y Cartesian axis.");
57 532 : params.addCoupledVar("lm_z",
58 : "Mechanical contact Lagrange multiplier along the z Cartesian axis.");
59 266 : params.set<bool>("interpolate_normals") = false;
60 532 : params.addParam<bool>(
61 : "normalize_c",
62 532 : false,
63 : "Whether to normalize c by weighting function norm. When unnormalized "
64 : "the value of c effectively depends on element size since in the constraint we compare nodal "
65 : "Lagrange Multiplier values to integrated gap values (LM nodal value is independent of "
66 : "element size, where integrated values are dependent on element size).");
67 266 : params.set<bool>("use_displaced_mesh") = true;
68 266 : return params;
69 0 : }
70 :
71 133 : ComputeWeightedGapCartesianLMMechanicalContact::ComputeWeightedGapCartesianLMMechanicalContact(
72 133 : const InputParameters & parameters)
73 : : ADMortarConstraint(assignVarsInParamsCartesianWeightedGap(parameters)),
74 133 : _secondary_disp_x(adCoupledValue("disp_x")),
75 133 : _primary_disp_x(adCoupledNeighborValue("disp_x")),
76 133 : _secondary_disp_y(adCoupledValue("disp_y")),
77 133 : _primary_disp_y(adCoupledNeighborValue("disp_y")),
78 133 : _has_disp_z(isCoupled("disp_z")),
79 133 : _secondary_disp_z(_has_disp_z ? &adCoupledValue("disp_z") : nullptr),
80 133 : _primary_disp_z(_has_disp_z ? &adCoupledNeighborValue("disp_z") : nullptr),
81 266 : _c(getParam<Real>("c")),
82 266 : _normalize_c(getParam<bool>("normalize_c")),
83 133 : _nodal(getVar("disp_x", 0)->feType().family == LAGRANGE),
84 133 : _disp_x_var(getVar("disp_x", 0)),
85 133 : _disp_y_var(getVar("disp_y", 0)),
86 415 : _disp_z_var(_has_disp_z ? getVar("disp_z", 0) : nullptr)
87 : {
88 133 : _lm_vars.push_back(getVar("lm_x", 0));
89 133 : _lm_vars.push_back(getVar("lm_y", 0));
90 :
91 266 : if (isParamValid("lm_z") ^ _has_disp_z)
92 0 : paramError("lm_z",
93 : "In three-dimensions, both the Z Lagrange multiplier and the Z displacement need to "
94 : "be provided");
95 :
96 133 : if (_has_disp_z)
97 32 : _lm_vars.push_back(getVar("lm_z", 0));
98 :
99 266 : if (!getParam<bool>("use_displaced_mesh"))
100 0 : paramError(
101 : "use_displaced_mesh",
102 : "'use_displaced_mesh' must be true for the ComputeWeightedGapLMMechanicalContact object");
103 :
104 415 : for (const auto i : index_range(_lm_vars))
105 282 : if (!_lm_vars[i]->isNodal())
106 0 : if (_lm_vars[i]->feType().order != static_cast<Order>(0))
107 0 : mooseError("Normal contact constraints only support elemental variables of CONSTANT order");
108 133 : }
109 :
110 0 : ADReal ComputeWeightedGapCartesianLMMechanicalContact::computeQpResidual(Moose::MortarType)
111 : {
112 0 : mooseError("We should never call computeQpResidual for ComputeWeightedGapLMMechanicalContact");
113 : }
114 :
115 : void
116 798836 : ComputeWeightedGapCartesianLMMechanicalContact::computeQpProperties()
117 : {
118 : // Trim interior node variable derivatives
119 : const auto & primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
120 798836 : *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
121 : const auto & secondary_ip_lowerd_map =
122 798836 : amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
123 :
124 798836 : std::array<const MooseVariable *, 3> var_array{{_disp_x_var, _disp_y_var, _disp_z_var}};
125 : std::array<ADReal, 3> primary_disp{
126 798836 : {_primary_disp_x[_qp], _primary_disp_y[_qp], _has_disp_z ? (*_primary_disp_z)[_qp] : 0}};
127 798836 : std::array<ADReal, 3> secondary_disp{{_secondary_disp_x[_qp],
128 798836 : _secondary_disp_y[_qp],
129 798836 : _has_disp_z ? (*_secondary_disp_z)[_qp] : 0}};
130 :
131 798836 : trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp, false);
132 798836 : trimInteriorNodeDerivatives(secondary_ip_lowerd_map, var_array, secondary_disp, true);
133 :
134 : const ADReal & prim_x = primary_disp[0];
135 : const ADReal & prim_y = primary_disp[1];
136 : const ADReal * prim_z = nullptr;
137 798836 : if (_has_disp_z)
138 : prim_z = &primary_disp[2];
139 :
140 : const ADReal & sec_x = secondary_disp[0];
141 : const ADReal & sec_y = secondary_disp[1];
142 : const ADReal * sec_z = nullptr;
143 798836 : if (_has_disp_z)
144 : sec_z = &secondary_disp[2];
145 :
146 : // Compute gap vector
147 798836 : ADRealVectorValue gap_vec = _phys_points_primary[_qp] - _phys_points_secondary[_qp];
148 :
149 798836 : gap_vec(0).derivatives() = prim_x.derivatives() - sec_x.derivatives();
150 798836 : gap_vec(1).derivatives() = prim_y.derivatives() - sec_y.derivatives();
151 798836 : if (_has_disp_z)
152 282624 : gap_vec(2).derivatives() = prim_z->derivatives() - sec_z->derivatives();
153 :
154 : // Compute integration point quantities
155 798836 : _qp_gap_nodal = gap_vec * (_JxW_msm[_qp] * _coord[_qp]);
156 :
157 : // To do normalization of constraint coefficient (c_n)
158 798836 : _qp_factor = _JxW_msm[_qp] * _coord[_qp];
159 798836 : }
160 :
161 : void
162 2162920 : ComputeWeightedGapCartesianLMMechanicalContact::computeQpIProperties()
163 : {
164 : mooseAssert(_normals.size() == _lower_secondary_elem->n_nodes(),
165 : "Making sure that _normals is the expected size");
166 :
167 : // Get the _dof_to_weighted_gap map
168 2162920 : const DofObject * dof = _lm_vars[0]->isNodal()
169 2162920 : ? static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i))
170 0 : : static_cast<const DofObject *>(_lower_secondary_elem);
171 :
172 4325840 : _dof_to_weighted_gap[dof].first += _test[_i][_qp] * _qp_gap_nodal * _normals[_i];
173 :
174 2162920 : if (_normalize_c)
175 0 : _dof_to_weighted_gap[dof].second += _test[_i][_qp] * _qp_factor;
176 2162920 : }
177 :
178 : void
179 1712 : ComputeWeightedGapCartesianLMMechanicalContact::timestepSetup()
180 : {
181 : _dof_to_old_normal_vector.clear();
182 :
183 12208 : for (auto & map_pr : _dof_to_normal_vector)
184 : _dof_to_old_normal_vector.emplace(map_pr);
185 1712 : }
186 :
187 : void
188 24508 : ComputeWeightedGapCartesianLMMechanicalContact::residualSetup()
189 : {
190 : _dof_to_weighted_gap.clear();
191 : _dof_to_normal_vector.clear();
192 : _dof_to_tangent_vectors.clear();
193 24508 : }
194 :
195 : void
196 6678 : ComputeWeightedGapCartesianLMMechanicalContact::jacobianSetup()
197 : {
198 6678 : residualSetup();
199 6678 : }
200 :
201 : void
202 986286 : ComputeWeightedGapCartesianLMMechanicalContact::computeResidual(const Moose::MortarType mortar_type)
203 : {
204 986286 : if (mortar_type != Moose::MortarType::Lower)
205 : return;
206 :
207 : for (const auto i : index_range(_lm_vars))
208 : {
209 : mooseAssert(_lm_vars[i], "LM variable is null");
210 : libmesh_ignore(i);
211 : }
212 :
213 1127598 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
214 : {
215 798836 : computeQpProperties();
216 2961756 : for (_i = 0; _i < _test.size(); ++_i)
217 : {
218 2162920 : computeQpIProperties();
219 :
220 : // Get the _dof_to_weighted_gap map
221 : const DofObject * dof =
222 2162920 : _lm_vars[0]->isNodal()
223 2162920 : ? static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i))
224 0 : : static_cast<const DofObject *>(_lower_secondary_elem);
225 : // We do not interpolate geometry, so just match the local node _i with the corresponding _i
226 2162920 : _dof_to_normal_vector[dof] = _normals[_i];
227 2162920 : const auto & nodal_tangents = amg().getNodalTangents(*_lower_secondary_elem);
228 2162920 : _dof_to_tangent_vectors[dof][0] = nodal_tangents[0][_i];
229 2162920 : _dof_to_tangent_vectors[dof][1] = nodal_tangents[1][_i];
230 : }
231 : }
232 : }
233 :
234 : void
235 390066 : ComputeWeightedGapCartesianLMMechanicalContact::computeJacobian(const Moose::MortarType mortar_type)
236 : {
237 : // During "computeResidual" and "computeJacobian" we are actually just computing properties on the
238 : // mortar segment element mesh. We are *not* actually assembling into the residual/Jacobian. For
239 : // the zero-penetration constraint, the property of interest is the map from node to weighted gap.
240 : // Computation of the properties proceeds identically for residual and Jacobian evaluation hence
241 : // why we simply call computeResidual here. We will assemble into the residual/Jacobian later from
242 : // the post() method
243 390066 : computeResidual(mortar_type);
244 390066 : }
245 :
246 : void
247 966 : ComputeWeightedGapCartesianLMMechanicalContact::post()
248 : {
249 966 : Moose::Mortar::Contact::communicateGaps(
250 966 : _dof_to_weighted_gap, _mesh, _nodal, _normalize_c, _communicator, false);
251 :
252 7576 : for (const auto & pr : _dof_to_weighted_gap)
253 : {
254 6610 : if (pr.first->processor_id() != this->processor_id())
255 212 : continue;
256 :
257 6398 : _weighted_gap_ptr = &pr.second.first;
258 6398 : _normalization_ptr = &pr.second.second;
259 :
260 6398 : enforceConstraintOnDof(pr.first);
261 : }
262 966 : }
263 :
264 : void
265 18208 : ComputeWeightedGapCartesianLMMechanicalContact::incorrectEdgeDroppingPost(
266 : const std::unordered_set<const Node *> & inactive_lm_nodes)
267 : {
268 18208 : Moose::Mortar::Contact::communicateGaps(
269 18208 : _dof_to_weighted_gap, _mesh, _nodal, _normalize_c, _communicator, false);
270 :
271 68280 : for (const auto & pr : _dof_to_weighted_gap)
272 : {
273 50072 : if ((inactive_lm_nodes.find(static_cast<const Node *>(pr.first)) != inactive_lm_nodes.end()) ||
274 : (pr.first->processor_id() != this->processor_id()))
275 0 : continue;
276 :
277 50072 : _weighted_gap_ptr = &pr.second.first;
278 50072 : _normalization_ptr = &pr.second.second;
279 :
280 50072 : enforceConstraintOnDof(pr.first);
281 : }
282 18208 : }
283 :
284 : void
285 56470 : ComputeWeightedGapCartesianLMMechanicalContact::enforceConstraintOnDof(const DofObject * const dof)
286 : {
287 56470 : const auto & weighted_gap = *_weighted_gap_ptr;
288 56470 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
289 :
290 56470 : const auto dof_index_x = dof->dof_number(_sys.number(), _lm_vars[0]->number(), 0);
291 56470 : const auto dof_index_y = dof->dof_number(_sys.number(), _lm_vars[1]->number(), 0);
292 56470 : const Real scaling_factor_x = _lm_vars[0]->scalingFactor();
293 56470 : const Real scaling_factor_y = _lm_vars[1]->scalingFactor();
294 : Real scaling_factor_z = 1;
295 :
296 56470 : ADReal lm_x = (*_sys.currentSolution())(dof_index_x);
297 56470 : ADReal lm_y = (*_sys.currentSolution())(dof_index_y);
298 :
299 : Moose::derivInsert(lm_x.derivatives(), dof_index_x, 1.);
300 : Moose::derivInsert(lm_y.derivatives(), dof_index_y, 1.);
301 :
302 : dof_id_type dof_index_z(-1);
303 : ADReal lm_z;
304 56470 : if (_has_disp_z)
305 : {
306 2688 : dof_index_z = dof->dof_number(_sys.number(), _lm_vars[2]->number(), 0);
307 2688 : lm_z = (*_sys.currentSolution())(dof_index_z);
308 : Moose::derivInsert(lm_z.derivatives(), dof_index_z, 1.);
309 2688 : scaling_factor_z = _lm_vars[2]->scalingFactor();
310 : }
311 :
312 : ADReal normal_pressure_value =
313 56470 : lm_x * _dof_to_normal_vector[dof](0) + lm_y * _dof_to_normal_vector[dof](1);
314 : ADReal tangential_pressure_value =
315 56470 : lm_x * _dof_to_tangent_vectors[dof][0](0) + lm_y * _dof_to_tangent_vectors[dof][0](1);
316 :
317 : ADReal tangential_pressure_value_dir;
318 :
319 56470 : if (_has_disp_z)
320 : {
321 2688 : normal_pressure_value += lm_z * _dof_to_normal_vector[dof](2);
322 2688 : tangential_pressure_value += lm_z * _dof_to_tangent_vectors[dof][0](2);
323 2688 : tangential_pressure_value_dir = lm_x * _dof_to_tangent_vectors[dof][1](0) +
324 2688 : lm_y * _dof_to_tangent_vectors[dof][1](1) +
325 2688 : lm_z * _dof_to_tangent_vectors[dof][1](2);
326 : }
327 :
328 56470 : ADReal normal_dof_residual = std::min(normal_pressure_value, weighted_gap * c);
329 56470 : ADReal tangential_dof_residual = tangential_pressure_value;
330 :
331 : // Get index for normal constraint.
332 : // We do this to get a decent Jacobian structure, which is key for the use of iterative solvers.
333 : // Using old normal vector to avoid changes in the Jacobian structure within one time step
334 :
335 : Real ny, nz;
336 : // Intially, use the current normal vector
337 56470 : if (_dof_to_old_normal_vector[dof].norm() < TOLERANCE)
338 : {
339 3087 : ny = _dof_to_normal_vector[dof](1);
340 3087 : nz = _dof_to_normal_vector[dof](2);
341 : }
342 : else
343 : {
344 53383 : ny = _dof_to_old_normal_vector[dof](1);
345 53383 : nz = _dof_to_old_normal_vector[dof](2);
346 : }
347 : unsigned int component_normal = 0;
348 56470 : if (std::abs(ny) > 0.57735)
349 : component_normal = 1;
350 6398 : else if (std::abs(nz) > 0.57735)
351 : component_normal = 2;
352 :
353 : libmesh_ignore(component_normal);
354 :
355 112940 : addResidualsAndJacobian(
356 : _assembly,
357 56470 : std::array<ADReal, 1>{{normal_dof_residual}},
358 112940 : std::array<dof_id_type, 1>{{component_normal == 0
359 52760 : ? dof_index_x
360 : : (component_normal == 1 ? dof_index_y : dof_index_z)}},
361 : component_normal == 0 ? scaling_factor_x
362 : : (component_normal == 1 ? scaling_factor_y : scaling_factor_z));
363 :
364 56470 : addResidualsAndJacobian(
365 : _assembly,
366 56470 : std::array<ADReal, 1>{{tangential_dof_residual}},
367 112940 : std::array<dof_id_type, 1>{
368 56470 : {(component_normal == 0 || component_normal == 2) ? dof_index_y : dof_index_x}},
369 56470 : (component_normal == 0 || component_normal == 2) ? scaling_factor_y : scaling_factor_x);
370 :
371 56470 : if (_has_disp_z)
372 5376 : addResidualsAndJacobian(
373 : _assembly,
374 2688 : std::array<ADReal, 1>{{tangential_pressure_value_dir}},
375 2688 : std::array<dof_id_type, 1>{
376 : {(component_normal == 0 || component_normal == 1) ? dof_index_z : dof_index_x}},
377 : (component_normal == 0 || component_normal == 1) ? scaling_factor_z : scaling_factor_x);
378 172098 : }
|
