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 : // MOOSE includes
11 : #include "MechanicalContactConstraint.h"
12 : #include "FEProblem.h"
13 : #include "DisplacedProblem.h"
14 : #include "AuxiliarySystem.h"
15 : #include "PenetrationLocator.h"
16 : #include "NearestNodeLocator.h"
17 : #include "SystemBase.h"
18 : #include "Assembly.h"
19 : #include "MooseMesh.h"
20 : #include "MathUtils.h"
21 : #include "AugmentedLagrangianContactProblem.h"
22 : #include "Executioner.h"
23 : #include "AddVariableAction.h"
24 : #include "ContactLineSearchBase.h"
25 : #include "ContactAction.h"
26 :
27 : #include "libmesh/string_to_enum.h"
28 : #include "libmesh/sparse_matrix.h"
29 :
30 : registerMooseObject("ContactApp", MechanicalContactConstraint);
31 :
32 : const unsigned int MechanicalContactConstraint::_no_iterations = 0;
33 :
34 : InputParameters
35 7072 : MechanicalContactConstraint::validParams()
36 : {
37 7072 : InputParameters params = NodeFaceConstraint::validParams();
38 7072 : params += ContactAction::commonParameters();
39 :
40 14144 : params.addRequiredParam<BoundaryName>("boundary", "The primary boundary");
41 14144 : params.addParam<BoundaryName>("secondary", "The secondary boundary");
42 14144 : params.addRequiredParam<unsigned int>("component",
43 : "An integer corresponding to the direction "
44 : "the variable this constraint acts on. (0 for x, "
45 : "1 for y, 2 for z)");
46 :
47 14144 : params.addCoupledVar(
48 : "displacements",
49 : "The displacements appropriate for the simulation geometry and coordinate system");
50 :
51 14144 : params.addCoupledVar("secondary_gap_offset", "offset to the gap distance from secondary side");
52 14144 : params.addCoupledVar("mapped_primary_gap_offset",
53 : "offset to the gap distance mapped from primary side");
54 14144 : params.addRequiredCoupledVar("nodal_area", "The nodal area");
55 :
56 7072 : params.set<bool>("use_displaced_mesh") = true;
57 14144 : params.addParam<Real>(
58 : "penalty",
59 14144 : 1e8,
60 : "The penalty to apply. This can vary depending on the stiffness of your materials");
61 14144 : params.addParam<Real>("penalty_multiplier",
62 14144 : 1.0,
63 : "The growth factor for the penalty applied at the end of each augmented "
64 : "Lagrange update iteration");
65 14144 : params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
66 14144 : params.addParam<Real>("tangential_tolerance",
67 : "Tangential distance to extend edges of contact surfaces");
68 14144 : params.addParam<Real>(
69 14144 : "capture_tolerance", 0, "Normal distance from surface within which nodes are captured");
70 :
71 14144 : params.addParam<Real>("tension_release",
72 14144 : 0.0,
73 : "Tension release threshold. A node in contact "
74 : "will not be released if its tensile load is below "
75 : "this value. No tension release if negative.");
76 :
77 14144 : params.addParam<bool>(
78 : "normalize_penalty",
79 14144 : false,
80 : "Whether to normalize the penalty parameter with the nodal area for penalty contact.");
81 14144 : params.addParam<bool>(
82 : "primary_secondary_jacobian",
83 14144 : true,
84 : "Whether to include jacobian entries coupling primary and secondary nodes.");
85 14144 : params.addParam<bool>(
86 : "connected_secondary_nodes_jacobian",
87 14144 : true,
88 : "Whether to include jacobian entries coupling nodes connected to secondary nodes.");
89 14144 : params.addParam<bool>("non_displacement_variables_jacobian",
90 14144 : true,
91 : "Whether to include jacobian entries coupling with variables that are not "
92 : "displacement variables.");
93 14144 : params.addParam<unsigned int>("stick_lock_iterations",
94 14144 : std::numeric_limits<unsigned int>::max(),
95 : "Number of times permitted to switch between sticking and slipping "
96 : "in a solution before locking node in a sticked state.");
97 14144 : params.addParam<Real>("stick_unlock_factor",
98 14144 : 1.5,
99 : "Factor by which frictional capacity must be "
100 : "exceeded to permit stick-locked node to slip "
101 : "again.");
102 14144 : params.addParam<Real>("al_penetration_tolerance",
103 : "The tolerance of the penetration for augmented Lagrangian method.");
104 14144 : params.addParam<Real>("al_incremental_slip_tolerance",
105 : "The tolerance of the incremental slip for augmented Lagrangian method.");
106 :
107 14144 : params.addParam<Real>("al_frictional_force_tolerance",
108 : "The tolerance of the frictional force for augmented Lagrangian method.");
109 14144 : params.addParam<bool>(
110 14144 : "print_contact_nodes", false, "Whether to print the number of nodes in contact.");
111 :
112 7072 : params.addClassDescription(
113 : "Apply non-penetration constraints on the mechanical deformation "
114 : "using a node on face, primary/secondary algorithm, and multiple options "
115 : "for the physical behavior on the interface and the mathematical "
116 : "formulation for constraint enforcement");
117 :
118 7072 : return params;
119 0 : }
120 :
121 : Threads::spin_mutex MechanicalContactConstraint::_contact_set_mutex;
122 :
123 5011 : MechanicalContactConstraint::MechanicalContactConstraint(const InputParameters & parameters)
124 : : NodeFaceConstraint(parameters),
125 5011 : _displaced_problem(parameters.get<FEProblemBase *>("_fe_problem_base")->getDisplacedProblem()),
126 10022 : _component(getParam<unsigned int>("component")),
127 10022 : _model(getParam<MooseEnum>("model").getEnum<ContactModel>()),
128 10022 : _formulation(getParam<MooseEnum>("formulation").getEnum<ContactFormulation>()),
129 10022 : _normalize_penalty(getParam<bool>("normalize_penalty")),
130 10022 : _penalty(getParam<Real>("penalty")),
131 10022 : _penalty_multiplier(getParam<Real>("penalty_multiplier")),
132 10022 : _friction_coefficient(getParam<Real>("friction_coefficient")),
133 10022 : _tension_release(getParam<Real>("tension_release")),
134 10022 : _capture_tolerance(getParam<Real>("capture_tolerance")),
135 10022 : _stick_lock_iterations(getParam<unsigned int>("stick_lock_iterations")),
136 10022 : _stick_unlock_factor(getParam<Real>("stick_unlock_factor")),
137 5011 : _update_stateful_data(true),
138 5011 : _residual_copy(_sys.residualGhosted()),
139 5011 : _mesh_dimension(_mesh.dimension()),
140 5011 : _vars(3, libMesh::invalid_uint),
141 5011 : _var_objects(3, nullptr),
142 5011 : _has_secondary_gap_offset(isCoupled("secondary_gap_offset")),
143 5011 : _secondary_gap_offset_var(_has_secondary_gap_offset ? getVar("secondary_gap_offset", 0)
144 : : nullptr),
145 5011 : _has_mapped_primary_gap_offset(isCoupled("mapped_primary_gap_offset")),
146 5011 : _mapped_primary_gap_offset_var(
147 5011 : _has_mapped_primary_gap_offset ? getVar("mapped_primary_gap_offset", 0) : nullptr),
148 5011 : _nodal_area_var(getVar("nodal_area", 0)),
149 5011 : _aux_system(_nodal_area_var->sys()),
150 5011 : _aux_solution(_aux_system.currentSolution()),
151 10022 : _primary_secondary_jacobian(getParam<bool>("primary_secondary_jacobian")),
152 10022 : _connected_secondary_nodes_jacobian(getParam<bool>("connected_secondary_nodes_jacobian")),
153 10022 : _non_displacement_vars_jacobian(getParam<bool>("non_displacement_variables_jacobian")),
154 5011 : _contact_linesearch(dynamic_cast<ContactLineSearchBase *>(_subproblem.getLineSearch())),
155 10022 : _print_contact_nodes(getParam<bool>("print_contact_nodes")),
156 5011 : _augmented_lagrange_problem(
157 5011 : dynamic_cast<AugmentedLagrangianContactProblemInterface *>(&_fe_problem)),
158 5011 : _lagrangian_iteration_number(_augmented_lagrange_problem
159 5011 : ? _augmented_lagrange_problem->getLagrangianIterationNumber()
160 10022 : : _no_iterations)
161 : {
162 5011 : _overwrite_secondary_residual = false;
163 :
164 10022 : if (isParamValid("displacements"))
165 : {
166 : // modern parameter scheme for displacements
167 32964 : for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
168 : {
169 11471 : _vars[i] = coupled("displacements", i);
170 11471 : _var_objects[i] = getVar("displacements", i);
171 : }
172 : }
173 :
174 10022 : if (parameters.isParamValid("tangential_tolerance"))
175 10206 : _penetration_locator.setTangentialTolerance(getParam<Real>("tangential_tolerance"));
176 :
177 10022 : if (parameters.isParamValid("normal_smoothing_distance"))
178 1425 : _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
179 :
180 10022 : if (parameters.isParamValid("normal_smoothing_method"))
181 0 : _penetration_locator.setNormalSmoothingMethod(
182 : parameters.get<std::string>("normal_smoothing_method"));
183 :
184 5011 : if (_formulation == ContactFormulation::TANGENTIAL_PENALTY && _model != ContactModel::COULOMB)
185 0 : mooseError("The 'tangential_penalty' formulation can only be used with the 'coulomb' model");
186 :
187 5011 : if (_model == ContactModel::GLUED)
188 644 : _penetration_locator.setUpdate(false);
189 :
190 5011 : if (_friction_coefficient < 0)
191 0 : mooseError("The friction coefficient must be nonnegative");
192 :
193 : // set _penalty_tangential to the value of _penalty for now
194 5011 : _penalty_tangential = _penalty;
195 :
196 5011 : if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
197 : {
198 564 : if (_model == ContactModel::GLUED)
199 0 : mooseError("The Augmented Lagrangian contact formulation does not support GLUED case.");
200 :
201 564 : if (!_augmented_lagrange_problem)
202 0 : mooseError("The Augmented Lagrangian contact formulation must use "
203 : "AugmentedLagrangianContactProblem.");
204 :
205 1128 : if (!parameters.isParamValid("al_penetration_tolerance"))
206 0 : mooseError("For Augmented Lagrangian contact, al_penetration_tolerance must be provided.");
207 : else
208 564 : _al_penetration_tolerance = parameters.get<Real>("al_penetration_tolerance");
209 :
210 564 : if (_model != ContactModel::FRICTIONLESS)
211 : {
212 720 : if (!parameters.isParamValid("al_incremental_slip_tolerance") ||
213 720 : !parameters.isParamValid("al_frictional_force_tolerance"))
214 : {
215 0 : mooseError("For the Augmented Lagrangian frictional contact formualton, "
216 : "al_incremental_slip_tolerance and "
217 : "al_frictional_force_tolerance must be provided.");
218 : }
219 : else
220 : {
221 240 : _al_incremental_slip_tolerance = parameters.get<Real>("al_incremental_slip_tolerance");
222 240 : _al_frictional_force_tolerance = parameters.get<Real>("al_frictional_force_tolerance");
223 : }
224 : }
225 : }
226 5011 : }
227 :
228 : void
229 20450 : MechanicalContactConstraint::timestepSetup()
230 : {
231 20450 : if (_component == 0)
232 : {
233 9801 : updateContactStatefulData(/* beginning_of_step = */ true);
234 9801 : if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
235 213 : updateAugmentedLagrangianMultiplier(/* beginning_of_step = */ true);
236 :
237 9801 : _update_stateful_data = false;
238 :
239 9801 : if (_contact_linesearch)
240 17 : _contact_linesearch->reset();
241 : }
242 20450 : }
243 :
244 : void
245 59402 : MechanicalContactConstraint::jacobianSetup()
246 : {
247 59402 : if (_component == 0)
248 : {
249 28735 : if (_update_stateful_data)
250 19620 : updateContactStatefulData(/* beginning_of_step = */ false);
251 28735 : _update_stateful_data = true;
252 : }
253 59402 : }
254 :
255 : void
256 429 : MechanicalContactConstraint::updateAugmentedLagrangianMultiplier(bool beginning_of_step)
257 : {
258 3242 : for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
259 : {
260 2813 : if (!pinfo)
261 175 : continue;
262 :
263 2638 : const Node & node = _mesh.nodeRef(secondary_node_num);
264 2638 : if (node.n_comp(_sys.number(), _vars[_component]) < 1)
265 0 : continue;
266 :
267 2638 : const Real distance = pinfo->_normal * (pinfo->_closest_point - node) - gapOffset(node);
268 :
269 2638 : if (beginning_of_step && _model == ContactModel::COULOMB)
270 : {
271 420 : pinfo->_lagrange_multiplier_slip.zero();
272 420 : if (pinfo->isCaptured())
273 0 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
274 : }
275 :
276 2638 : if (pinfo->isCaptured())
277 : {
278 1516 : if (_model == ContactModel::FRICTIONLESS)
279 1026 : pinfo->_lagrange_multiplier -= getPenalty(node) * distance;
280 :
281 1516 : if (_model == ContactModel::COULOMB)
282 : {
283 490 : if (!beginning_of_step)
284 : {
285 490 : Real penalty = getPenalty(node);
286 : RealVectorValue pen_force_normal =
287 490 : penalty * (-distance) * pinfo->_normal + pinfo->_lagrange_multiplier * pinfo->_normal;
288 :
289 : // update normal lagrangian multiplier
290 490 : pinfo->_lagrange_multiplier += penalty * (-distance);
291 :
292 : // Frictional capacity
293 889 : const Real capacity(_friction_coefficient * (pen_force_normal * pinfo->_normal < 0
294 490 : ? -pen_force_normal * pinfo->_normal
295 : : 0));
296 :
297 : RealVectorValue tangential_inc_slip =
298 : pinfo->_incremental_slip -
299 : (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
300 :
301 490 : Real penalty_slip = getTangentialPenalty(node);
302 :
303 : RealVectorValue inc_pen_force_tangential =
304 490 : pinfo->_lagrange_multiplier_slip + penalty_slip * tangential_inc_slip;
305 :
306 : RealVectorValue tau_old = pinfo->_contact_force_old -
307 : pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old);
308 :
309 : RealVectorValue contact_force_tangential = inc_pen_force_tangential + tau_old;
310 490 : const Real tan_mag(contact_force_tangential.norm());
311 :
312 490 : if (tan_mag > capacity * (_al_frictional_force_tolerance + 1.0))
313 : {
314 63 : pinfo->_lagrange_multiplier_slip =
315 : -tau_old + capacity * contact_force_tangential / tan_mag;
316 63 : if (MooseUtils::absoluteFuzzyEqual(capacity, 0.0))
317 63 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
318 : else
319 0 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
320 : }
321 : else
322 : {
323 427 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
324 : pinfo->_lagrange_multiplier_slip += penalty_slip * tangential_inc_slip;
325 : }
326 : }
327 : }
328 : }
329 : }
330 429 : }
331 :
332 : bool
333 708 : MechanicalContactConstraint::AugmentedLagrangianContactConverged()
334 : {
335 : Real contactResidual = 0.0;
336 708 : unsigned int converged = 0;
337 :
338 4130 : for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
339 : {
340 3638 : if (!pinfo)
341 0 : continue;
342 :
343 3638 : const auto & node = _mesh.nodeRef(secondary_node_num);
344 :
345 : // Skip this pinfo if there are no DOFs on this node.
346 3638 : if (node.n_comp(_sys.number(), _vars[_component]) < 1)
347 0 : continue;
348 :
349 3638 : const Real distance = pinfo->_normal * (pinfo->_closest_point - node) - gapOffset(node);
350 :
351 3638 : if (pinfo->isCaptured())
352 : {
353 3638 : if (contactResidual < std::abs(distance))
354 : contactResidual = std::abs(distance);
355 :
356 : // penetration < tol
357 3638 : if (contactResidual > _al_penetration_tolerance)
358 : {
359 216 : converged = 1;
360 216 : break;
361 : }
362 :
363 3422 : if (_model == ContactModel::COULOMB)
364 : {
365 : RealVectorValue contact_force_normal((pinfo->_contact_force * pinfo->_normal) *
366 : pinfo->_normal);
367 : RealVectorValue contact_force_tangential(pinfo->_contact_force - contact_force_normal);
368 :
369 : RealVectorValue tangential_inc_slip =
370 : pinfo->_incremental_slip - (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
371 :
372 1274 : const Real tan_mag(contact_force_tangential.norm());
373 1274 : const Real tangential_inc_slip_mag = tangential_inc_slip.norm();
374 :
375 : RealVectorValue distance_vec =
376 1274 : (pinfo->_normal * (node - pinfo->_closest_point) + gapOffset(node)) * pinfo->_normal;
377 :
378 1274 : Real penalty = getPenalty(node);
379 : RealVectorValue pen_force_normal =
380 1274 : penalty * distance_vec + pinfo->_lagrange_multiplier * pinfo->_normal;
381 :
382 : // Frictional capacity
383 2275 : Real capacity(_friction_coefficient * (pen_force_normal * pinfo->_normal < 0
384 1274 : ? -pen_force_normal * pinfo->_normal
385 : : 0.0));
386 :
387 : // incremental slip <= tol for all pinfo_pair such that tan_mag < capacity
388 1274 : if (MooseUtils::absoluteFuzzyLessThan(tan_mag, capacity) &&
389 1001 : pinfo->_mech_status == PenetrationInfo::MS_STICKING)
390 : {
391 1001 : if (MooseUtils::absoluteFuzzyGreaterThan(tangential_inc_slip_mag,
392 : _al_incremental_slip_tolerance))
393 : {
394 0 : converged = 2;
395 0 : break;
396 : }
397 : }
398 :
399 : // for all pinfo_pair, tag_mag should be less than (1 + tol) * (capacity + tol)
400 1274 : if (tan_mag >
401 1274 : (1 + _al_frictional_force_tolerance) * (capacity + _al_frictional_force_tolerance))
402 : {
403 0 : converged = 3;
404 0 : break;
405 : }
406 : }
407 : }
408 : }
409 :
410 708 : _communicator.max(converged);
411 :
412 708 : if (converged == 1)
413 216 : _console << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied "
414 216 : << std::endl;
415 492 : else if (converged == 2)
416 0 : _console << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied "
417 0 : << std::endl;
418 492 : else if (converged == 3)
419 0 : _console << "The Augmented Lagrangian contact frictional force enforcement is NOT satisfied "
420 0 : << std::endl;
421 : else
422 : return true;
423 :
424 : return false;
425 : }
426 :
427 : void
428 29421 : MechanicalContactConstraint::updateContactStatefulData(bool beginning_of_step)
429 : {
430 253784 : for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
431 : {
432 224363 : if (!pinfo)
433 8355 : continue;
434 :
435 216008 : const Node & node = _mesh.nodeRef(secondary_node_num);
436 216008 : if (node.n_comp(_sys.number(), _vars[_component]) < 1)
437 0 : continue;
438 :
439 216008 : if (beginning_of_step)
440 : {
441 68439 : if (_app.getExecutioner()->lastSolveConverged())
442 : {
443 66577 : pinfo->_contact_force_old = pinfo->_contact_force;
444 66577 : pinfo->_accumulated_slip_old = pinfo->_accumulated_slip;
445 66577 : pinfo->_frictional_energy_old = pinfo->_frictional_energy;
446 66577 : pinfo->_mech_status_old = pinfo->_mech_status;
447 : }
448 1862 : else if (pinfo->_mech_status_old == PenetrationInfo::MS_NO_CONTACT &&
449 1675 : pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
450 : {
451 : // The penetration info object could be based on a bad state so delete it
452 1217 : delete pinfo;
453 1217 : pinfo = NULL;
454 1217 : continue;
455 : }
456 :
457 67222 : pinfo->_locked_this_step = 0;
458 67222 : pinfo->_stick_locked_this_step = 0;
459 67222 : pinfo->_starting_elem = pinfo->_elem;
460 67222 : pinfo->_starting_side_num = pinfo->_side_num;
461 67222 : pinfo->_starting_closest_point_ref = pinfo->_closest_point_ref;
462 : }
463 214791 : pinfo->_incremental_slip_prev_iter = pinfo->_incremental_slip;
464 : }
465 29421 : }
466 :
467 : bool
468 3291286 : MechanicalContactConstraint::shouldApply()
469 : {
470 : bool in_contact = false;
471 :
472 : std::map<dof_id_type, PenetrationInfo *>::iterator found =
473 3291286 : _penetration_locator._penetration_info.find(_current_node->id());
474 3291286 : if (found != _penetration_locator._penetration_info.end())
475 : {
476 3291286 : PenetrationInfo * pinfo = found->second;
477 3291286 : if (pinfo != NULL)
478 : {
479 3291286 : bool is_nonlinear = _subproblem.computingNonlinearResid();
480 :
481 : // This computes the contact force once per constraint, rather than once per quad point
482 : // and for both primary and secondary cases.
483 3291286 : if (_component == 0)
484 1403672 : computeContactForce(*_current_node, pinfo, is_nonlinear);
485 :
486 3291286 : if (pinfo->isCaptured())
487 : {
488 : in_contact = true;
489 2165165 : if (is_nonlinear)
490 : {
491 : Threads::spin_mutex::scoped_lock lock(_contact_set_mutex);
492 253107 : _current_contact_state.insert(_current_node->id());
493 : }
494 : }
495 : }
496 : }
497 :
498 3291286 : return in_contact;
499 : }
500 :
501 : void
502 1403672 : MechanicalContactConstraint::computeContactForce(const Node & node,
503 : PenetrationInfo * pinfo,
504 : bool update_contact_set)
505 : {
506 : // Build up residual vector
507 : RealVectorValue res_vec;
508 4694958 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
509 : {
510 3291286 : dof_id_type dof_number = node.dof_number(0, _vars[i], 0);
511 3291286 : res_vec(i) = _residual_copy(dof_number) / _var_objects[i]->scalingFactor();
512 : }
513 :
514 : RealVectorValue distance_vec(node - pinfo->_closest_point);
515 1403672 : if (distance_vec.norm() != 0)
516 1347896 : distance_vec += gapOffset(node) * pinfo->_normal * distance_vec.unit() * distance_vec.unit();
517 :
518 : const Real gap_size = -1.0 * pinfo->_normal * distance_vec;
519 :
520 : // This is for preventing an increment of pinfo->_locked_this_step for nodes that are
521 : // captured and released in this function
522 : bool newly_captured = false;
523 :
524 : // Capture nodes that are newly in contact
525 1403672 : if (update_contact_set && !pinfo->isCaptured() &&
526 : MooseUtils::absoluteFuzzyGreaterEqual(gap_size, 0.0, _capture_tolerance))
527 : {
528 : newly_captured = true;
529 : pinfo->capture();
530 :
531 : // Increment the lock count every time the node comes back into contact from not being in
532 : // contact.
533 14208 : if (_formulation == ContactFormulation::KINEMATIC ||
534 : _formulation == ContactFormulation::TANGENTIAL_PENALTY)
535 7720 : ++pinfo->_locked_this_step;
536 : }
537 :
538 1403672 : if (!pinfo->isCaptured())
539 501625 : return;
540 :
541 902047 : const Real penalty = getPenalty(node);
542 902047 : const Real penalty_slip = getTangentialPenalty(node);
543 :
544 : RealVectorValue pen_force(penalty * distance_vec);
545 :
546 902047 : switch (_model)
547 : {
548 627985 : case ContactModel::FRICTIONLESS:
549 627985 : switch (_formulation)
550 : {
551 350657 : case ContactFormulation::KINEMATIC:
552 350657 : pinfo->_contact_force = -pinfo->_normal * (pinfo->_normal * res_vec);
553 350657 : break;
554 :
555 249635 : case ContactFormulation::PENALTY:
556 249635 : pinfo->_contact_force = pinfo->_normal * (pinfo->_normal * pen_force);
557 249635 : break;
558 :
559 27693 : case ContactFormulation::AUGMENTED_LAGRANGE:
560 27693 : pinfo->_contact_force =
561 : (pinfo->_normal *
562 : (pinfo->_normal * (pen_force + pinfo->_lagrange_multiplier * pinfo->_normal)));
563 27693 : break;
564 :
565 0 : default:
566 0 : mooseError("Invalid contact formulation");
567 : break;
568 : }
569 627985 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
570 627985 : break;
571 106989 : case ContactModel::COULOMB:
572 106989 : switch (_formulation)
573 : {
574 24152 : case ContactFormulation::KINEMATIC:
575 : {
576 : // Frictional capacity
577 45512 : const Real capacity(_friction_coefficient *
578 24152 : (res_vec * pinfo->_normal > 0 ? res_vec * pinfo->_normal : 0));
579 :
580 : // Normal and tangential components of predictor force
581 24152 : pinfo->_contact_force = -res_vec;
582 : RealVectorValue contact_force_normal((pinfo->_contact_force * pinfo->_normal) *
583 : pinfo->_normal);
584 : RealVectorValue contact_force_tangential(pinfo->_contact_force - contact_force_normal);
585 :
586 : RealVectorValue tangential_inc_slip =
587 : pinfo->_incremental_slip -
588 : (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
589 :
590 : // Magnitude of tangential predictor force
591 24152 : const Real tan_mag(contact_force_tangential.norm());
592 24152 : const Real tangential_inc_slip_mag = tangential_inc_slip.norm();
593 24152 : const Real slip_tol = capacity / penalty;
594 24152 : pinfo->_slip_tol = slip_tol;
595 :
596 24152 : if ((tangential_inc_slip_mag > slip_tol || tan_mag > capacity) &&
597 22412 : (pinfo->_stick_locked_this_step < _stick_lock_iterations ||
598 0 : tan_mag > capacity * _stick_unlock_factor))
599 : {
600 22412 : if (pinfo->_stick_locked_this_step >= _stick_lock_iterations)
601 0 : pinfo->_stick_locked_this_step = 0;
602 :
603 : // Check for scenario where node has slipped too far in a previous iteration
604 : // for special treatment (only done if the incremental slip is non-trivial)
605 : bool slipped_too_far = false;
606 : RealVectorValue slip_inc_direction;
607 22412 : if (tangential_inc_slip_mag > slip_tol)
608 : {
609 : slip_inc_direction = tangential_inc_slip / tangential_inc_slip_mag;
610 : Real slip_dot_tang_force = slip_inc_direction * contact_force_tangential;
611 21650 : if (slip_dot_tang_force < capacity)
612 : slipped_too_far = true;
613 : }
614 :
615 : if (slipped_too_far) // slip back along slip increment
616 12520 : pinfo->_contact_force = contact_force_normal + capacity * slip_inc_direction;
617 : else
618 : {
619 9892 : if (tan_mag > 0) // slip along tangential force direction
620 9869 : pinfo->_contact_force =
621 : contact_force_normal + capacity * contact_force_tangential / tan_mag;
622 : else // treat as frictionless
623 23 : pinfo->_contact_force = contact_force_normal;
624 : }
625 22412 : if (capacity == 0)
626 2809 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
627 : else
628 19603 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
629 : }
630 : else
631 : {
632 1740 : if (pinfo->_mech_status != PenetrationInfo::MS_STICKING &&
633 : pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
634 449 : ++pinfo->_stick_locked_this_step;
635 1740 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
636 : }
637 : break;
638 : }
639 :
640 17783 : case ContactFormulation::PENALTY:
641 : {
642 17783 : distance_vec =
643 : pinfo->_incremental_slip +
644 17783 : (pinfo->_normal * (node - pinfo->_closest_point) + gapOffset(node)) * pinfo->_normal;
645 : pen_force = penalty * distance_vec;
646 :
647 : // Frictional capacity
648 30076 : const Real capacity(_friction_coefficient *
649 17783 : (pen_force * pinfo->_normal < 0 ? -pen_force * pinfo->_normal : 0));
650 :
651 : // Elastic predictor
652 17783 : pinfo->_contact_force =
653 : pen_force + (pinfo->_contact_force_old -
654 : pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old));
655 : RealVectorValue contact_force_normal((pinfo->_contact_force * pinfo->_normal) *
656 : pinfo->_normal);
657 : RealVectorValue contact_force_tangential(pinfo->_contact_force - contact_force_normal);
658 :
659 : // Tangential magnitude of elastic predictor
660 17783 : const Real tan_mag(contact_force_tangential.norm());
661 :
662 17783 : if (tan_mag > capacity)
663 : {
664 3181 : pinfo->_contact_force =
665 : contact_force_normal + capacity * contact_force_tangential / tan_mag;
666 3181 : if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
667 2461 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
668 : else
669 720 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
670 : }
671 : else
672 14602 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
673 : break;
674 : }
675 :
676 17598 : case ContactFormulation::AUGMENTED_LAGRANGE:
677 : {
678 17598 : distance_vec =
679 17598 : (pinfo->_normal * (node - pinfo->_closest_point) + gapOffset(node)) * pinfo->_normal;
680 :
681 : RealVectorValue contact_force_normal =
682 : penalty * distance_vec + pinfo->_lagrange_multiplier * pinfo->_normal;
683 :
684 : RealVectorValue tangential_inc_slip =
685 : pinfo->_incremental_slip -
686 : (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
687 :
688 : RealVectorValue contact_force_tangential =
689 : pinfo->_lagrange_multiplier_slip +
690 : (pinfo->_contact_force_old -
691 : pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old));
692 :
693 : RealVectorValue inc_pen_force_tangential = penalty_slip * tangential_inc_slip;
694 :
695 17598 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
696 7441 : pinfo->_contact_force =
697 : contact_force_normal + contact_force_tangential + inc_pen_force_tangential;
698 : else
699 10157 : pinfo->_contact_force = contact_force_normal + contact_force_tangential;
700 :
701 : break;
702 : }
703 :
704 47456 : case ContactFormulation::TANGENTIAL_PENALTY:
705 : {
706 : // Frictional capacity (kinematic formulation)
707 51614 : const Real capacity = _friction_coefficient * std::max(res_vec * pinfo->_normal, 0.0);
708 :
709 : // Normal component of contact force (kinematic formulation)
710 : RealVectorValue contact_force_normal((-res_vec * pinfo->_normal) * pinfo->_normal);
711 :
712 : // Predictor tangential component of contact force (penalty formulation)
713 : RealVectorValue inc_pen_force_tangential = penalty * pinfo->_incremental_slip;
714 : RealVectorValue contact_force_tangential =
715 : inc_pen_force_tangential +
716 : (pinfo->_contact_force_old -
717 : pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old));
718 :
719 : // Magnitude of tangential predictor
720 47456 : const Real tan_mag(contact_force_tangential.norm());
721 :
722 47456 : if (tan_mag > capacity)
723 : {
724 38201 : pinfo->_contact_force =
725 : contact_force_normal + capacity * contact_force_tangential / tan_mag;
726 38201 : if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
727 5521 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
728 : else
729 32680 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
730 : }
731 : else
732 : {
733 9255 : pinfo->_contact_force = contact_force_normal + contact_force_tangential;
734 9255 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
735 : }
736 : break;
737 : }
738 :
739 0 : default:
740 0 : mooseError("Invalid contact formulation");
741 : break;
742 : }
743 : break;
744 :
745 167073 : case ContactModel::GLUED:
746 167073 : switch (_formulation)
747 : {
748 : case ContactFormulation::KINEMATIC:
749 116033 : pinfo->_contact_force = -res_vec;
750 116033 : break;
751 :
752 51040 : case ContactFormulation::PENALTY:
753 51040 : pinfo->_contact_force = pen_force;
754 51040 : break;
755 :
756 0 : case ContactFormulation::AUGMENTED_LAGRANGE:
757 0 : pinfo->_contact_force =
758 0 : pen_force + pinfo->_lagrange_multiplier * distance_vec / distance_vec.norm();
759 0 : break;
760 :
761 0 : default:
762 0 : mooseError("Invalid contact formulation");
763 : break;
764 : }
765 167073 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
766 167073 : break;
767 :
768 0 : default:
769 0 : mooseError("Invalid or unavailable contact model");
770 : break;
771 : }
772 :
773 : // Release
774 109105 : if (update_contact_set && _model != ContactModel::GLUED && pinfo->isCaptured() &&
775 982272 : !newly_captured && _tension_release >= 0.0 &&
776 77183 : (_contact_linesearch ? true : pinfo->_locked_this_step < 2))
777 : {
778 77120 : const Real contact_pressure = -(pinfo->_normal * pinfo->_contact_force) / nodalArea(node);
779 77120 : if (-contact_pressure >= _tension_release)
780 : {
781 : pinfo->release();
782 : pinfo->_contact_force.zero();
783 : }
784 : }
785 : }
786 :
787 : Real
788 33113 : MechanicalContactConstraint::computeQpSecondaryValue()
789 : {
790 33113 : return _u_secondary[_qp];
791 : }
792 :
793 : Real
794 20620598 : MechanicalContactConstraint::computeQpResidual(Moose::ConstraintType type)
795 : {
796 20620598 : PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
797 20620598 : Real resid = pinfo->_contact_force(_component);
798 20620598 : switch (type)
799 : {
800 1943405 : case Moose::Secondary:
801 1943405 : if (_formulation == ContactFormulation::KINEMATIC)
802 : {
803 1112652 : RealVectorValue distance_vec(*_current_node - pinfo->_closest_point);
804 1112652 : if (distance_vec.norm() != 0)
805 1054878 : distance_vec += gapOffset(*_current_node) * pinfo->_normal * distance_vec.unit() *
806 2109756 : distance_vec.unit();
807 :
808 1112652 : const Real penalty = getPenalty(*_current_node);
809 : RealVectorValue pen_force(penalty * distance_vec);
810 :
811 1112652 : if (_model == ContactModel::FRICTIONLESS)
812 786753 : resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
813 325899 : else if (_model == ContactModel::COULOMB)
814 : {
815 : distance_vec = distance_vec - pinfo->_incremental_slip;
816 40214 : pen_force = penalty * distance_vec;
817 40214 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
818 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
819 36846 : resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
820 : else
821 3368 : resid += pen_force(_component);
822 : }
823 285685 : else if (_model == ContactModel::GLUED)
824 285685 : resid += pen_force(_component);
825 : }
826 830753 : else if (_formulation == ContactFormulation::TANGENTIAL_PENALTY &&
827 80958 : _model == ContactModel::COULOMB)
828 : {
829 80958 : RealVectorValue distance_vec = (pinfo->_normal * (*_current_node - pinfo->_closest_point) +
830 80958 : gapOffset(*_current_node)) *
831 : pinfo->_normal;
832 :
833 80958 : const Real penalty = getPenalty(*_current_node);
834 : RealVectorValue pen_force(penalty * distance_vec);
835 80958 : resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
836 : }
837 1943405 : return _test_secondary[_i][_qp] * resid;
838 :
839 18677193 : case Moose::Primary:
840 18677193 : return _test_primary[_i][_qp] * -resid;
841 : }
842 :
843 : return 0.0;
844 : }
845 :
846 : Real
847 73522813 : MechanicalContactConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
848 : {
849 73522813 : PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
850 :
851 73522813 : const Real penalty = getPenalty(*_current_node);
852 73522813 : const Real penalty_slip = getTangentialPenalty(*_current_node);
853 :
854 73522813 : switch (type)
855 : {
856 0 : default:
857 0 : mooseError("Unhandled ConstraintJacobianType");
858 :
859 2564905 : case Moose::SecondarySecondary:
860 2564905 : switch (_model)
861 : {
862 2017787 : case ContactModel::FRICTIONLESS:
863 2017787 : switch (_formulation)
864 : {
865 : case ContactFormulation::KINEMATIC:
866 : {
867 : RealVectorValue jac_vec;
868 5357460 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
869 : {
870 3972467 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
871 3972467 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
872 3972467 : _var_objects[i]->scalingFactor();
873 : }
874 1384993 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
875 1384993 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
876 1384993 : pinfo->_normal(_component) * pinfo->_normal(_component);
877 : }
878 :
879 632794 : case ContactFormulation::PENALTY:
880 : case ContactFormulation::AUGMENTED_LAGRANGE:
881 632794 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
882 632794 : pinfo->_normal(_component) * pinfo->_normal(_component);
883 :
884 0 : default:
885 0 : mooseError("Invalid contact formulation");
886 : }
887 :
888 277238 : case ContactModel::COULOMB:
889 277238 : switch (_formulation)
890 : {
891 38712 : case ContactFormulation::KINEMATIC:
892 : {
893 38712 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
894 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
895 : {
896 : RealVectorValue jac_vec;
897 114300 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
898 : {
899 76200 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
900 76200 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
901 76200 : _var_objects[i]->scalingFactor();
902 : }
903 38100 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
904 38100 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
905 38100 : pinfo->_normal(_component) * pinfo->_normal(_component);
906 : }
907 : else
908 : {
909 : const Real curr_jac =
910 612 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
911 612 : _connected_dof_indices[_j]) /
912 612 : _var.scalingFactor();
913 612 : return (-curr_jac + _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
914 : }
915 : }
916 :
917 68769 : case ContactFormulation::PENALTY:
918 : {
919 68769 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
920 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
921 14282 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
922 14282 : pinfo->_normal(_component) * pinfo->_normal(_component);
923 : else
924 54487 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp];
925 : }
926 106341 : case ContactFormulation::AUGMENTED_LAGRANGE:
927 : {
928 106341 : Real normal_comp = _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
929 106341 : pinfo->_normal(_component) * pinfo->_normal(_component);
930 :
931 : Real tang_comp = 0.0;
932 106341 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
933 39930 : tang_comp = _phi_secondary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
934 39930 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
935 106341 : return normal_comp + tang_comp;
936 : }
937 :
938 : case ContactFormulation::TANGENTIAL_PENALTY:
939 : {
940 : RealVectorValue jac_vec;
941 190248 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
942 : {
943 126832 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
944 126832 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
945 126832 : _var_objects[i]->scalingFactor();
946 : }
947 63416 : Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
948 63416 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
949 63416 : pinfo->_normal(_component) * pinfo->_normal(_component);
950 :
951 : Real tang_comp = 0.0;
952 63416 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
953 576 : tang_comp = _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
954 576 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
955 :
956 63416 : return normal_comp + tang_comp;
957 : }
958 :
959 0 : default:
960 0 : mooseError("Invalid contact formulation");
961 : }
962 :
963 269880 : case ContactModel::GLUED:
964 269880 : switch (_formulation)
965 : {
966 170883 : case ContactFormulation::KINEMATIC:
967 : {
968 170883 : const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
969 170883 : _connected_dof_indices[_j]) /
970 170883 : _var.scalingFactor();
971 170883 : return (-curr_jac + _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
972 : }
973 :
974 98997 : case ContactFormulation::PENALTY:
975 : case ContactFormulation::AUGMENTED_LAGRANGE:
976 98997 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp];
977 :
978 0 : default:
979 0 : mooseError("Invalid contact formulation");
980 : }
981 0 : default:
982 0 : mooseError("Invalid or unavailable contact model");
983 : }
984 :
985 1829048 : case Moose::SecondaryPrimary:
986 1829048 : switch (_model)
987 : {
988 1408288 : case ContactModel::FRICTIONLESS:
989 1408288 : switch (_formulation)
990 : {
991 949544 : case ContactFormulation::KINEMATIC:
992 : {
993 949544 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
994 :
995 : RealVectorValue jac_vec;
996 3656544 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
997 : {
998 2707000 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
999 5414000 : jac_vec(i) = (*_jacobian)(dof_number,
1000 2707000 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1001 2707000 : _var_objects[i]->scalingFactor();
1002 : }
1003 949544 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1004 949544 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1005 949544 : pinfo->_normal(_component) * pinfo->_normal(_component);
1006 : }
1007 :
1008 458744 : case ContactFormulation::PENALTY:
1009 : case ContactFormulation::AUGMENTED_LAGRANGE:
1010 458744 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1011 458744 : pinfo->_normal(_component) * pinfo->_normal(_component);
1012 :
1013 0 : default:
1014 0 : mooseError("Invalid contact formulation");
1015 : }
1016 :
1017 222424 : case ContactModel::COULOMB:
1018 222424 : switch (_formulation)
1019 : {
1020 29376 : case ContactFormulation::KINEMATIC:
1021 : {
1022 29376 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1023 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1024 : {
1025 28944 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1026 :
1027 : RealVectorValue jac_vec;
1028 86832 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1029 : {
1030 57888 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1031 57888 : jac_vec(i) =
1032 57888 : (*_jacobian)(dof_number,
1033 57888 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1034 57888 : _var_objects[i]->scalingFactor();
1035 : }
1036 28944 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1037 28944 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1038 28944 : pinfo->_normal(_component) * pinfo->_normal(_component);
1039 : }
1040 : else
1041 : {
1042 432 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1043 : const Real curr_jac =
1044 432 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1045 432 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1046 432 : _var.scalingFactor();
1047 432 : return (-curr_jac - _phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
1048 : }
1049 : }
1050 :
1051 55756 : case ContactFormulation::PENALTY:
1052 : {
1053 55756 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1054 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1055 10700 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1056 10700 : pinfo->_normal(_component) * pinfo->_normal(_component);
1057 : else
1058 45056 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp];
1059 : }
1060 87500 : case ContactFormulation::AUGMENTED_LAGRANGE:
1061 : {
1062 87500 : Real normal_comp = -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1063 87500 : pinfo->_normal(_component) * pinfo->_normal(_component);
1064 :
1065 : Real tang_comp = 0.0;
1066 87500 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1067 34720 : tang_comp = -_phi_primary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
1068 34720 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1069 87500 : return normal_comp + tang_comp;
1070 : }
1071 :
1072 49792 : case ContactFormulation::TANGENTIAL_PENALTY:
1073 : {
1074 49792 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1075 :
1076 : RealVectorValue jac_vec;
1077 149376 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1078 : {
1079 99584 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1080 199168 : jac_vec(i) = (*_jacobian)(dof_number,
1081 99584 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1082 99584 : _var_objects[i]->scalingFactor();
1083 : }
1084 49792 : Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1085 49792 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1086 49792 : pinfo->_normal(_component) * pinfo->_normal(_component);
1087 :
1088 : Real tang_comp = 0.0;
1089 49792 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1090 432 : tang_comp = -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1091 432 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1092 :
1093 49792 : return normal_comp + tang_comp;
1094 : }
1095 :
1096 0 : default:
1097 0 : mooseError("Invalid contact formulation");
1098 : }
1099 198336 : case ContactModel::GLUED:
1100 198336 : switch (_formulation)
1101 : {
1102 119580 : case ContactFormulation::KINEMATIC:
1103 : {
1104 119580 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1105 : const Real curr_jac =
1106 119580 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1107 119580 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1108 119580 : _var.scalingFactor();
1109 119580 : return (-curr_jac - _phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
1110 : }
1111 :
1112 78756 : case ContactFormulation::PENALTY:
1113 : case ContactFormulation::AUGMENTED_LAGRANGE:
1114 78756 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp];
1115 :
1116 0 : default:
1117 0 : mooseError("Invalid contact formulation");
1118 : }
1119 :
1120 0 : default:
1121 0 : mooseError("Invalid or unavailable contact model");
1122 : }
1123 :
1124 40565148 : case Moose::PrimarySecondary:
1125 40565148 : switch (_model)
1126 : {
1127 34401220 : case ContactModel::FRICTIONLESS:
1128 34401220 : switch (_formulation)
1129 : {
1130 : case ContactFormulation::KINEMATIC:
1131 : {
1132 : RealVectorValue jac_vec;
1133 105044016 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1134 : {
1135 78541748 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1136 78541748 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1137 78541748 : _var_objects[i]->scalingFactor();
1138 : }
1139 26502268 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1140 26502268 : _test_primary[_i][_qp];
1141 : }
1142 :
1143 7898952 : case ContactFormulation::PENALTY:
1144 : case ContactFormulation::AUGMENTED_LAGRANGE:
1145 7898952 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1146 7898952 : pinfo->_normal(_component) * pinfo->_normal(_component);
1147 :
1148 0 : default:
1149 0 : mooseError("Invalid contact formulation");
1150 : }
1151 3250904 : case ContactModel::COULOMB:
1152 3250904 : switch (_formulation)
1153 : {
1154 154848 : case ContactFormulation::KINEMATIC:
1155 : {
1156 154848 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1157 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1158 : {
1159 : RealVectorValue jac_vec;
1160 457200 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1161 : {
1162 304800 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1163 304800 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1164 304800 : _var_objects[i]->scalingFactor();
1165 : }
1166 152400 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1167 152400 : _test_primary[_i][_qp];
1168 : }
1169 : else
1170 : {
1171 : const Real secondary_jac =
1172 2448 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1173 2448 : _connected_dof_indices[_j]) /
1174 2448 : _var.scalingFactor();
1175 2448 : return secondary_jac * _test_primary[_i][_qp];
1176 : }
1177 : }
1178 :
1179 1037796 : case ContactFormulation::PENALTY:
1180 : {
1181 1037796 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1182 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1183 265224 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1184 265224 : pinfo->_normal(_component) * pinfo->_normal(_component);
1185 : else
1186 772572 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp];
1187 : }
1188 1804596 : case ContactFormulation::AUGMENTED_LAGRANGE:
1189 : {
1190 1804596 : Real normal_comp = -_phi_secondary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1191 1804596 : pinfo->_normal(_component) * pinfo->_normal(_component);
1192 :
1193 : Real tang_comp = 0.0;
1194 1804596 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1195 617328 : tang_comp = -_phi_secondary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1196 617328 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1197 1804596 : return normal_comp + tang_comp;
1198 : }
1199 :
1200 : case ContactFormulation::TANGENTIAL_PENALTY:
1201 : {
1202 : RealVectorValue jac_vec;
1203 760992 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1204 : {
1205 507328 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1206 507328 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1207 507328 : _var_objects[i]->scalingFactor();
1208 : }
1209 : Real normal_comp =
1210 253664 : pinfo->_normal(_component) * (pinfo->_normal * jac_vec) * _test_primary[_i][_qp];
1211 :
1212 : Real tang_comp = 0.0;
1213 253664 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1214 2304 : tang_comp = -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1215 2304 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1216 :
1217 253664 : return normal_comp + tang_comp;
1218 : }
1219 :
1220 0 : default:
1221 0 : mooseError("Invalid contact formulation");
1222 : }
1223 :
1224 2913024 : case ContactModel::GLUED:
1225 2913024 : switch (_formulation)
1226 : {
1227 1684980 : case ContactFormulation::KINEMATIC:
1228 : {
1229 : const Real secondary_jac =
1230 1684980 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1231 1684980 : _connected_dof_indices[_j]) /
1232 1684980 : _var.scalingFactor();
1233 1684980 : return secondary_jac * _test_primary[_i][_qp];
1234 : }
1235 :
1236 1228044 : case ContactFormulation::PENALTY:
1237 : case ContactFormulation::AUGMENTED_LAGRANGE:
1238 1228044 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp];
1239 :
1240 0 : default:
1241 0 : mooseError("Invalid contact formulation");
1242 : }
1243 :
1244 0 : default:
1245 0 : mooseError("Invalid or unavailable contact model");
1246 : }
1247 :
1248 28563712 : case Moose::PrimaryPrimary:
1249 28563712 : switch (_model)
1250 : {
1251 23756768 : case ContactModel::FRICTIONLESS:
1252 23756768 : switch (_formulation)
1253 : {
1254 : case ContactFormulation::KINEMATIC:
1255 : return 0.0;
1256 :
1257 5838752 : case ContactFormulation::PENALTY:
1258 : case ContactFormulation::AUGMENTED_LAGRANGE:
1259 5838752 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1260 5838752 : pinfo->_normal(_component) * pinfo->_normal(_component);
1261 :
1262 0 : default:
1263 0 : mooseError("Invalid contact formulation");
1264 : }
1265 :
1266 4806944 : case ContactModel::COULOMB:
1267 : case ContactModel::GLUED:
1268 4806944 : switch (_formulation)
1269 : {
1270 : case ContactFormulation::KINEMATIC:
1271 : return 0.0;
1272 :
1273 1818880 : case ContactFormulation::PENALTY:
1274 : {
1275 1818880 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1276 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1277 197808 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1278 197808 : pinfo->_normal(_component) * pinfo->_normal(_component);
1279 : else
1280 1621072 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp];
1281 : }
1282 :
1283 199168 : case ContactFormulation::TANGENTIAL_PENALTY:
1284 : {
1285 : Real tang_comp = 0.0;
1286 199168 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1287 1728 : tang_comp = _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1288 1728 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1289 : return tang_comp; // normal component is zero
1290 : }
1291 :
1292 1467760 : case ContactFormulation::AUGMENTED_LAGRANGE:
1293 : {
1294 1467760 : Real normal_comp = _phi_primary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1295 1467760 : pinfo->_normal(_component) * pinfo->_normal(_component);
1296 :
1297 : Real tang_comp = 0.0;
1298 1467760 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1299 533120 : tang_comp = _phi_primary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1300 533120 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1301 1467760 : return normal_comp + tang_comp;
1302 : }
1303 :
1304 0 : default:
1305 0 : mooseError("Invalid contact formulation");
1306 : }
1307 :
1308 0 : default:
1309 0 : mooseError("Invalid or unavailable contact model");
1310 : }
1311 : }
1312 :
1313 : return 0.0;
1314 : }
1315 :
1316 : Real
1317 13095784 : MechanicalContactConstraint::computeQpOffDiagJacobian(Moose::ConstraintJacobianType type,
1318 : unsigned int jvar)
1319 : {
1320 13095784 : PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
1321 :
1322 13095784 : const Real penalty = getPenalty(*_current_node);
1323 13095784 : const Real penalty_slip = getTangentialPenalty(*_current_node);
1324 :
1325 : unsigned int coupled_component;
1326 : Real normal_component_in_coupled_var_dir = 1.0;
1327 13095784 : if (getCoupledVarComponent(jvar, coupled_component))
1328 13095784 : normal_component_in_coupled_var_dir = pinfo->_normal(coupled_component);
1329 :
1330 13095784 : switch (type)
1331 : {
1332 0 : default:
1333 0 : mooseError("Unhandled ConstraintJacobianType");
1334 :
1335 615816 : case Moose::SecondarySecondary:
1336 615816 : switch (_model)
1337 : {
1338 517408 : case ContactModel::FRICTIONLESS:
1339 517408 : switch (_formulation)
1340 : {
1341 : case ContactFormulation::KINEMATIC:
1342 : {
1343 : RealVectorValue jac_vec;
1344 1128600 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1345 : {
1346 839304 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1347 839304 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1348 839304 : _var_objects[i]->scalingFactor();
1349 : }
1350 289296 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1351 289296 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1352 289296 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1353 : }
1354 :
1355 228112 : case ContactFormulation::PENALTY:
1356 : case ContactFormulation::AUGMENTED_LAGRANGE:
1357 228112 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1358 228112 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1359 :
1360 0 : default:
1361 0 : mooseError("Invalid contact formulation");
1362 : }
1363 :
1364 98408 : case ContactModel::COULOMB:
1365 : {
1366 98408 : if ((_formulation == ContactFormulation::KINEMATIC ||
1367 89152 : _formulation == ContactFormulation::TANGENTIAL_PENALTY) &&
1368 89152 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1369 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1370 : {
1371 : RealVectorValue jac_vec;
1372 264864 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1373 : {
1374 176576 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1375 176576 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1376 176576 : _var_objects[i]->scalingFactor();
1377 : }
1378 :
1379 88288 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1380 88288 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1381 88288 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1382 : }
1383 10120 : else if ((_formulation == ContactFormulation::PENALTY) &&
1384 9256 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1385 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1386 792 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1387 792 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1388 9328 : else if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
1389 : {
1390 0 : Real normal_comp = _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1391 0 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1392 :
1393 : Real tang_comp = 0.0;
1394 0 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1395 0 : tang_comp = _phi_secondary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
1396 0 : (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1397 0 : return normal_comp + tang_comp;
1398 : }
1399 : else
1400 : {
1401 9328 : const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1402 9328 : _connected_dof_indices[_j]);
1403 9328 : return -curr_jac;
1404 : }
1405 : }
1406 :
1407 0 : case ContactModel::GLUED:
1408 : {
1409 0 : const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1410 0 : _connected_dof_indices[_j]);
1411 0 : return -curr_jac;
1412 : }
1413 0 : default:
1414 0 : mooseError("Invalid or unavailable contact model");
1415 : }
1416 :
1417 454144 : case Moose::SecondaryPrimary:
1418 454144 : switch (_model)
1419 : {
1420 377840 : case ContactModel::FRICTIONLESS:
1421 377840 : switch (_formulation)
1422 : {
1423 218176 : case ContactFormulation::KINEMATIC:
1424 : {
1425 218176 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1426 :
1427 : RealVectorValue jac_vec;
1428 849648 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1429 : {
1430 631472 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1431 1262944 : jac_vec(i) = (*_jacobian)(dof_number,
1432 631472 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1433 631472 : _var_objects[i]->scalingFactor();
1434 : }
1435 218176 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1436 218176 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1437 218176 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1438 : }
1439 :
1440 159664 : case ContactFormulation::PENALTY:
1441 : case ContactFormulation::AUGMENTED_LAGRANGE:
1442 159664 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1443 159664 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1444 :
1445 0 : default:
1446 0 : mooseError("Invalid contact formulation");
1447 : }
1448 :
1449 76304 : case ContactModel::COULOMB:
1450 76304 : if ((_formulation == ContactFormulation::KINEMATIC ||
1451 70240 : _formulation == ContactFormulation::TANGENTIAL_PENALTY) &&
1452 70240 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1453 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1454 : {
1455 69592 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1456 :
1457 : RealVectorValue jac_vec;
1458 208776 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1459 : {
1460 139184 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1461 139184 : jac_vec(i) =
1462 139184 : (*_jacobian)(dof_number, curr_primary_node->dof_number(0, _vars[_component], 0)) /
1463 139184 : _var_objects[i]->scalingFactor();
1464 : }
1465 :
1466 69592 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1467 69592 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1468 69592 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1469 : }
1470 6712 : else if ((_formulation == ContactFormulation::PENALTY) &&
1471 6064 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1472 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1473 528 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1474 528 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1475 6184 : else if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
1476 : {
1477 0 : Real normal_comp = -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1478 0 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1479 :
1480 : Real tang_comp = 0.0;
1481 0 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1482 0 : tang_comp = -_phi_primary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
1483 0 : (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1484 0 : return normal_comp + tang_comp;
1485 : }
1486 : else
1487 : return 0.0;
1488 :
1489 : case ContactModel::GLUED:
1490 : return 0;
1491 :
1492 0 : default:
1493 0 : mooseError("Invalid or unavailable contact model");
1494 : }
1495 :
1496 6931360 : case Moose::PrimarySecondary:
1497 6931360 : switch (_model)
1498 : {
1499 6537728 : case ContactModel::FRICTIONLESS:
1500 6537728 : switch (_formulation)
1501 : {
1502 : case ContactFormulation::KINEMATIC:
1503 : {
1504 : RealVectorValue jac_vec;
1505 20496000 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1506 : {
1507 15331520 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1508 15331520 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1509 15331520 : _var_objects[i]->scalingFactor();
1510 : }
1511 5164480 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1512 5164480 : _test_primary[_i][_qp];
1513 : }
1514 :
1515 1373248 : case ContactFormulation::PENALTY:
1516 : case ContactFormulation::AUGMENTED_LAGRANGE:
1517 1373248 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1518 1373248 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1519 :
1520 0 : default:
1521 0 : mooseError("Invalid contact formulation");
1522 : }
1523 393632 : case ContactModel::COULOMB:
1524 393632 : switch (_formulation)
1525 : {
1526 128336 : case ContactFormulation::KINEMATIC:
1527 : {
1528 128336 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1529 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1530 : {
1531 : RealVectorValue jac_vec;
1532 380256 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1533 : {
1534 253504 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1535 253504 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1536 253504 : _var_objects[i]->scalingFactor();
1537 : }
1538 126752 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1539 126752 : _test_primary[_i][_qp];
1540 : }
1541 : else
1542 : {
1543 : const Real secondary_jac =
1544 1584 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1545 1584 : _connected_dof_indices[_j]) /
1546 1584 : _var.scalingFactor();
1547 1584 : return secondary_jac * _test_primary[_i][_qp];
1548 : }
1549 : }
1550 :
1551 37024 : case ContactFormulation::PENALTY:
1552 : {
1553 37024 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1554 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1555 3168 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1556 3168 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1557 : else
1558 : return 0.0;
1559 : }
1560 0 : case ContactFormulation::AUGMENTED_LAGRANGE:
1561 : {
1562 0 : Real normal_comp = -_phi_secondary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1563 0 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1564 :
1565 : Real tang_comp = 0.0;
1566 0 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1567 0 : tang_comp = -_phi_secondary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1568 0 : (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1569 0 : return normal_comp + tang_comp;
1570 : }
1571 228272 : case ContactFormulation::TANGENTIAL_PENALTY:
1572 : {
1573 228272 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1574 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1575 : {
1576 : RealVectorValue jac_vec;
1577 679200 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1578 : {
1579 452800 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1580 452800 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1581 452800 : _var_objects[i]->scalingFactor();
1582 : }
1583 226400 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1584 226400 : _test_primary[_i][_qp];
1585 : }
1586 : else
1587 : return 0.0;
1588 : }
1589 :
1590 0 : default:
1591 0 : mooseError("Invalid contact formulation");
1592 : }
1593 :
1594 0 : case ContactModel::GLUED:
1595 0 : switch (_formulation)
1596 : {
1597 0 : case ContactFormulation::KINEMATIC:
1598 : {
1599 : const Real secondary_jac =
1600 0 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1601 0 : _connected_dof_indices[_j]) /
1602 0 : _var.scalingFactor();
1603 0 : return secondary_jac * _test_primary[_i][_qp];
1604 : }
1605 :
1606 : case ContactFormulation::PENALTY:
1607 : case ContactFormulation::AUGMENTED_LAGRANGE:
1608 : return 0.0;
1609 :
1610 0 : default:
1611 0 : mooseError("Invalid contact formulation");
1612 : }
1613 :
1614 0 : default:
1615 0 : mooseError("Invalid or unavailable contact model");
1616 : }
1617 :
1618 5094464 : case Moose::PrimaryPrimary:
1619 5094464 : switch (_model)
1620 : {
1621 4789248 : case ContactModel::FRICTIONLESS:
1622 4789248 : switch (_formulation)
1623 : {
1624 : case ContactFormulation::KINEMATIC:
1625 : return 0.0;
1626 :
1627 915136 : case ContactFormulation::PENALTY:
1628 : case ContactFormulation::AUGMENTED_LAGRANGE:
1629 915136 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1630 915136 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1631 :
1632 0 : default:
1633 0 : mooseError("Invalid contact formulation");
1634 : }
1635 :
1636 305216 : case ContactModel::COULOMB:
1637 : {
1638 305216 : if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
1639 : {
1640 0 : Real normal_comp = _phi_primary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1641 0 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1642 :
1643 0 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1644 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1645 : return normal_comp;
1646 : else
1647 0 : return 0.0;
1648 : }
1649 305216 : else if (_formulation == ContactFormulation::PENALTY &&
1650 24256 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1651 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1652 2112 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1653 2112 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1654 : else
1655 : return 0.0;
1656 : }
1657 :
1658 0 : case ContactModel::GLUED:
1659 0 : if (_formulation == ContactFormulation::PENALTY &&
1660 0 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1661 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1662 0 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1663 0 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1664 0 : else if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
1665 : {
1666 0 : Real normal_comp = _phi_primary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1667 0 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1668 :
1669 : Real tang_comp = 0.0;
1670 0 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1671 0 : tang_comp = _phi_primary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1672 0 : (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1673 0 : return normal_comp + tang_comp;
1674 : }
1675 : else
1676 : return 0.0;
1677 :
1678 0 : default:
1679 0 : mooseError("Invalid or unavailable contact model");
1680 : }
1681 : }
1682 :
1683 : return 0.0;
1684 : }
1685 :
1686 : Real
1687 2526663 : MechanicalContactConstraint::gapOffset(const Node & node)
1688 : {
1689 : Real val = 0;
1690 :
1691 2526663 : if (_has_secondary_gap_offset)
1692 95202 : val += _secondary_gap_offset_var->getNodalValue(node);
1693 :
1694 2526663 : if (_has_mapped_primary_gap_offset)
1695 95202 : val += _mapped_primary_gap_offset_var->getNodalValue(node);
1696 :
1697 2526663 : return val;
1698 : }
1699 :
1700 : Real
1701 38608901 : MechanicalContactConstraint::nodalArea(const Node & node)
1702 : {
1703 38608901 : dof_id_type dof = node.dof_number(_aux_system.number(), _nodal_area_var->number(), 0);
1704 :
1705 38608901 : Real area = (*_aux_solution)(dof);
1706 38608901 : if (area == 0.0)
1707 : {
1708 0 : if (_t_step > 1)
1709 0 : mooseError("Zero nodal area found");
1710 : else
1711 : area = 1.0; // Avoid divide by zero during initialization
1712 : }
1713 :
1714 38608901 : return area;
1715 : }
1716 :
1717 : Real
1718 88717044 : MechanicalContactConstraint::getPenalty(const Node & node)
1719 : {
1720 88717044 : Real penalty = _penalty;
1721 88717044 : if (_normalize_penalty)
1722 19308376 : penalty *= nodalArea(node);
1723 88717044 : return penalty * MathUtils::pow(_penalty_multiplier, _lagrangian_iteration_number);
1724 : }
1725 :
1726 : Real
1727 87521134 : MechanicalContactConstraint::getTangentialPenalty(const Node & node)
1728 : {
1729 87521134 : Real penalty = _penalty_tangential;
1730 87521134 : if (_normalize_penalty)
1731 19223405 : penalty *= nodalArea(node);
1732 :
1733 87521134 : return penalty * MathUtils::pow(_penalty_multiplier, _lagrangian_iteration_number);
1734 : }
1735 :
1736 : void
1737 188647 : MechanicalContactConstraint::computeJacobian()
1738 : {
1739 188647 : getConnectedDofIndices(_var.number());
1740 :
1741 188647 : prepareMatrixTagNeighbor(
1742 188647 : _assembly, _primary_var.number(), _var.number(), Moose::NeighborNeighbor, _Knn);
1743 :
1744 188647 : _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
1745 :
1746 377294 : for (_i = 0; _i < _test_secondary.size(); _i++)
1747 : // Loop over the connected dof indices so we can get all the jacobian contributions
1748 2753552 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1749 2564905 : _Kee(_i, _j) += computeQpJacobian(Moose::SecondarySecondary);
1750 :
1751 188647 : if (_primary_secondary_jacobian)
1752 : {
1753 185617 : prepareMatrixTagNeighbor(_assembly, _var.number(), _var.number(), Moose::ElementNeighbor, _Ken);
1754 185617 : if (_Ken.m() && _Ken.n())
1755 : {
1756 371234 : for (_i = 0; _i < _test_secondary.size(); _i++)
1757 2014665 : for (_j = 0; _j < _phi_primary.size(); _j++)
1758 1829048 : _Ken(_i, _j) += computeQpJacobian(Moose::SecondaryPrimary);
1759 185617 : accumulateTaggedLocalMatrix(
1760 185617 : _assembly, _var.number(), _var.number(), Moose::ElementNeighbor, _Ken);
1761 : }
1762 :
1763 185617 : _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
1764 2014665 : for (_i = 0; _i < _test_primary.size(); _i++)
1765 : // Loop over the connected dof indices so we can get all the jacobian contributions
1766 42394196 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1767 40565148 : _Kne(_i, _j) += computeQpJacobian(Moose::PrimarySecondary);
1768 : }
1769 :
1770 188647 : if (_Knn.m() && _Knn.n())
1771 : {
1772 2041935 : for (_i = 0; _i < _test_primary.size(); _i++)
1773 30417000 : for (_j = 0; _j < _phi_primary.size(); _j++)
1774 28563712 : _Knn(_i, _j) += computeQpJacobian(Moose::PrimaryPrimary);
1775 188647 : accumulateTaggedLocalMatrix(
1776 188647 : _assembly, _primary_var.number(), _var.number(), Moose::NeighborNeighbor, _Knn);
1777 : }
1778 188647 : }
1779 :
1780 : void
1781 65750 : MechanicalContactConstraint::computeOffDiagJacobian(const unsigned int jvar_num)
1782 : {
1783 65750 : getConnectedDofIndices(jvar_num);
1784 :
1785 65750 : _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
1786 :
1787 65750 : prepareMatrixTagNeighbor(
1788 65750 : _assembly, _primary_var.number(), jvar_num, Moose::NeighborNeighbor, _Knn);
1789 :
1790 131500 : for (_i = 0; _i < _test_secondary.size(); _i++)
1791 : // Loop over the connected dof indices so we can get all the jacobian contributions
1792 681566 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1793 615816 : _Kee(_i, _j) += computeQpOffDiagJacobian(Moose::SecondarySecondary, jvar_num);
1794 :
1795 65750 : if (_primary_secondary_jacobian)
1796 : {
1797 65750 : prepareMatrixTagNeighbor(_assembly, _var.number(), jvar_num, Moose::ElementNeighbor, _Ken);
1798 131500 : for (_i = 0; _i < _test_secondary.size(); _i++)
1799 519894 : for (_j = 0; _j < _phi_primary.size(); _j++)
1800 454144 : _Ken(_i, _j) += computeQpOffDiagJacobian(Moose::SecondaryPrimary, jvar_num);
1801 65750 : accumulateTaggedLocalMatrix(_assembly, _var.number(), jvar_num, Moose::ElementNeighbor, _Ken);
1802 :
1803 65750 : _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
1804 65750 : if (_Kne.m() && _Kne.n())
1805 519894 : for (_i = 0; _i < _test_primary.size(); _i++)
1806 : // Loop over the connected dof indices so we can get all the jacobian contributions
1807 7385504 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1808 6931360 : _Kne(_i, _j) += computeQpOffDiagJacobian(Moose::PrimarySecondary, jvar_num);
1809 : }
1810 :
1811 519894 : for (_i = 0; _i < _test_primary.size(); _i++)
1812 5548608 : for (_j = 0; _j < _phi_primary.size(); _j++)
1813 5094464 : _Knn(_i, _j) += computeQpOffDiagJacobian(Moose::PrimaryPrimary, jvar_num);
1814 65750 : accumulateTaggedLocalMatrix(
1815 65750 : _assembly, _primary_var.number(), jvar_num, Moose::NeighborNeighbor, _Knn);
1816 65750 : }
1817 :
1818 : void
1819 254397 : MechanicalContactConstraint::getConnectedDofIndices(unsigned int var_num)
1820 : {
1821 : unsigned int component;
1822 254397 : if (getCoupledVarComponent(var_num, component) || _non_displacement_vars_jacobian)
1823 : {
1824 254397 : if (_primary_secondary_jacobian && _connected_secondary_nodes_jacobian)
1825 251367 : NodeFaceConstraint::getConnectedDofIndices(var_num);
1826 : else
1827 : {
1828 3030 : _connected_dof_indices.clear();
1829 3030 : MooseVariableFEBase & var = _sys.getVariable(0, var_num);
1830 3030 : _connected_dof_indices.push_back(var.nodalDofIndex());
1831 : }
1832 : }
1833 :
1834 254397 : _phi_secondary.resize(_connected_dof_indices.size());
1835 :
1836 254397 : dof_id_type current_node_var_dof_index = _sys.getVariable(0, var_num).nodalDofIndex();
1837 254397 : _qp = 0;
1838 :
1839 : // Fill up _phi_secondary so that it is 1 when j corresponds to the dof associated with this node
1840 : // and 0 for every other dof
1841 : // This corresponds to evaluating all of the connected shape functions at _this_ node
1842 3435118 : for (unsigned int j = 0; j < _connected_dof_indices.size(); j++)
1843 : {
1844 3180721 : _phi_secondary[j].resize(1);
1845 :
1846 3180721 : if (_connected_dof_indices[j] == current_node_var_dof_index)
1847 254397 : _phi_secondary[j][_qp] = 1.0;
1848 : else
1849 2926324 : _phi_secondary[j][_qp] = 0.0;
1850 : }
1851 254397 : }
1852 :
1853 : bool
1854 13350181 : MechanicalContactConstraint::getCoupledVarComponent(unsigned int var_num, unsigned int & component)
1855 : {
1856 13350181 : component = std::numeric_limits<unsigned int>::max();
1857 : bool coupled_var_is_disp_var = false;
1858 25509339 : for (const auto i : make_range(Moose::dim))
1859 : {
1860 25509339 : if (var_num == _vars[i])
1861 : {
1862 : coupled_var_is_disp_var = true;
1863 13350181 : component = i;
1864 13350181 : break;
1865 : }
1866 : }
1867 :
1868 13350181 : return coupled_var_is_disp_var;
1869 : }
1870 :
1871 : void
1872 79814 : MechanicalContactConstraint::residualEnd()
1873 : {
1874 79814 : if (_component == 0 && (_print_contact_nodes || _contact_linesearch))
1875 : {
1876 86 : _communicator.set_union(_current_contact_state);
1877 86 : if (_print_contact_nodes)
1878 : {
1879 0 : if (_current_contact_state == _old_contact_state)
1880 0 : _console << "Unchanged contact state. " << _current_contact_state.size()
1881 0 : << " nodes in contact.\n";
1882 : else
1883 0 : _console << "Changed contact state!!! " << _current_contact_state.size()
1884 0 : << " nodes in contact.\n";
1885 : }
1886 86 : if (_contact_linesearch)
1887 86 : _contact_linesearch->insertSet(_current_contact_state);
1888 : _old_contact_state.swap(_current_contact_state);
1889 : _current_contact_state.clear();
1890 : }
1891 79814 : }
|
