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 6433 : MechanicalContactConstraint::validParams()
36 : {
37 6433 : InputParameters params = NodeFaceConstraint::validParams();
38 6433 : params += ContactAction::commonParameters();
39 :
40 12866 : params.addRequiredParam<BoundaryName>("boundary", "The primary boundary");
41 12866 : params.addParam<BoundaryName>("secondary", "The secondary boundary");
42 12866 : 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 12866 : params.addCoupledVar(
48 : "displacements",
49 : "The displacements appropriate for the simulation geometry and coordinate system");
50 :
51 12866 : params.addCoupledVar("secondary_gap_offset", "offset to the gap distance from secondary side");
52 12866 : params.addCoupledVar("mapped_primary_gap_offset",
53 : "offset to the gap distance mapped from primary side");
54 12866 : params.addRequiredCoupledVar("nodal_area", "The nodal area");
55 :
56 6433 : params.set<bool>("use_displaced_mesh") = true;
57 12866 : params.addParam<Real>(
58 : "penalty",
59 12866 : 1e8,
60 : "The penalty to apply. This can vary depending on the stiffness of your materials");
61 12866 : params.addParam<Real>("penalty_multiplier",
62 12866 : 1.0,
63 : "The growth factor for the penalty applied at the end of each augmented "
64 : "Lagrange update iteration");
65 12866 : params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
66 12866 : params.addParam<Real>("tangential_tolerance",
67 : "Tangential distance to extend edges of contact surfaces");
68 12866 : params.addParam<Real>(
69 12866 : "capture_tolerance", 0, "Normal distance from surface within which nodes are captured");
70 :
71 12866 : params.addParam<Real>("tension_release",
72 12866 : 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 12866 : params.addParam<bool>(
78 : "normalize_penalty",
79 12866 : false,
80 : "Whether to normalize the penalty parameter with the nodal area for penalty contact.");
81 12866 : params.addParam<bool>(
82 : "primary_secondary_jacobian",
83 12866 : true,
84 : "Whether to include jacobian entries coupling primary and secondary nodes.");
85 12866 : params.addParam<bool>(
86 : "connected_secondary_nodes_jacobian",
87 12866 : true,
88 : "Whether to include jacobian entries coupling nodes connected to secondary nodes.");
89 12866 : params.addParam<bool>("non_displacement_variables_jacobian",
90 12866 : true,
91 : "Whether to include jacobian entries coupling with variables that are not "
92 : "displacement variables.");
93 12866 : params.addParam<unsigned int>("stick_lock_iterations",
94 12866 : 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 12866 : params.addParam<Real>("stick_unlock_factor",
98 12866 : 1.5,
99 : "Factor by which frictional capacity must be "
100 : "exceeded to permit stick-locked node to slip "
101 : "again.");
102 12866 : params.addParam<Real>("al_penetration_tolerance",
103 : "The tolerance of the penetration for augmented Lagrangian method.");
104 12866 : params.addParam<Real>("al_incremental_slip_tolerance",
105 : "The tolerance of the incremental slip for augmented Lagrangian method.");
106 :
107 12866 : params.addParam<Real>("al_frictional_force_tolerance",
108 : "The tolerance of the frictional force for augmented Lagrangian method.");
109 12866 : params.addParam<bool>(
110 12866 : "print_contact_nodes", false, "Whether to print the number of nodes in contact.");
111 :
112 6433 : 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 6433 : return params;
119 0 : }
120 :
121 : Threads::spin_mutex MechanicalContactConstraint::_contact_set_mutex;
122 :
123 4562 : MechanicalContactConstraint::MechanicalContactConstraint(const InputParameters & parameters)
124 : : NodeFaceConstraint(parameters),
125 4562 : _displaced_problem(parameters.get<FEProblemBase *>("_fe_problem_base")->getDisplacedProblem()),
126 9124 : _component(getParam<unsigned int>("component")),
127 9124 : _model(getParam<MooseEnum>("model").getEnum<ContactModel>()),
128 9124 : _formulation(getParam<MooseEnum>("formulation").getEnum<ContactFormulation>()),
129 9124 : _normalize_penalty(getParam<bool>("normalize_penalty")),
130 9124 : _penalty(getParam<Real>("penalty")),
131 9124 : _penalty_multiplier(getParam<Real>("penalty_multiplier")),
132 9124 : _friction_coefficient(getParam<Real>("friction_coefficient")),
133 9124 : _tension_release(getParam<Real>("tension_release")),
134 9124 : _capture_tolerance(getParam<Real>("capture_tolerance")),
135 9124 : _stick_lock_iterations(getParam<unsigned int>("stick_lock_iterations")),
136 9124 : _stick_unlock_factor(getParam<Real>("stick_unlock_factor")),
137 4562 : _update_stateful_data(true),
138 4562 : _residual_copy(_sys.residualGhosted()),
139 4562 : _mesh_dimension(_mesh.dimension()),
140 4562 : _vars(3, libMesh::invalid_uint),
141 4562 : _var_objects(3, nullptr),
142 4562 : _has_secondary_gap_offset(isCoupled("secondary_gap_offset")),
143 4562 : _secondary_gap_offset_var(_has_secondary_gap_offset ? getVar("secondary_gap_offset", 0)
144 : : nullptr),
145 4562 : _has_mapped_primary_gap_offset(isCoupled("mapped_primary_gap_offset")),
146 4562 : _mapped_primary_gap_offset_var(
147 4562 : _has_mapped_primary_gap_offset ? getVar("mapped_primary_gap_offset", 0) : nullptr),
148 4562 : _nodal_area_var(getVar("nodal_area", 0)),
149 4562 : _aux_system(_nodal_area_var->sys()),
150 4562 : _aux_solution(_aux_system.currentSolution()),
151 9124 : _primary_secondary_jacobian(getParam<bool>("primary_secondary_jacobian")),
152 9124 : _connected_secondary_nodes_jacobian(getParam<bool>("connected_secondary_nodes_jacobian")),
153 9124 : _non_displacement_vars_jacobian(getParam<bool>("non_displacement_variables_jacobian")),
154 4562 : _contact_linesearch(dynamic_cast<ContactLineSearchBase *>(_subproblem.getLineSearch())),
155 9124 : _print_contact_nodes(getParam<bool>("print_contact_nodes")),
156 4562 : _augmented_lagrange_problem(
157 4562 : dynamic_cast<AugmentedLagrangianContactProblemInterface *>(&_fe_problem)),
158 4562 : _lagrangian_iteration_number(_augmented_lagrange_problem
159 4562 : ? _augmented_lagrange_problem->getLagrangianIterationNumber()
160 9124 : : _no_iterations)
161 : {
162 4562 : _overwrite_secondary_residual = false;
163 :
164 9124 : if (isParamValid("displacements"))
165 : {
166 : // modern parameter scheme for displacements
167 30008 : for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
168 : {
169 10442 : _vars[i] = coupled("displacements", i);
170 10442 : _var_objects[i] = getVar("displacements", i);
171 : }
172 : }
173 :
174 9124 : if (parameters.isParamValid("tangential_tolerance"))
175 9189 : _penetration_locator.setTangentialTolerance(getParam<Real>("tangential_tolerance"));
176 :
177 9124 : if (parameters.isParamValid("normal_smoothing_distance"))
178 1350 : _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
179 :
180 9124 : if (parameters.isParamValid("normal_smoothing_method"))
181 0 : _penetration_locator.setNormalSmoothingMethod(
182 : parameters.get<std::string>("normal_smoothing_method"));
183 :
184 4562 : 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 4562 : if (_model == ContactModel::GLUED)
188 567 : _penetration_locator.setUpdate(false);
189 :
190 4562 : 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 4562 : _penalty_tangential = _penalty;
195 :
196 4562 : if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
197 : {
198 495 : if (_model == ContactModel::GLUED)
199 0 : mooseError("The Augmented Lagrangian contact formulation does not support GLUED case.");
200 :
201 495 : if (!_augmented_lagrange_problem)
202 0 : mooseError("The Augmented Lagrangian contact formulation must use "
203 : "AugmentedLagrangianContactProblem.");
204 :
205 990 : if (!parameters.isParamValid("al_penetration_tolerance"))
206 0 : mooseError("For Augmented Lagrangian contact, al_penetration_tolerance must be provided.");
207 : else
208 495 : _al_penetration_tolerance = parameters.get<Real>("al_penetration_tolerance");
209 :
210 495 : if (_model != ContactModel::FRICTIONLESS)
211 : {
212 630 : if (!parameters.isParamValid("al_incremental_slip_tolerance") ||
213 630 : !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 210 : _al_incremental_slip_tolerance = parameters.get<Real>("al_incremental_slip_tolerance");
222 210 : _al_frictional_force_tolerance = parameters.get<Real>("al_frictional_force_tolerance");
223 : }
224 : }
225 : }
226 4562 : }
227 :
228 : void
229 18040 : MechanicalContactConstraint::timestepSetup()
230 : {
231 18040 : if (_component == 0)
232 : {
233 8661 : updateContactStatefulData(/* beginning_of_step = */ true);
234 8661 : if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
235 153 : updateAugmentedLagrangianMultiplier(/* beginning_of_step = */ true);
236 :
237 8661 : _update_stateful_data = false;
238 :
239 8661 : if (_contact_linesearch)
240 14 : _contact_linesearch->reset();
241 : }
242 18040 : }
243 :
244 : void
245 52667 : MechanicalContactConstraint::jacobianSetup()
246 : {
247 52667 : if (_component == 0)
248 : {
249 25538 : if (_update_stateful_data)
250 17483 : updateContactStatefulData(/* beginning_of_step = */ false);
251 25538 : _update_stateful_data = true;
252 : }
253 52667 : }
254 :
255 : void
256 309 : MechanicalContactConstraint::updateAugmentedLagrangianMultiplier(bool beginning_of_step)
257 : {
258 2380 : for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
259 : {
260 2071 : if (!pinfo)
261 0 : continue;
262 :
263 2071 : const Node & node = _mesh.nodeRef(secondary_node_num);
264 2071 : if (node.n_comp(_sys.number(), _vars[_component]) < 1)
265 0 : continue;
266 :
267 2071 : const Real distance = pinfo->_normal * (pinfo->_closest_point - node) - gapOffset(node);
268 :
269 2071 : if (beginning_of_step && _model == ContactModel::COULOMB)
270 : {
271 350 : pinfo->_lagrange_multiplier_slip.zero();
272 350 : if (pinfo->isCaptured())
273 0 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
274 : }
275 :
276 2071 : if (pinfo->isCaptured())
277 : {
278 1124 : if (_model == ContactModel::FRICTIONLESS)
279 774 : pinfo->_lagrange_multiplier -= getPenalty(node) * distance;
280 :
281 1124 : if (_model == ContactModel::COULOMB)
282 : {
283 350 : if (!beginning_of_step)
284 : {
285 350 : Real penalty = getPenalty(node);
286 : RealVectorValue pen_force_normal =
287 350 : penalty * (-distance) * pinfo->_normal + pinfo->_lagrange_multiplier * pinfo->_normal;
288 :
289 : // update normal lagrangian multiplier
290 350 : pinfo->_lagrange_multiplier += penalty * (-distance);
291 :
292 : // Frictional capacity
293 635 : const Real capacity(_friction_coefficient * (pen_force_normal * pinfo->_normal < 0
294 350 : ? -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 350 : Real penalty_slip = getTangentialPenalty(node);
302 :
303 : RealVectorValue inc_pen_force_tangential =
304 350 : 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 350 : const Real tan_mag(contact_force_tangential.norm());
311 :
312 350 : if (tan_mag > capacity * (_al_frictional_force_tolerance + 1.0))
313 : {
314 45 : pinfo->_lagrange_multiplier_slip =
315 : -tau_old + capacity * contact_force_tangential / tan_mag;
316 45 : if (MooseUtils::absoluteFuzzyEqual(capacity, 0.0))
317 45 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
318 : else
319 0 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
320 : }
321 : else
322 : {
323 305 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
324 : pinfo->_lagrange_multiplier_slip += penalty_slip * tangential_inc_slip;
325 : }
326 : }
327 : }
328 : }
329 : }
330 309 : }
331 :
332 : bool
333 510 : MechanicalContactConstraint::AugmentedLagrangianContactConverged()
334 : {
335 : Real contactResidual = 0.0;
336 510 : unsigned int converged = 0;
337 :
338 3016 : for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
339 : {
340 2662 : if (!pinfo)
341 0 : continue;
342 :
343 2662 : const auto & node = _mesh.nodeRef(secondary_node_num);
344 :
345 : // Skip this pinfo if there are no DOFs on this node.
346 2662 : if (node.n_comp(_sys.number(), _vars[_component]) < 1)
347 0 : continue;
348 :
349 2662 : const Real distance = pinfo->_normal * (pinfo->_closest_point - node) - gapOffset(node);
350 :
351 2662 : if (pinfo->isCaptured())
352 : {
353 2662 : if (contactResidual < std::abs(distance))
354 : contactResidual = std::abs(distance);
355 :
356 : // penetration < tol
357 2662 : if (contactResidual > _al_penetration_tolerance)
358 : {
359 156 : converged = 1;
360 156 : break;
361 : }
362 :
363 2506 : 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 910 : const Real tan_mag(contact_force_tangential.norm());
373 910 : const Real tangential_inc_slip_mag = tangential_inc_slip.norm();
374 :
375 : RealVectorValue distance_vec =
376 910 : (pinfo->_normal * (node - pinfo->_closest_point) + gapOffset(node)) * pinfo->_normal;
377 :
378 910 : Real penalty = getPenalty(node);
379 : RealVectorValue pen_force_normal =
380 910 : penalty * distance_vec + pinfo->_lagrange_multiplier * pinfo->_normal;
381 :
382 : // Frictional capacity
383 1625 : Real capacity(_friction_coefficient * (pen_force_normal * pinfo->_normal < 0
384 910 : ? -pen_force_normal * pinfo->_normal
385 : : 0.0));
386 :
387 : // incremental slip <= tol for all pinfo_pair such that tan_mag < capacity
388 910 : if (MooseUtils::absoluteFuzzyLessThan(tan_mag, capacity) &&
389 715 : pinfo->_mech_status == PenetrationInfo::MS_STICKING)
390 : {
391 715 : 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 910 : if (tan_mag >
401 910 : (1 + _al_frictional_force_tolerance) * (capacity + _al_frictional_force_tolerance))
402 : {
403 0 : converged = 3;
404 0 : break;
405 : }
406 : }
407 : }
408 : }
409 :
410 510 : _communicator.max(converged);
411 :
412 510 : if (converged == 1)
413 156 : _console << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied "
414 156 : << std::endl;
415 354 : else if (converged == 2)
416 0 : _console << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied "
417 0 : << std::endl;
418 354 : 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 26144 : MechanicalContactConstraint::updateContactStatefulData(bool beginning_of_step)
429 : {
430 222794 : for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
431 : {
432 196650 : if (!pinfo)
433 7355 : continue;
434 :
435 189295 : const Node & node = _mesh.nodeRef(secondary_node_num);
436 189295 : if (node.n_comp(_sys.number(), _vars[_component]) < 1)
437 0 : continue;
438 :
439 189295 : if (beginning_of_step)
440 : {
441 59880 : if (_app.getExecutioner()->lastSolveConverged())
442 : {
443 59880 : pinfo->_contact_force_old = pinfo->_contact_force;
444 59880 : pinfo->_accumulated_slip_old = pinfo->_accumulated_slip;
445 59880 : pinfo->_frictional_energy_old = pinfo->_frictional_energy;
446 59880 : pinfo->_mech_status_old = pinfo->_mech_status;
447 : }
448 0 : else if (pinfo->_mech_status_old == PenetrationInfo::MS_NO_CONTACT &&
449 0 : pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
450 : {
451 : // The penetration info object could be based on a bad state so delete it
452 0 : delete pinfo;
453 0 : pinfo = NULL;
454 0 : continue;
455 : }
456 :
457 59880 : pinfo->_locked_this_step = 0;
458 59880 : pinfo->_stick_locked_this_step = 0;
459 59880 : pinfo->_starting_elem = pinfo->_elem;
460 59880 : pinfo->_starting_side_num = pinfo->_side_num;
461 59880 : pinfo->_starting_closest_point_ref = pinfo->_closest_point_ref;
462 : }
463 189295 : pinfo->_incremental_slip_prev_iter = pinfo->_incremental_slip;
464 : }
465 26144 : }
466 :
467 : bool
468 2774272 : MechanicalContactConstraint::shouldApply()
469 : {
470 : bool in_contact = false;
471 :
472 : std::map<dof_id_type, PenetrationInfo *>::iterator found =
473 2774272 : _penetration_locator._penetration_info.find(_current_node->id());
474 2774272 : if (found != _penetration_locator._penetration_info.end())
475 : {
476 2774272 : PenetrationInfo * pinfo = found->second;
477 2774272 : if (pinfo != NULL)
478 : {
479 2774272 : 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 2774272 : if (_component == 0)
484 1189014 : computeContactForce(*_current_node, pinfo, is_nonlinear);
485 :
486 2774272 : if (pinfo->isCaptured())
487 : {
488 : in_contact = true;
489 1805776 : if (is_nonlinear)
490 : {
491 : Threads::spin_mutex::scoped_lock lock(_contact_set_mutex);
492 209148 : _current_contact_state.insert(_current_node->id());
493 : }
494 : }
495 : }
496 : }
497 :
498 2774272 : return in_contact;
499 : }
500 :
501 : void
502 1189014 : MechanicalContactConstraint::computeContactForce(const Node & node,
503 : PenetrationInfo * pinfo,
504 : bool update_contact_set)
505 : {
506 : // Build up residual vector
507 : RealVectorValue res_vec;
508 3963286 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
509 : {
510 2774272 : dof_id_type dof_number = node.dof_number(0, _vars[i], 0);
511 2774272 : res_vec(i) = _residual_copy(dof_number) / _var_objects[i]->scalingFactor();
512 : }
513 :
514 : RealVectorValue distance_vec(node - pinfo->_closest_point);
515 1189014 : if (distance_vec.norm() != 0)
516 1144710 : 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 1189014 : 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 11327 : if (_formulation == ContactFormulation::KINEMATIC ||
534 : _formulation == ContactFormulation::TANGENTIAL_PENALTY)
535 6236 : ++pinfo->_locked_this_step;
536 : }
537 :
538 1189014 : if (!pinfo->isCaptured())
539 431630 : return;
540 :
541 757384 : const Real penalty = getPenalty(node);
542 757384 : const Real penalty_slip = getTangentialPenalty(node);
543 :
544 : RealVectorValue pen_force(penalty * distance_vec);
545 :
546 757384 : switch (_model)
547 : {
548 526224 : case ContactModel::FRICTIONLESS:
549 526224 : switch (_formulation)
550 : {
551 293608 : case ContactFormulation::KINEMATIC:
552 293608 : pinfo->_contact_force = -pinfo->_normal * (pinfo->_normal * res_vec);
553 293608 : break;
554 :
555 212033 : case ContactFormulation::PENALTY:
556 212033 : pinfo->_contact_force = pinfo->_normal * (pinfo->_normal * pen_force);
557 212033 : break;
558 :
559 20583 : case ContactFormulation::AUGMENTED_LAGRANGE:
560 20583 : pinfo->_contact_force =
561 : (pinfo->_normal *
562 : (pinfo->_normal * (pen_force + pinfo->_lagrange_multiplier * pinfo->_normal)));
563 20583 : break;
564 :
565 0 : default:
566 0 : mooseError("Invalid contact formulation");
567 : break;
568 : }
569 526224 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
570 526224 : break;
571 88612 : case ContactModel::COULOMB:
572 88612 : switch (_formulation)
573 : {
574 21571 : case ContactFormulation::KINEMATIC:
575 : {
576 : // Frictional capacity
577 40667 : const Real capacity(_friction_coefficient *
578 21571 : (res_vec * pinfo->_normal > 0 ? res_vec * pinfo->_normal : 0));
579 :
580 : // Normal and tangential components of predictor force
581 21571 : 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 21571 : const Real tan_mag(contact_force_tangential.norm());
592 21571 : const Real tangential_inc_slip_mag = tangential_inc_slip.norm();
593 21571 : const Real slip_tol = capacity / penalty;
594 21571 : pinfo->_slip_tol = slip_tol;
595 :
596 21571 : if ((tangential_inc_slip_mag > slip_tol || tan_mag > capacity) &&
597 20036 : (pinfo->_stick_locked_this_step < _stick_lock_iterations ||
598 0 : tan_mag > capacity * _stick_unlock_factor))
599 : {
600 20036 : 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 20036 : 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 19370 : if (slip_dot_tang_force < capacity)
612 : slipped_too_far = true;
613 : }
614 :
615 : if (slipped_too_far) // slip back along slip increment
616 11165 : pinfo->_contact_force = contact_force_normal + capacity * slip_inc_direction;
617 : else
618 : {
619 8871 : if (tan_mag > 0) // slip along tangential force direction
620 8848 : 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 20036 : if (capacity == 0)
626 2471 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
627 : else
628 17565 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
629 : }
630 : else
631 : {
632 1535 : if (pinfo->_mech_status != PenetrationInfo::MS_STICKING &&
633 : pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
634 387 : ++pinfo->_stick_locked_this_step;
635 1535 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
636 : }
637 : break;
638 : }
639 :
640 13670 : case ContactFormulation::PENALTY:
641 : {
642 13670 : distance_vec =
643 : pinfo->_incremental_slip +
644 13670 : (pinfo->_normal * (node - pinfo->_closest_point) + gapOffset(node)) * pinfo->_normal;
645 : pen_force = penalty * distance_vec;
646 :
647 : // Frictional capacity
648 23311 : const Real capacity(_friction_coefficient *
649 13670 : (pen_force * pinfo->_normal < 0 ? -pen_force * pinfo->_normal : 0));
650 :
651 : // Elastic predictor
652 13670 : 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 13670 : const Real tan_mag(contact_force_tangential.norm());
661 :
662 13670 : if (tan_mag > capacity)
663 : {
664 2384 : pinfo->_contact_force =
665 : contact_force_normal + capacity * contact_force_tangential / tan_mag;
666 2384 : if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
667 1784 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
668 : else
669 600 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
670 : }
671 : else
672 11286 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
673 : break;
674 : }
675 :
676 12540 : case ContactFormulation::AUGMENTED_LAGRANGE:
677 : {
678 12540 : distance_vec =
679 12540 : (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 12540 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
696 5295 : pinfo->_contact_force =
697 : contact_force_normal + contact_force_tangential + inc_pen_force_tangential;
698 : else
699 7245 : pinfo->_contact_force = contact_force_normal + contact_force_tangential;
700 :
701 : break;
702 : }
703 :
704 40831 : case ContactFormulation::TANGENTIAL_PENALTY:
705 : {
706 : // Frictional capacity (kinematic formulation)
707 44339 : 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 40831 : const Real tan_mag(contact_force_tangential.norm());
721 :
722 40831 : if (tan_mag > capacity)
723 : {
724 32514 : pinfo->_contact_force =
725 : contact_force_normal + capacity * contact_force_tangential / tan_mag;
726 32514 : if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
727 4668 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
728 : else
729 27846 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
730 : }
731 : else
732 : {
733 8317 : pinfo->_contact_force = contact_force_normal + contact_force_tangential;
734 8317 : 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 142548 : case ContactModel::GLUED:
746 142548 : switch (_formulation)
747 : {
748 : case ContactFormulation::KINEMATIC:
749 100540 : pinfo->_contact_force = -res_vec;
750 100540 : break;
751 :
752 42008 : case ContactFormulation::PENALTY:
753 42008 : pinfo->_contact_force = pen_force;
754 42008 : 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 142548 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
766 142548 : break;
767 :
768 0 : default:
769 0 : mooseError("Invalid or unavailable contact model");
770 : break;
771 : }
772 :
773 : // Release
774 90842 : if (update_contact_set && _model != ContactModel::GLUED && pinfo->isCaptured() &&
775 824617 : !newly_captured && _tension_release >= 0.0 &&
776 64681 : (_contact_linesearch ? true : pinfo->_locked_this_step < 2))
777 : {
778 64625 : const Real contact_pressure = -(pinfo->_normal * pinfo->_contact_force) / nodalArea(node);
779 64625 : if (-contact_pressure >= _tension_release)
780 : {
781 : pinfo->release();
782 : pinfo->_contact_force.zero();
783 : }
784 : }
785 : }
786 :
787 : Real
788 29023 : MechanicalContactConstraint::computeQpSecondaryValue()
789 : {
790 29023 : return _u_secondary[_qp];
791 : }
792 :
793 : Real
794 16771335 : MechanicalContactConstraint::computeQpResidual(Moose::ConstraintType type)
795 : {
796 16771335 : PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
797 16771335 : Real resid = pinfo->_contact_force(_component);
798 16771335 : switch (type)
799 : {
800 1621090 : case Moose::Secondary:
801 1621090 : if (_formulation == ContactFormulation::KINEMATIC)
802 : {
803 937072 : RealVectorValue distance_vec(*_current_node - pinfo->_closest_point);
804 937072 : if (distance_vec.norm() != 0)
805 891622 : distance_vec += gapOffset(*_current_node) * pinfo->_normal * distance_vec.unit() *
806 1783244 : distance_vec.unit();
807 :
808 937072 : const Real penalty = getPenalty(*_current_node);
809 : RealVectorValue pen_force(penalty * distance_vec);
810 :
811 937072 : if (_model == ContactModel::FRICTIONLESS)
812 653331 : resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
813 283741 : else if (_model == ContactModel::COULOMB)
814 : {
815 : distance_vec = distance_vec - pinfo->_incremental_slip;
816 35916 : pen_force = penalty * distance_vec;
817 35916 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
818 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
819 32946 : resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
820 : else
821 2970 : resid += pen_force(_component);
822 : }
823 247825 : else if (_model == ContactModel::GLUED)
824 247825 : resid += pen_force(_component);
825 : }
826 684018 : else if (_formulation == ContactFormulation::TANGENTIAL_PENALTY &&
827 69832 : _model == ContactModel::COULOMB)
828 : {
829 69832 : RealVectorValue distance_vec = (pinfo->_normal * (*_current_node - pinfo->_closest_point) +
830 69832 : gapOffset(*_current_node)) *
831 : pinfo->_normal;
832 :
833 69832 : const Real penalty = getPenalty(*_current_node);
834 : RealVectorValue pen_force(penalty * distance_vec);
835 69832 : resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
836 : }
837 1621090 : return _test_secondary[_i][_qp] * resid;
838 :
839 15150245 : case Moose::Primary:
840 15150245 : return _test_primary[_i][_qp] * -resid;
841 : }
842 :
843 : return 0.0;
844 : }
845 :
846 : Real
847 58164767 : MechanicalContactConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
848 : {
849 58164767 : PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
850 :
851 58164767 : const Real penalty = getPenalty(*_current_node);
852 58164767 : const Real penalty_slip = getTangentialPenalty(*_current_node);
853 :
854 58164767 : switch (type)
855 : {
856 0 : default:
857 0 : mooseError("Unhandled ConstraintJacobianType");
858 :
859 2059607 : case Moose::SecondarySecondary:
860 2059607 : switch (_model)
861 : {
862 1630243 : case ContactModel::FRICTIONLESS:
863 1630243 : switch (_formulation)
864 : {
865 : case ContactFormulation::KINEMATIC:
866 : {
867 : RealVectorValue jac_vec;
868 4336464 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
869 : {
870 3213371 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
871 3213371 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
872 3213371 : _var_objects[i]->scalingFactor();
873 : }
874 1123093 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
875 1123093 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
876 1123093 : pinfo->_normal(_component) * pinfo->_normal(_component);
877 : }
878 :
879 507150 : case ContactFormulation::PENALTY:
880 : case ContactFormulation::AUGMENTED_LAGRANGE:
881 507150 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
882 507150 : pinfo->_normal(_component) * pinfo->_normal(_component);
883 :
884 0 : default:
885 0 : mooseError("Invalid contact formulation");
886 : }
887 :
888 214170 : case ContactModel::COULOMB:
889 214170 : switch (_formulation)
890 : {
891 34692 : case ContactFormulation::KINEMATIC:
892 : {
893 34692 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
894 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
895 : {
896 : RealVectorValue jac_vec;
897 102444 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
898 : {
899 68296 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
900 68296 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
901 68296 : _var_objects[i]->scalingFactor();
902 : }
903 34148 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
904 34148 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
905 34148 : pinfo->_normal(_component) * pinfo->_normal(_component);
906 : }
907 : else
908 : {
909 : const Real curr_jac =
910 544 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
911 544 : _connected_dof_indices[_j]) /
912 544 : _var.scalingFactor();
913 544 : return (-curr_jac + _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
914 : }
915 : }
916 :
917 50553 : case ContactFormulation::PENALTY:
918 : {
919 50553 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
920 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
921 10342 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
922 10342 : pinfo->_normal(_component) * pinfo->_normal(_component);
923 : else
924 40211 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp];
925 : }
926 75489 : case ContactFormulation::AUGMENTED_LAGRANGE:
927 : {
928 75489 : Real normal_comp = _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
929 75489 : pinfo->_normal(_component) * pinfo->_normal(_component);
930 :
931 : Real tang_comp = 0.0;
932 75489 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
933 28230 : tang_comp = _phi_secondary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
934 28230 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
935 75489 : return normal_comp + tang_comp;
936 : }
937 :
938 : case ContactFormulation::TANGENTIAL_PENALTY:
939 : {
940 : RealVectorValue jac_vec;
941 160308 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
942 : {
943 106872 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
944 106872 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
945 106872 : _var_objects[i]->scalingFactor();
946 : }
947 53436 : Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
948 53436 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
949 53436 : pinfo->_normal(_component) * pinfo->_normal(_component);
950 :
951 : Real tang_comp = 0.0;
952 53436 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
953 512 : tang_comp = _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
954 512 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
955 :
956 53436 : return normal_comp + tang_comp;
957 : }
958 :
959 0 : default:
960 0 : mooseError("Invalid contact formulation");
961 : }
962 :
963 215194 : case ContactModel::GLUED:
964 215194 : switch (_formulation)
965 : {
966 141305 : case ContactFormulation::KINEMATIC:
967 : {
968 141305 : const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
969 141305 : _connected_dof_indices[_j]) /
970 141305 : _var.scalingFactor();
971 141305 : return (-curr_jac + _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
972 : }
973 :
974 73889 : case ContactFormulation::PENALTY:
975 : case ContactFormulation::AUGMENTED_LAGRANGE:
976 73889 : 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 1465948 : case Moose::SecondaryPrimary:
986 1465948 : switch (_model)
987 : {
988 1138156 : case ContactModel::FRICTIONLESS:
989 1138156 : switch (_formulation)
990 : {
991 770828 : case ContactFormulation::KINEMATIC:
992 : {
993 770828 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
994 :
995 : RealVectorValue jac_vec;
996 2962632 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
997 : {
998 2191804 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
999 4383608 : jac_vec(i) = (*_jacobian)(dof_number,
1000 2191804 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1001 2191804 : _var_objects[i]->scalingFactor();
1002 : }
1003 770828 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1004 770828 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1005 770828 : pinfo->_normal(_component) * pinfo->_normal(_component);
1006 : }
1007 :
1008 367328 : case ContactFormulation::PENALTY:
1009 : case ContactFormulation::AUGMENTED_LAGRANGE:
1010 367328 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1011 367328 : pinfo->_normal(_component) * pinfo->_normal(_component);
1012 :
1013 0 : default:
1014 0 : mooseError("Invalid contact formulation");
1015 : }
1016 :
1017 171240 : case ContactModel::COULOMB:
1018 171240 : switch (_formulation)
1019 : {
1020 26312 : case ContactFormulation::KINEMATIC:
1021 : {
1022 26312 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1023 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1024 : {
1025 25928 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1026 :
1027 : RealVectorValue jac_vec;
1028 77784 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1029 : {
1030 51856 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1031 51856 : jac_vec(i) =
1032 51856 : (*_jacobian)(dof_number,
1033 51856 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1034 51856 : _var_objects[i]->scalingFactor();
1035 : }
1036 25928 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1037 25928 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1038 25928 : pinfo->_normal(_component) * pinfo->_normal(_component);
1039 : }
1040 : else
1041 : {
1042 384 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1043 : const Real curr_jac =
1044 384 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1045 384 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1046 384 : _var.scalingFactor();
1047 384 : return (-curr_jac - _phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
1048 : }
1049 : }
1050 :
1051 40780 : case ContactFormulation::PENALTY:
1052 : {
1053 40780 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1054 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1055 7748 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1056 7748 : pinfo->_normal(_component) * pinfo->_normal(_component);
1057 : else
1058 33032 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp];
1059 : }
1060 62140 : case ContactFormulation::AUGMENTED_LAGRANGE:
1061 : {
1062 62140 : Real normal_comp = -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1063 62140 : pinfo->_normal(_component) * pinfo->_normal(_component);
1064 :
1065 : Real tang_comp = 0.0;
1066 62140 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1067 24560 : tang_comp = -_phi_primary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
1068 24560 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1069 62140 : return normal_comp + tang_comp;
1070 : }
1071 :
1072 42008 : case ContactFormulation::TANGENTIAL_PENALTY:
1073 : {
1074 42008 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1075 :
1076 : RealVectorValue jac_vec;
1077 126024 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1078 : {
1079 84016 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1080 168032 : jac_vec(i) = (*_jacobian)(dof_number,
1081 84016 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1082 84016 : _var_objects[i]->scalingFactor();
1083 : }
1084 42008 : Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1085 42008 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1086 42008 : pinfo->_normal(_component) * pinfo->_normal(_component);
1087 :
1088 : Real tang_comp = 0.0;
1089 42008 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1090 384 : tang_comp = -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1091 384 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1092 :
1093 42008 : return normal_comp + tang_comp;
1094 : }
1095 :
1096 0 : default:
1097 0 : mooseError("Invalid contact formulation");
1098 : }
1099 156552 : case ContactModel::GLUED:
1100 156552 : switch (_formulation)
1101 : {
1102 98028 : case ContactFormulation::KINEMATIC:
1103 : {
1104 98028 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1105 : const Real curr_jac =
1106 98028 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1107 98028 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1108 98028 : _var.scalingFactor();
1109 98028 : return (-curr_jac - _phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
1110 : }
1111 :
1112 58524 : case ContactFormulation::PENALTY:
1113 : case ContactFormulation::AUGMENTED_LAGRANGE:
1114 58524 : 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 32098236 : case Moose::PrimarySecondary:
1125 32098236 : switch (_model)
1126 : {
1127 27493772 : case ContactModel::FRICTIONLESS:
1128 27493772 : switch (_formulation)
1129 : {
1130 : case ContactFormulation::KINEMATIC:
1131 : {
1132 : RealVectorValue jac_vec;
1133 84513264 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1134 : {
1135 63180724 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1136 63180724 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1137 63180724 : _var_objects[i]->scalingFactor();
1138 : }
1139 21332540 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1140 21332540 : _test_primary[_i][_qp];
1141 : }
1142 :
1143 6161232 : case ContactFormulation::PENALTY:
1144 : case ContactFormulation::AUGMENTED_LAGRANGE:
1145 6161232 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1146 6161232 : pinfo->_normal(_component) * pinfo->_normal(_component);
1147 :
1148 0 : default:
1149 0 : mooseError("Invalid contact formulation");
1150 : }
1151 2379144 : case ContactModel::COULOMB:
1152 2379144 : switch (_formulation)
1153 : {
1154 138768 : case ContactFormulation::KINEMATIC:
1155 : {
1156 138768 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1157 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1158 : {
1159 : RealVectorValue jac_vec;
1160 409776 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1161 : {
1162 273184 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1163 273184 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1164 273184 : _var_objects[i]->scalingFactor();
1165 : }
1166 136592 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1167 136592 : _test_primary[_i][_qp];
1168 : }
1169 : else
1170 : {
1171 : const Real secondary_jac =
1172 2176 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1173 2176 : _connected_dof_indices[_j]) /
1174 2176 : _var.scalingFactor();
1175 2176 : return secondary_jac * _test_primary[_i][_qp];
1176 : }
1177 : }
1178 :
1179 747012 : case ContactFormulation::PENALTY:
1180 : {
1181 747012 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1182 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1183 190008 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1184 190008 : pinfo->_normal(_component) * pinfo->_normal(_component);
1185 : else
1186 557004 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp];
1187 : }
1188 1279620 : case ContactFormulation::AUGMENTED_LAGRANGE:
1189 : {
1190 1279620 : Real normal_comp = -_phi_secondary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1191 1279620 : pinfo->_normal(_component) * pinfo->_normal(_component);
1192 :
1193 : Real tang_comp = 0.0;
1194 1279620 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1195 435120 : tang_comp = -_phi_secondary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1196 435120 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1197 1279620 : return normal_comp + tang_comp;
1198 : }
1199 :
1200 : case ContactFormulation::TANGENTIAL_PENALTY:
1201 : {
1202 : RealVectorValue jac_vec;
1203 641232 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1204 : {
1205 427488 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1206 427488 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1207 427488 : _var_objects[i]->scalingFactor();
1208 : }
1209 : Real normal_comp =
1210 213744 : pinfo->_normal(_component) * (pinfo->_normal * jac_vec) * _test_primary[_i][_qp];
1211 :
1212 : Real tang_comp = 0.0;
1213 213744 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1214 2048 : tang_comp = -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1215 2048 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1216 :
1217 213744 : return normal_comp + tang_comp;
1218 : }
1219 :
1220 0 : default:
1221 0 : mooseError("Invalid contact formulation");
1222 : }
1223 :
1224 2225320 : case ContactModel::GLUED:
1225 2225320 : switch (_formulation)
1226 : {
1227 1342940 : case ContactFormulation::KINEMATIC:
1228 : {
1229 : const Real secondary_jac =
1230 1342940 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1231 1342940 : _connected_dof_indices[_j]) /
1232 1342940 : _var.scalingFactor();
1233 1342940 : return secondary_jac * _test_primary[_i][_qp];
1234 : }
1235 :
1236 882380 : case ContactFormulation::PENALTY:
1237 : case ContactFormulation::AUGMENTED_LAGRANGE:
1238 882380 : 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 22540976 : case Moose::PrimaryPrimary:
1249 22540976 : switch (_model)
1250 : {
1251 18970032 : case ContactModel::FRICTIONLESS:
1252 18970032 : switch (_formulation)
1253 : {
1254 : case ContactFormulation::KINEMATIC:
1255 : return 0.0;
1256 :
1257 4541056 : case ContactFormulation::PENALTY:
1258 : case ContactFormulation::AUGMENTED_LAGRANGE:
1259 4541056 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1260 4541056 : pinfo->_normal(_component) * pinfo->_normal(_component);
1261 :
1262 0 : default:
1263 0 : mooseError("Invalid contact formulation");
1264 : }
1265 :
1266 3570944 : case ContactModel::COULOMB:
1267 : case ContactModel::GLUED:
1268 3570944 : switch (_formulation)
1269 : {
1270 : case ContactFormulation::KINEMATIC:
1271 : return 0.0;
1272 :
1273 1306336 : case ContactFormulation::PENALTY:
1274 : {
1275 1306336 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1276 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1277 141712 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1278 141712 : pinfo->_normal(_component) * pinfo->_normal(_component);
1279 : else
1280 1164624 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp];
1281 : }
1282 :
1283 168032 : case ContactFormulation::TANGENTIAL_PENALTY:
1284 : {
1285 : Real tang_comp = 0.0;
1286 168032 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1287 1536 : tang_comp = _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1288 1536 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1289 : return tang_comp; // normal component is zero
1290 : }
1291 :
1292 1041200 : case ContactFormulation::AUGMENTED_LAGRANGE:
1293 : {
1294 1041200 : Real normal_comp = _phi_primary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1295 1041200 : pinfo->_normal(_component) * pinfo->_normal(_component);
1296 :
1297 : Real tang_comp = 0.0;
1298 1041200 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1299 376000 : tang_comp = _phi_primary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1300 376000 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1301 1041200 : 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 11169190 : MechanicalContactConstraint::computeQpOffDiagJacobian(Moose::ConstraintJacobianType type,
1318 : unsigned int jvar)
1319 : {
1320 11169190 : PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
1321 :
1322 11169190 : const Real penalty = getPenalty(*_current_node);
1323 11169190 : const Real penalty_slip = getTangentialPenalty(*_current_node);
1324 :
1325 : unsigned int coupled_component;
1326 : Real normal_component_in_coupled_var_dir = 1.0;
1327 11169190 : if (getCoupledVarComponent(jvar, coupled_component))
1328 11169190 : normal_component_in_coupled_var_dir = pinfo->_normal(coupled_component);
1329 :
1330 11169190 : switch (type)
1331 : {
1332 0 : default:
1333 0 : mooseError("Unhandled ConstraintJacobianType");
1334 :
1335 524806 : case Moose::SecondarySecondary:
1336 524806 : switch (_model)
1337 : {
1338 440000 : case ContactModel::FRICTIONLESS:
1339 440000 : switch (_formulation)
1340 : {
1341 : case ContactFormulation::KINEMATIC:
1342 : {
1343 : RealVectorValue jac_vec;
1344 966672 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1345 : {
1346 718784 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1347 718784 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1348 718784 : _var_objects[i]->scalingFactor();
1349 : }
1350 247888 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1351 247888 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1352 247888 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1353 : }
1354 :
1355 192112 : case ContactFormulation::PENALTY:
1356 : case ContactFormulation::AUGMENTED_LAGRANGE:
1357 192112 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1358 192112 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1359 :
1360 0 : default:
1361 0 : mooseError("Invalid contact formulation");
1362 : }
1363 :
1364 84806 : case ContactModel::COULOMB:
1365 : {
1366 84806 : if ((_formulation == ContactFormulation::KINEMATIC ||
1367 76808 : _formulation == ContactFormulation::TANGENTIAL_PENALTY) &&
1368 76808 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1369 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1370 : {
1371 : RealVectorValue jac_vec;
1372 228120 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1373 : {
1374 152080 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1375 152080 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1376 152080 : _var_objects[i]->scalingFactor();
1377 : }
1378 :
1379 76040 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1380 76040 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1381 76040 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1382 : }
1383 8766 : else if ((_formulation == ContactFormulation::PENALTY) &&
1384 7998 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1385 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1386 672 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1387 672 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1388 8094 : 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 8094 : const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1402 8094 : _connected_dof_indices[_j]);
1403 8094 : 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 388264 : case Moose::SecondaryPrimary:
1418 388264 : switch (_model)
1419 : {
1420 322464 : case ContactModel::FRICTIONLESS:
1421 322464 : switch (_formulation)
1422 : {
1423 186944 : case ContactFormulation::KINEMATIC:
1424 : {
1425 186944 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1426 :
1427 : RealVectorValue jac_vec;
1428 727680 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1429 : {
1430 540736 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1431 1081472 : jac_vec(i) = (*_jacobian)(dof_number,
1432 540736 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1433 540736 : _var_objects[i]->scalingFactor();
1434 : }
1435 186944 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1436 186944 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1437 186944 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1438 : }
1439 :
1440 135520 : case ContactFormulation::PENALTY:
1441 : case ContactFormulation::AUGMENTED_LAGRANGE:
1442 135520 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1443 135520 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1444 :
1445 0 : default:
1446 0 : mooseError("Invalid contact formulation");
1447 : }
1448 :
1449 65800 : case ContactModel::COULOMB:
1450 65800 : if ((_formulation == ContactFormulation::KINEMATIC ||
1451 60560 : _formulation == ContactFormulation::TANGENTIAL_PENALTY) &&
1452 60560 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1453 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1454 : {
1455 59984 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1456 :
1457 : RealVectorValue jac_vec;
1458 179952 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1459 : {
1460 119968 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1461 119968 : jac_vec(i) =
1462 119968 : (*_jacobian)(dof_number, curr_primary_node->dof_number(0, _vars[_component], 0)) /
1463 119968 : _var_objects[i]->scalingFactor();
1464 : }
1465 :
1466 59984 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1467 59984 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1468 59984 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1469 : }
1470 5816 : else if ((_formulation == ContactFormulation::PENALTY) &&
1471 5240 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1472 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1473 448 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1474 448 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1475 5368 : 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 5907032 : case Moose::PrimarySecondary:
1497 5907032 : switch (_model)
1498 : {
1499 5567808 : case ContactModel::FRICTIONLESS:
1500 5567808 : switch (_formulation)
1501 : {
1502 : case ContactFormulation::KINEMATIC:
1503 : {
1504 : RealVectorValue jac_vec;
1505 17582976 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1506 : {
1507 13152256 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1508 13152256 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1509 13152256 : _var_objects[i]->scalingFactor();
1510 : }
1511 4430720 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1512 4430720 : _test_primary[_i][_qp];
1513 : }
1514 :
1515 1137088 : case ContactFormulation::PENALTY:
1516 : case ContactFormulation::AUGMENTED_LAGRANGE:
1517 1137088 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1518 1137088 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1519 :
1520 0 : default:
1521 0 : mooseError("Invalid contact formulation");
1522 : }
1523 339224 : case ContactModel::COULOMB:
1524 339224 : switch (_formulation)
1525 : {
1526 115856 : case ContactFormulation::KINEMATIC:
1527 : {
1528 115856 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1529 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1530 : {
1531 : RealVectorValue jac_vec;
1532 343344 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1533 : {
1534 228896 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1535 228896 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1536 228896 : _var_objects[i]->scalingFactor();
1537 : }
1538 114448 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1539 114448 : _test_primary[_i][_qp];
1540 : }
1541 : else
1542 : {
1543 : const Real secondary_jac =
1544 1408 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1545 1408 : _connected_dof_indices[_j]) /
1546 1408 : _var.scalingFactor();
1547 1408 : return secondary_jac * _test_primary[_i][_qp];
1548 : }
1549 : }
1550 :
1551 31992 : case ContactFormulation::PENALTY:
1552 : {
1553 31992 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1554 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1555 2688 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1556 2688 : 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 191376 : case ContactFormulation::TANGENTIAL_PENALTY:
1572 : {
1573 191376 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1574 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1575 : {
1576 : RealVectorValue jac_vec;
1577 569136 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1578 : {
1579 379424 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1580 379424 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1581 379424 : _var_objects[i]->scalingFactor();
1582 : }
1583 189712 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1584 189712 : _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 4349088 : case Moose::PrimaryPrimary:
1619 4349088 : switch (_model)
1620 : {
1621 4085888 : case ContactModel::FRICTIONLESS:
1622 4085888 : switch (_formulation)
1623 : {
1624 : case ContactFormulation::KINEMATIC:
1625 : return 0.0;
1626 :
1627 763264 : case ContactFormulation::PENALTY:
1628 : case ContactFormulation::AUGMENTED_LAGRANGE:
1629 763264 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1630 763264 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1631 :
1632 0 : default:
1633 0 : mooseError("Invalid contact formulation");
1634 : }
1635 :
1636 263200 : case ContactModel::COULOMB:
1637 : {
1638 263200 : 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 263200 : else if (_formulation == ContactFormulation::PENALTY &&
1650 20960 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1651 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1652 1792 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1653 1792 : 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 2138017 : MechanicalContactConstraint::gapOffset(const Node & node)
1688 : {
1689 : Real val = 0;
1690 :
1691 2138017 : if (_has_secondary_gap_offset)
1692 83808 : val += _secondary_gap_offset_var->getNodalValue(node);
1693 :
1694 2138017 : if (_has_mapped_primary_gap_offset)
1695 83808 : val += _mapped_primary_gap_offset_var->getNodalValue(node);
1696 :
1697 2138017 : return val;
1698 : }
1699 :
1700 : Real
1701 28352946 : MechanicalContactConstraint::nodalArea(const Node & node)
1702 : {
1703 28352946 : dof_id_type dof = node.dof_number(_aux_system.number(), _nodal_area_var->number(), 0);
1704 :
1705 28352946 : Real area = (*_aux_solution)(dof);
1706 28352946 : 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 28352946 : return area;
1715 : }
1716 :
1717 : Real
1718 71100279 : MechanicalContactConstraint::getPenalty(const Node & node)
1719 : {
1720 71100279 : Real penalty = _penalty;
1721 71100279 : if (_normalize_penalty)
1722 14176840 : penalty *= nodalArea(node);
1723 71100279 : return penalty * MathUtils::pow(_penalty_multiplier, _lagrangian_iteration_number);
1724 : }
1725 :
1726 : Real
1727 70091691 : MechanicalContactConstraint::getTangentialPenalty(const Node & node)
1728 : {
1729 70091691 : Real penalty = _penalty_tangential;
1730 70091691 : if (_normalize_penalty)
1731 14111481 : penalty *= nodalArea(node);
1732 :
1733 70091691 : return penalty * MathUtils::pow(_penalty_multiplier, _lagrangian_iteration_number);
1734 : }
1735 :
1736 : void
1737 155663 : MechanicalContactConstraint::computeJacobian()
1738 : {
1739 155663 : getConnectedDofIndices(_var.number());
1740 :
1741 155663 : prepareMatrixTagNeighbor(
1742 155663 : _assembly, _primary_var.number(), _var.number(), Moose::NeighborNeighbor, _Knn);
1743 :
1744 155663 : _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
1745 :
1746 311326 : for (_i = 0; _i < _test_secondary.size(); _i++)
1747 : // Loop over the connected dof indices so we can get all the jacobian contributions
1748 2215270 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1749 2059607 : _Kee(_i, _j) += computeQpJacobian(Moose::SecondarySecondary);
1750 :
1751 155663 : if (_primary_secondary_jacobian)
1752 : {
1753 153239 : prepareMatrixTagNeighbor(_assembly, _var.number(), _var.number(), Moose::ElementNeighbor, _Ken);
1754 153239 : if (_Ken.m() && _Ken.n())
1755 : {
1756 306478 : for (_i = 0; _i < _test_secondary.size(); _i++)
1757 1619187 : for (_j = 0; _j < _phi_primary.size(); _j++)
1758 1465948 : _Ken(_i, _j) += computeQpJacobian(Moose::SecondaryPrimary);
1759 153239 : accumulateTaggedLocalMatrix(
1760 153239 : _assembly, _var.number(), _var.number(), Moose::ElementNeighbor, _Ken);
1761 : }
1762 :
1763 153239 : _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
1764 1619187 : for (_i = 0; _i < _test_primary.size(); _i++)
1765 : // Loop over the connected dof indices so we can get all the jacobian contributions
1766 33564184 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1767 32098236 : _Kne(_i, _j) += computeQpJacobian(Moose::PrimarySecondary);
1768 : }
1769 :
1770 155663 : if (_Knn.m() && _Knn.n())
1771 : {
1772 1641003 : for (_i = 0; _i < _test_primary.size(); _i++)
1773 24026316 : for (_j = 0; _j < _phi_primary.size(); _j++)
1774 22540976 : _Knn(_i, _j) += computeQpJacobian(Moose::PrimaryPrimary);
1775 155663 : accumulateTaggedLocalMatrix(
1776 155663 : _assembly, _primary_var.number(), _var.number(), Moose::NeighborNeighbor, _Knn);
1777 : }
1778 155663 : }
1779 :
1780 : void
1781 56634 : MechanicalContactConstraint::computeOffDiagJacobian(const unsigned int jvar_num)
1782 : {
1783 56634 : getConnectedDofIndices(jvar_num);
1784 :
1785 56634 : _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
1786 :
1787 56634 : prepareMatrixTagNeighbor(
1788 56634 : _assembly, _primary_var.number(), jvar_num, Moose::NeighborNeighbor, _Knn);
1789 :
1790 113268 : for (_i = 0; _i < _test_secondary.size(); _i++)
1791 : // Loop over the connected dof indices so we can get all the jacobian contributions
1792 581440 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1793 524806 : _Kee(_i, _j) += computeQpOffDiagJacobian(Moose::SecondarySecondary, jvar_num);
1794 :
1795 56634 : if (_primary_secondary_jacobian)
1796 : {
1797 56634 : prepareMatrixTagNeighbor(_assembly, _var.number(), jvar_num, Moose::ElementNeighbor, _Ken);
1798 113268 : for (_i = 0; _i < _test_secondary.size(); _i++)
1799 444898 : for (_j = 0; _j < _phi_primary.size(); _j++)
1800 388264 : _Ken(_i, _j) += computeQpOffDiagJacobian(Moose::SecondaryPrimary, jvar_num);
1801 56634 : accumulateTaggedLocalMatrix(_assembly, _var.number(), jvar_num, Moose::ElementNeighbor, _Ken);
1802 :
1803 56634 : _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
1804 56634 : if (_Kne.m() && _Kne.n())
1805 444898 : for (_i = 0; _i < _test_primary.size(); _i++)
1806 : // Loop over the connected dof indices so we can get all the jacobian contributions
1807 6295296 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1808 5907032 : _Kne(_i, _j) += computeQpOffDiagJacobian(Moose::PrimarySecondary, jvar_num);
1809 : }
1810 :
1811 444898 : for (_i = 0; _i < _test_primary.size(); _i++)
1812 4737352 : for (_j = 0; _j < _phi_primary.size(); _j++)
1813 4349088 : _Knn(_i, _j) += computeQpOffDiagJacobian(Moose::PrimaryPrimary, jvar_num);
1814 56634 : accumulateTaggedLocalMatrix(
1815 56634 : _assembly, _primary_var.number(), jvar_num, Moose::NeighborNeighbor, _Knn);
1816 56634 : }
1817 :
1818 : void
1819 212297 : MechanicalContactConstraint::getConnectedDofIndices(unsigned int var_num)
1820 : {
1821 : unsigned int component;
1822 212297 : if (getCoupledVarComponent(var_num, component) || _non_displacement_vars_jacobian)
1823 : {
1824 212297 : if (_primary_secondary_jacobian && _connected_secondary_nodes_jacobian)
1825 209873 : NodeFaceConstraint::getConnectedDofIndices(var_num);
1826 : else
1827 : {
1828 2424 : _connected_dof_indices.clear();
1829 2424 : MooseVariableFEBase & var = _sys.getVariable(0, var_num);
1830 2424 : _connected_dof_indices.push_back(var.nodalDofIndex());
1831 : }
1832 : }
1833 :
1834 212297 : _phi_secondary.resize(_connected_dof_indices.size());
1835 :
1836 212297 : dof_id_type current_node_var_dof_index = _sys.getVariable(0, var_num).nodalDofIndex();
1837 212297 : _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 2796710 : for (unsigned int j = 0; j < _connected_dof_indices.size(); j++)
1843 : {
1844 2584413 : _phi_secondary[j].resize(1);
1845 :
1846 2584413 : if (_connected_dof_indices[j] == current_node_var_dof_index)
1847 212297 : _phi_secondary[j][_qp] = 1.0;
1848 : else
1849 2372116 : _phi_secondary[j][_qp] = 0.0;
1850 : }
1851 212297 : }
1852 :
1853 : bool
1854 11381487 : MechanicalContactConstraint::getCoupledVarComponent(unsigned int var_num, unsigned int & component)
1855 : {
1856 11381487 : component = std::numeric_limits<unsigned int>::max();
1857 : bool coupled_var_is_disp_var = false;
1858 21723906 : for (const auto i : make_range(Moose::dim))
1859 : {
1860 21723906 : if (var_num == _vars[i])
1861 : {
1862 : coupled_var_is_disp_var = true;
1863 11381487 : component = i;
1864 11381487 : break;
1865 : }
1866 : }
1867 :
1868 11381487 : return coupled_var_is_disp_var;
1869 : }
1870 :
1871 : void
1872 70671 : MechanicalContactConstraint::residualEnd()
1873 : {
1874 70671 : if (_component == 0 && (_print_contact_nodes || _contact_linesearch))
1875 : {
1876 70 : _communicator.set_union(_current_contact_state);
1877 70 : 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 70 : if (_contact_linesearch)
1887 70 : _contact_linesearch->insertSet(_current_contact_state);
1888 : _old_contact_state.swap(_current_contact_state);
1889 : _current_contact_state.clear();
1890 : }
1891 70671 : }
|
