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 "RANFSNormalMechanicalContact.h"
11 : #include "PenetrationLocator.h"
12 : #include "PenetrationInfo.h"
13 : #include "SystemBase.h"
14 : #include "Assembly.h"
15 : #include "NearestNodeLocator.h"
16 : #include "MooseVariableFE.h"
17 :
18 : #include "libmesh/numeric_vector.h"
19 : #include "libmesh/enum_fe_family.h"
20 :
21 : registerMooseObject("ContactApp", RANFSNormalMechanicalContact);
22 :
23 : InputParameters
24 374 : RANFSNormalMechanicalContact::validParams()
25 : {
26 374 : InputParameters params = NodeFaceConstraint::validParams();
27 374 : params.set<bool>("use_displaced_mesh") = true;
28 :
29 748 : MooseEnum component("x=0 y=1 z=2");
30 748 : params.addRequiredParam<MooseEnum>(
31 : "component", component, "The force component constraint that this object is supplying");
32 748 : params.addRequiredCoupledVar(
33 : "displacements",
34 : "The displacements appropriate for the simulation geometry and coordinate system");
35 748 : params.addParam<bool>("ping_pong_protection",
36 748 : false,
37 : "Whether to protect against ping-ponging, e.g. the oscillation of the "
38 : "secondary node between two "
39 : "different primary faces, by tying the secondary node to the "
40 : "edge between the involved primary faces");
41 748 : params.addParam<Real>(
42 : "normal_smoothing_distance",
43 : "Distance from edge in parametric coordinates over which to smooth contact normal");
44 374 : params.addClassDescription("Applies the Reduced Active Nonlinear Function Set scheme in which "
45 : "the secondary node's non-linear residual function is replaced by the "
46 : "zero penetration constraint equation when the constraint is active");
47 374 : return params;
48 374 : }
49 :
50 196 : RANFSNormalMechanicalContact::RANFSNormalMechanicalContact(const InputParameters & parameters)
51 : : NodeFaceConstraint(parameters),
52 196 : _component(getParam<MooseEnum>("component")),
53 196 : _mesh_dimension(_mesh.dimension()),
54 196 : _residual_copy(_sys.residualGhosted()),
55 196 : _dof_number_to_value(coupledComponents("displacements")),
56 196 : _disp_coupling(coupledComponents("displacements")),
57 784 : _ping_pong_protection(getParam<bool>("ping_pong_protection"))
58 : {
59 : // modern parameter scheme for displacements
60 1176 : for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
61 : {
62 392 : _vars.push_back(coupled("displacements", i));
63 784 : _var_objects.push_back(getVar("displacements", i));
64 : }
65 :
66 588 : for (auto & var : _var_objects)
67 392 : if (var->feType().family != LAGRANGE)
68 0 : mooseError("This object only works when the displacement variables use a Lagrange basis");
69 :
70 196 : if (_vars.size() != _mesh_dimension)
71 0 : mooseError("The number of displacement variables does not match the mesh dimension!");
72 :
73 392 : if (parameters.isParamValid("normal_smoothing_distance"))
74 216 : _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
75 196 : }
76 :
77 : void
78 196 : RANFSNormalMechanicalContact::initialSetup()
79 : {
80 196 : auto system_coupling_matrix = _subproblem.couplingMatrix(_sys.number());
81 :
82 588 : for (MooseIndex(_vars) i = 0; i < _vars.size(); ++i)
83 1176 : for (MooseIndex(_vars) j = 0; j < _vars.size(); ++j)
84 784 : _disp_coupling(i, j) = (*system_coupling_matrix)(_vars[i], _vars[j]);
85 196 : }
86 :
87 : void
88 1638 : RANFSNormalMechanicalContact::timestepSetup()
89 : {
90 : _node_to_primary_elem_sequence.clear();
91 : _ping_pong_secondary_node_to_primary_node.clear();
92 1638 : }
93 :
94 : void
95 44830 : RANFSNormalMechanicalContact::residualSetup()
96 : {
97 : _node_to_contact_lm.clear();
98 : _node_to_tied_lm.clear();
99 44830 : NodeFaceConstraint::residualSetup();
100 44830 : }
101 :
102 : bool
103 61954 : RANFSNormalMechanicalContact::overwriteSecondaryResidual()
104 : {
105 61954 : if (_tie_nodes)
106 : return true;
107 : else
108 61954 : return _largest_component == static_cast<unsigned int>(_component);
109 : }
110 :
111 : bool
112 93678 : RANFSNormalMechanicalContact::shouldApply()
113 : {
114 : std::map<dof_id_type, PenetrationInfo *>::iterator found =
115 93678 : _penetration_locator._penetration_info.find(_current_node->id());
116 93678 : if (found != _penetration_locator._penetration_info.end())
117 : {
118 93678 : _pinfo = found->second;
119 93678 : if (_pinfo)
120 : {
121 : // We overwrite the secondary residual when constraints are active so we cannot use the
122 : // residual copy for determining the Lagrange multiplier when computing the Jacobian
123 93678 : if (!_subproblem.currentlyComputingJacobian())
124 : {
125 : // Build up residual vector corresponding to contact forces
126 : _res_vec.zero();
127 249984 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
128 : {
129 166656 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
130 166656 : _res_vec(i) = _residual_copy(dof_number) / _var_objects[i]->scalingFactor();
131 : }
132 :
133 83328 : _node_to_contact_lm.insert(std::make_pair(_current_node->id(), _res_vec * _pinfo->_normal));
134 83328 : _node_to_tied_lm.insert(std::make_pair(_current_node->id(), _res_vec(_component)));
135 : }
136 : else
137 : {
138 : std::vector<dof_id_type> cols;
139 : std::vector<Number> values;
140 :
141 31050 : for (auto & d_to_v : _dof_number_to_value)
142 : d_to_v.clear();
143 :
144 : mooseAssert(_vars.size() == _dof_number_to_value.size() &&
145 : _vars.size() == _var_objects.size(),
146 : "Somehow the sizes of our variable containers got out of sync");
147 31050 : for (MooseIndex(_var_objects) i = 0; i < _var_objects.size(); ++i)
148 : {
149 20700 : auto secondary_dof_number = _current_node->dof_number(0, _vars[i], 0);
150 :
151 20700 : _jacobian->get_row(secondary_dof_number, cols, values);
152 : mooseAssert(cols.size() == values.size(),
153 : "The size of the dof container and value container are different");
154 :
155 1027564 : for (MooseIndex(cols) j = 0; j < cols.size(); ++j)
156 : _dof_number_to_value[i].insert(
157 1006864 : std::make_pair(cols[j], values[j] / _var_objects[i]->scalingFactor()));
158 : }
159 : }
160 :
161 : mooseAssert(_node_to_contact_lm.find(_current_node->id()) != _node_to_contact_lm.end(),
162 : "The node " << _current_node->id()
163 : << " should map to a contact lagrange multiplier");
164 : mooseAssert(_node_to_tied_lm.find(_current_node->id()) != _node_to_tied_lm.end(),
165 : "The node " << _current_node->id()
166 : << " should map to a tied lagrange multiplier");
167 93678 : _contact_lm = _node_to_contact_lm[_current_node->id()];
168 93678 : _tied_lm = _node_to_tied_lm[_current_node->id()];
169 :
170 : // Check to see whether we've locked a ping-ponging node
171 93678 : if (_ping_pong_secondary_node_to_primary_node.find(_current_node->id()) ==
172 : _ping_pong_secondary_node_to_primary_node.end())
173 : {
174 93678 : if (_contact_lm > -_pinfo->_distance)
175 : {
176 : // Ok, our math is telling us we should apply the constraint, but what if we are
177 : // ping-ponging back and forth between different primary faces?
178 :
179 : // This only works for a basic line search! Write assertion here
180 57412 : if (_subproblem.computingNonlinearResid())
181 : {
182 6794 : auto & primary_elem_sequence = _node_to_primary_elem_sequence[_current_node->id()];
183 : mooseAssert(
184 : _current_primary == _pinfo->_elem,
185 : "The current primary element and the PenetrationInfo object's element should "
186 : "be the same");
187 6794 : primary_elem_sequence.push_back(_pinfo->_elem);
188 :
189 : // 5 here is a heuristic choice. In testing, generally speaking, if a node ping-pongs
190 : // back and forth 5 times then it will ping pong indefinitely. However, if it goes 3
191 : // times for example, it is not guaranteed to ping-pong indefinitely and the Newton
192 : // iteration may naturally resolve the correct face dofs that need to be involved in the
193 : // constraint.
194 1632 : if (_ping_pong_protection && primary_elem_sequence.size() >= 5 &&
195 0 : _pinfo->_elem == *(primary_elem_sequence.rbegin() + 2) &&
196 0 : _pinfo->_elem == *(primary_elem_sequence.rbegin() + 4) &&
197 6794 : _pinfo->_elem != *(primary_elem_sequence.rbegin() + 1) &&
198 0 : *(primary_elem_sequence.rbegin() + 1) == *(primary_elem_sequence.rbegin() + 3))
199 : {
200 : // Ok we are ping-ponging. Determine the primary node
201 : auto primary_node =
202 0 : _penetration_locator._nearest_node.nearestNode(_current_node->id());
203 :
204 : // Sanity checks
205 : mooseAssert(_pinfo->_elem->get_node_index(primary_node) != libMesh::invalid_uint,
206 : "The primary node is not on the current element");
207 : mooseAssert((*(primary_elem_sequence.rbegin() + 1))->get_node_index(primary_node) !=
208 : libMesh::invalid_uint,
209 : "The primary node is not on the other ping-ponging element");
210 :
211 : _ping_pong_secondary_node_to_primary_node.insert(
212 0 : std::make_pair(_current_node->id(), primary_node));
213 : }
214 : }
215 : }
216 : else
217 : // We have not locked the node into contact nor is the gap smaller than the Lagrange
218 : // Multiplier so we should not apply
219 : return false;
220 : }
221 :
222 : // Determine whether we're going to apply the tied node equality constraint or the contact
223 : // inequality constraint
224 57412 : auto it = _ping_pong_secondary_node_to_primary_node.find(_current_node->id());
225 57412 : if (it != _ping_pong_secondary_node_to_primary_node.end())
226 : {
227 0 : _tie_nodes = true;
228 0 : _nearest_node = it->second;
229 0 : _primary_index = _current_primary->get_node_index(_nearest_node);
230 : mooseAssert(_primary_index != libMesh::invalid_uint,
231 : "nearest node not a node on the current primary element");
232 : }
233 : else
234 : {
235 57412 : _distance = _pinfo->_distance;
236 : // Do this to make sure constraint equation has a positive on the diagonal
237 57412 : if (_pinfo->_normal(_component) > 0)
238 42880 : _distance *= -1;
239 57412 : _tie_nodes = false;
240 :
241 : // The contact constraint is active -> we're going to use our linear solve to ensure that
242 : // the gap is driven to zero. We only have one zero-penetration constraint per node, so we
243 : // choose to apply the zero penetration constraint only to the displacement component with
244 : // the largest magnitude normal
245 57412 : auto largest_component_magnitude = std::abs(_pinfo->_normal(0));
246 57412 : _largest_component = 0;
247 114824 : for (MooseIndex(_mesh_dimension) i = 1; i < _mesh_dimension; ++i)
248 : {
249 57412 : auto component_magnitude = std::abs(_pinfo->_normal(i));
250 57412 : if (component_magnitude > largest_component_magnitude)
251 : {
252 : largest_component_magnitude = component_magnitude;
253 57412 : _largest_component = i;
254 : }
255 : }
256 : }
257 :
258 57412 : return true;
259 : }
260 : }
261 :
262 : // If we're not applying the constraint then we can clear the node to primary elem sequence for
263 : // this node
264 0 : if (_node_to_primary_elem_sequence.find(_current_node->id()) !=
265 : _node_to_primary_elem_sequence.end())
266 0 : _node_to_primary_elem_sequence[_current_node->id()].clear();
267 : return false;
268 : }
269 :
270 : Real
271 254270 : RANFSNormalMechanicalContact::computeQpResidual(Moose::ConstraintType type)
272 : {
273 254270 : switch (type)
274 : {
275 50854 : case Moose::ConstraintType::Secondary:
276 : {
277 50854 : if (_tie_nodes)
278 0 : return (*_current_node - *_nearest_node)(_component);
279 : else
280 : {
281 50854 : if (_largest_component == static_cast<unsigned int>(_component))
282 : {
283 : mooseAssert(_pinfo->_normal(_component) != 0,
284 : "We should be selecting the largest normal component, hence it should be "
285 : "impossible for this normal component to be zero");
286 :
287 25427 : return _distance;
288 : }
289 :
290 : else
291 : // The normal points out of the primary face
292 25427 : return _contact_lm * -_pinfo->_normal(_component);
293 : }
294 : }
295 :
296 203416 : case Moose::ConstraintType::Primary:
297 : {
298 203416 : if (_tie_nodes)
299 : {
300 0 : if (_i == _primary_index)
301 0 : return _tied_lm;
302 : else
303 : return 0;
304 : }
305 : else
306 203416 : return _test_primary[_i][_qp] * _contact_lm * _pinfo->_normal(_component);
307 : }
308 :
309 : default:
310 : return 0;
311 : }
312 : }
313 :
314 : Real
315 241580 : RANFSNormalMechanicalContact::computeQpJacobian(Moose::ConstraintJacobianType type)
316 : {
317 241580 : switch (type)
318 : {
319 26116 : case Moose::ConstraintJacobianType::SecondarySecondary:
320 : {
321 26116 : if (_tie_nodes)
322 0 : return _phi_secondary[_j][_qp];
323 :
324 : // doing contact
325 : else
326 : {
327 : // corresponds to gap equation
328 26116 : if (_largest_component == static_cast<unsigned int>(_component))
329 : // _phi_secondary has been set such that it is 1 when _j corresponds to the degree of
330 : // freedom associated with the _current node and 0 otherwise
331 13058 : return std::abs(_pinfo->_normal(_component)) * _phi_secondary[_j][_qp];
332 :
333 : // corresponds to regular residual with Lagrange Multiplier applied
334 : else
335 : {
336 : Real ret_val = 0;
337 39174 : for (MooseIndex(_disp_coupling) i = 0; i < _disp_coupling.size(); ++i)
338 26116 : if (_disp_coupling(_component, i))
339 : {
340 : mooseAssert(
341 : _dof_number_to_value[i].find(_connected_dof_indices[_j]) !=
342 : _dof_number_to_value[i].end(),
343 : "The connected dof index is not found in the _dof_number_to_value container. "
344 : "This must mean that insufficient sparsity was allocated");
345 25932 : ret_val += -_pinfo->_normal(_component) * _pinfo->_normal(i) *
346 25932 : _dof_number_to_value[i][_connected_dof_indices[_j]];
347 : }
348 13058 : return ret_val;
349 : }
350 : }
351 : }
352 :
353 22200 : case Moose::ConstraintJacobianType::SecondaryPrimary:
354 : {
355 22200 : if (_tie_nodes)
356 : {
357 0 : if (_primary_index == _j)
358 : return -1;
359 :
360 : // We're tying the secondary node to only one node on the primary side (specified by
361 : // _primary_index). If the current _j doesn't correspond to that tied primary node, then the
362 : // secondary residual doesn't depend on it
363 : else
364 0 : return 0;
365 : }
366 : else
367 : {
368 22200 : if (_largest_component == static_cast<unsigned int>(_component))
369 11100 : return -std::abs(_pinfo->_normal(_component)) * _phi_primary[_j][_qp];
370 :
371 : // If we're not applying the gap constraint equation on this _component, then we're
372 : // applying a Lagrange multiplier, and consequently there is no dependence of the secondary
373 : // residual on the primary dofs because the Lagrange multiplier is only a functon of the
374 : // secondary residuals
375 : else
376 : return 0;
377 : }
378 : }
379 :
380 104464 : case Moose::ConstraintJacobianType::PrimarySecondary:
381 : {
382 104464 : if (_tie_nodes)
383 : {
384 0 : if (_i == _primary_index)
385 : {
386 : mooseAssert(_dof_number_to_value[_component].find(_connected_dof_indices[_j]) !=
387 : _dof_number_to_value[_component].end(),
388 : "The connected dof index is not found in the _dof_number_to_value container. "
389 : "This must mean that insufficient sparsity was allocated");
390 0 : return _dof_number_to_value[_component][_connected_dof_indices[_j]];
391 : }
392 :
393 : // We only apply the tied node Lagrange multiplier to the closest primary node
394 : else
395 : return 0;
396 : }
397 : else
398 : {
399 : Real ret_val = 0;
400 313392 : for (MooseIndex(_disp_coupling) i = 0; i < _disp_coupling.size(); ++i)
401 208928 : if (_disp_coupling(_component, i))
402 : {
403 : mooseAssert(
404 : _dof_number_to_value[i].find(_connected_dof_indices[_j]) !=
405 : _dof_number_to_value[i].end(),
406 : "The connected dof index is not found in the _dof_number_to_value container. "
407 : "This must mean that insufficient sparsity was allocated");
408 207456 : ret_val += _test_primary[_i][_qp] * _pinfo->_normal(_component) * _pinfo->_normal(i) *
409 207456 : _dof_number_to_value[i][_connected_dof_indices[_j]];
410 : }
411 104464 : return ret_val;
412 : }
413 : }
414 :
415 : // The only primary-primary dependence would come from the dependence of the normal and also
416 : // the location of the integration (quadrature) points. We assume (valid or not) that this
417 : // dependence is weak
418 : // case MooseConstraintJacobianType::PrimaryPrimary
419 :
420 : default:
421 : return 0;
422 : }
423 : }
424 :
425 : void
426 1008 : RANFSNormalMechanicalContact::computeSecondaryValue(NumericVector<Number> &)
427 : {
428 1008 : }
429 :
430 : Real
431 0 : RANFSNormalMechanicalContact::computeQpSecondaryValue()
432 : {
433 0 : mooseError(
434 : "We overrode commputeSecondaryValue so computeQpSecondaryValue should never get called");
435 : }
|
