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