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