https://mooseframework.inl.gov
Loading...
Searching...
No Matches
MechanicalContactConstraint.C
Go to the documentation of this file.
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
30registerMooseObject("ContactApp", MechanicalContactConstraint);
31
32const unsigned int MechanicalContactConstraint::_no_iterations = 0;
33
34InputParameters
35MechanicalContactConstraint::validParams()
36{
37 InputParameters params = NodeFaceConstraint::validParams();
38 params += ContactAction::commonParameters();
39
40 params.addRequiredParam<BoundaryName>("boundary", "The primary boundary");
41 params.addParam<BoundaryName>("secondary", "The secondary boundary");
42 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 params.addCoupledVar(
48 "displacements",
49 "The displacements appropriate for the simulation geometry and coordinate system");
50
51 params.addCoupledVar("secondary_gap_offset", "offset to the gap distance from secondary side");
52 params.addCoupledVar("mapped_primary_gap_offset",
53 "offset to the gap distance mapped from primary side");
54 params.addRequiredCoupledVar("nodal_area", "The nodal area");
55
56 params.set<bool>("use_displaced_mesh") = true;
57 params.addParam<Real>(
58 "penalty",
59 1e8,
60 "The penalty to apply. This can vary depending on the stiffness of your materials");
61 params.addParam<Real>("penalty_multiplier",
62 1.0,
63 "The growth factor for the penalty applied at the end of each augmented "
64 "Lagrange update iteration");
65 params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
66 params.addParam<Real>("tangential_tolerance",
67 "Tangential distance to extend edges of contact surfaces");
68 params.addParam<Real>(
69 "capture_tolerance", 0, "Normal distance from surface within which nodes are captured");
70
71 params.addParam<Real>("tension_release",
72 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 params.addParam<bool>(
78 "normalize_penalty",
79 false,
80 "Whether to normalize the penalty parameter with the nodal area for penalty contact.");
81 params.addParam<bool>(
82 "primary_secondary_jacobian",
83 true,
84 "Whether to include jacobian entries coupling primary and secondary nodes.");
85 params.addParam<bool>(
86 "connected_secondary_nodes_jacobian",
87 true,
88 "Whether to include jacobian entries coupling nodes connected to secondary nodes.");
89 params.addParam<bool>("non_displacement_variables_jacobian",
90 true,
91 "Whether to include jacobian entries coupling with variables that are not "
92 "displacement variables.");
93 params.addParam<unsigned int>("stick_lock_iterations",
94 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 params.addParam<Real>("stick_unlock_factor",
98 1.5,
99 "Factor by which frictional capacity must be "
100 "exceeded to permit stick-locked node to slip "
101 "again.");
102 params.addParam<Real>("al_penetration_tolerance",
103 "The tolerance of the penetration for augmented Lagrangian method.");
104 params.addParam<Real>("al_incremental_slip_tolerance",
105 "The tolerance of the incremental slip for augmented Lagrangian method.");
106
107 params.addParam<Real>("al_frictional_force_tolerance",
108 "The tolerance of the frictional force for augmented Lagrangian method.");
109 params.addParam<bool>(
110 "print_contact_nodes", false, "Whether to print the number of nodes in contact.");
111
112 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 return params;
119}
120
121Threads::spin_mutex MechanicalContactConstraint::_contact_set_mutex;
122
123MechanicalContactConstraint::MechanicalContactConstraint(const InputParameters & parameters)
124 : NodeFaceConstraint(parameters),
125 _displaced_problem(parameters.get<FEProblemBase *>("_fe_problem_base")->getDisplacedProblem()),
126 _component(getParam<unsigned int>("component")),
127 _model(getParam<MooseEnum>("model").getEnum<ContactModel>()),
128 _formulation(getParam<MooseEnum>("formulation").getEnum<ContactFormulation>()),
129 _normalize_penalty(getParam<bool>("normalize_penalty")),
130 _penalty(getParam<Real>("penalty")),
131 _penalty_multiplier(getParam<Real>("penalty_multiplier")),
132 _friction_coefficient(getParam<Real>("friction_coefficient")),
133 _tension_release(getParam<Real>("tension_release")),
134 _capture_tolerance(getParam<Real>("capture_tolerance")),
135 _stick_lock_iterations(getParam<unsigned int>("stick_lock_iterations")),
136 _stick_unlock_factor(getParam<Real>("stick_unlock_factor")),
137 _update_stateful_data(true),
138 _residual_copy(_sys.residualGhosted()),
139 _mesh_dimension(_mesh.dimension()),
140 _vars(3, libMesh::invalid_uint),
141 _var_objects(3, nullptr),
142 _has_secondary_gap_offset(isCoupled("secondary_gap_offset")),
143 _secondary_gap_offset_var(_has_secondary_gap_offset ? getVar("secondary_gap_offset", 0)
144 : nullptr),
145 _has_mapped_primary_gap_offset(isCoupled("mapped_primary_gap_offset")),
146 _mapped_primary_gap_offset_var(
147 _has_mapped_primary_gap_offset ? getVar("mapped_primary_gap_offset", 0) : nullptr),
148 _nodal_area_var(getVar("nodal_area", 0)),
149 _aux_system(_nodal_area_var->sys()),
150 _aux_solution(_aux_system.currentSolution()),
151 _primary_secondary_jacobian(getParam<bool>("primary_secondary_jacobian")),
152 _connected_secondary_nodes_jacobian(getParam<bool>("connected_secondary_nodes_jacobian")),
153 _non_displacement_vars_jacobian(getParam<bool>("non_displacement_variables_jacobian")),
154 _contact_linesearch(dynamic_cast<ContactLineSearchBase *>(_subproblem.getLineSearch())),
155 _print_contact_nodes(getParam<bool>("print_contact_nodes")),
156 _augmented_lagrange_problem(
157 dynamic_cast<AugmentedLagrangianContactProblemInterface *>(&_fe_problem)),
158 _lagrangian_iteration_number(_augmented_lagrange_problem
159 ? _augmented_lagrange_problem->getLagrangianIterationNumber()
160 : _no_iterations)
161{
162 _overwrite_secondary_residual = false;
163
164 if (isParamValid("displacements"))
165 {
166 // modern parameter scheme for displacements
167 for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
168 {
169 _vars[i] = coupled("displacements", i);
170 _var_objects[i] = getVar("displacements", i);
171 }
172 }
173
174 if (parameters.isParamValid("tangential_tolerance"))
175 _penetration_locator.setTangentialTolerance(getParam<Real>("tangential_tolerance"));
176
177 if (parameters.isParamValid("normal_smoothing_distance"))
178 _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
179
180 if (parameters.isParamValid("normal_smoothing_method"))
181 _penetration_locator.setNormalSmoothingMethod(
182 parameters.get<std::string>("normal_smoothing_method"));
183
184 if (_formulation == ContactFormulation::TANGENTIAL_PENALTY && _model != ContactModel::COULOMB)
185 mooseError("The 'tangential_penalty' formulation can only be used with the 'coulomb' model");
186
187 if (_model == ContactModel::GLUED)
188 _penetration_locator.setUpdate(false);
189
190 if (_friction_coefficient < 0)
191 mooseError("The friction coefficient must be nonnegative");
192
193 // set _penalty_tangential to the value of _penalty for now
194 _penalty_tangential = _penalty;
195
196 if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
197 {
198 if (_model == ContactModel::GLUED)
199 mooseError("The Augmented Lagrangian contact formulation does not support GLUED case.");
200
201 if (!_augmented_lagrange_problem)
202 mooseError("The Augmented Lagrangian contact formulation must use "
203 "AugmentedLagrangianContactProblem.");
204
205 if (!parameters.isParamValid("al_penetration_tolerance"))
206 mooseError("For Augmented Lagrangian contact, al_penetration_tolerance must be provided.");
207 else
208 _al_penetration_tolerance = parameters.get<Real>("al_penetration_tolerance");
209
210 if (_model != ContactModel::FRICTIONLESS)
211 {
212 if (!parameters.isParamValid("al_incremental_slip_tolerance") ||
213 !parameters.isParamValid("al_frictional_force_tolerance"))
214 {
215 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 _al_incremental_slip_tolerance = parameters.get<Real>("al_incremental_slip_tolerance");
222 _al_frictional_force_tolerance = parameters.get<Real>("al_frictional_force_tolerance");
223 }
224 }
225 }
226}
227
228void
229MechanicalContactConstraint::timestepSetup()
230{
231 if (_component == 0)
232 {
233 updateContactStatefulData(/* beginning_of_step = */ true);
234 if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
235 updateAugmentedLagrangianMultiplier(/* beginning_of_step = */ true);
236
237 _update_stateful_data = false;
238
239 if (_contact_linesearch)
240 _contact_linesearch->reset();
241 }
242}
243
244void
245MechanicalContactConstraint::jacobianSetup()
246{
247 if (_component == 0)
248 {
249 if (_update_stateful_data)
250 updateContactStatefulData(/* beginning_of_step = */ false);
251 _update_stateful_data = true;
252 }
253}
254
255void
256MechanicalContactConstraint::updateAugmentedLagrangianMultiplier(bool beginning_of_step)
257{
258 for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
259 {
260 if (!pinfo)
261 continue;
262
263 const Node & node = _mesh.nodeRef(secondary_node_num);
264 if (node.n_comp(_sys.number(), _vars[_component]) < 1)
265 continue;
266
267 const Real distance = pinfo->_normal * (pinfo->_closest_point - node) - gapOffset(node);
268
269 if (beginning_of_step && _model == ContactModel::COULOMB)
270 {
271 pinfo->_lagrange_multiplier_slip.zero();
272 if (pinfo->isCaptured())
273 pinfo->_mech_status = PenetrationInfo::MS_STICKING;
274 }
275
276 if (pinfo->isCaptured())
277 {
278 if (_model == ContactModel::FRICTIONLESS)
279 pinfo->_lagrange_multiplier -= getPenalty(node) * distance;
280
281 if (_model == ContactModel::COULOMB)
282 {
283 if (!beginning_of_step)
284 {
285 Real penalty = getPenalty(node);
286 RealVectorValue pen_force_normal =
287 penalty * (-distance) * pinfo->_normal + pinfo->_lagrange_multiplier * pinfo->_normal;
288
289 // update normal lagrangian multiplier
290 pinfo->_lagrange_multiplier += penalty * (-distance);
291
292 // Frictional capacity
293 const Real capacity(_friction_coefficient * (pen_force_normal * pinfo->_normal < 0
294 ? -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 Real penalty_slip = getTangentialPenalty(node);
302
303 RealVectorValue inc_pen_force_tangential =
304 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 const Real tan_mag(contact_force_tangential.norm());
311
312 if (tan_mag > capacity * (_al_frictional_force_tolerance + 1.0))
313 {
314 pinfo->_lagrange_multiplier_slip =
315 -tau_old + capacity * contact_force_tangential / tan_mag;
316 if (MooseUtils::absoluteFuzzyEqual(capacity, 0.0))
317 pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
318 else
319 pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
320 }
321 else
322 {
323 pinfo->_mech_status = PenetrationInfo::MS_STICKING;
324 pinfo->_lagrange_multiplier_slip += penalty_slip * tangential_inc_slip;
325 }
326 }
327 }
328 }
329 }
330}
331
332bool
333MechanicalContactConstraint::AugmentedLagrangianContactConverged()
334{
335 Real contactResidual = 0.0;
336 unsigned int converged = 0;
337
338 for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
339 {
340 if (!pinfo)
341 continue;
342
343 const auto & node = _mesh.nodeRef(secondary_node_num);
344
345 // Skip this pinfo if there are no DOFs on this node.
346 if (node.n_comp(_sys.number(), _vars[_component]) < 1)
347 continue;
348
349 const Real distance = pinfo->_normal * (pinfo->_closest_point - node) - gapOffset(node);
350
351 if (pinfo->isCaptured())
352 {
353 if (contactResidual < std::abs(distance))
354 contactResidual = std::abs(distance);
355
356 // penetration < tol
357 if (contactResidual > _al_penetration_tolerance)
358 {
359 converged = 1;
360 break;
361 }
362
363 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 const Real tan_mag(contact_force_tangential.norm());
373 const Real tangential_inc_slip_mag = tangential_inc_slip.norm();
374
375 RealVectorValue distance_vec =
376 (pinfo->_normal * (node - pinfo->_closest_point) + gapOffset(node)) * pinfo->_normal;
377
378 Real penalty = getPenalty(node);
379 RealVectorValue pen_force_normal =
380 penalty * distance_vec + pinfo->_lagrange_multiplier * pinfo->_normal;
381
382 // Frictional capacity
383 Real capacity(_friction_coefficient * (pen_force_normal * pinfo->_normal < 0
384 ? -pen_force_normal * pinfo->_normal
385 : 0.0));
386
387 // incremental slip <= tol for all pinfo_pair such that tan_mag < capacity
388 if (MooseUtils::absoluteFuzzyLessThan(tan_mag, capacity) &&
389 pinfo->_mech_status == PenetrationInfo::MS_STICKING)
390 {
391 if (MooseUtils::absoluteFuzzyGreaterThan(tangential_inc_slip_mag,
392 _al_incremental_slip_tolerance))
393 {
394 converged = 2;
395 break;
396 }
397 }
398
399 // for all pinfo_pair, tag_mag should be less than (1 + tol) * (capacity + tol)
400 if (tan_mag >
401 (1 + _al_frictional_force_tolerance) * (capacity + _al_frictional_force_tolerance))
402 {
403 converged = 3;
404 break;
405 }
406 }
407 }
408 }
409
410 _communicator.max(converged);
411
412 if (converged == 1)
413 _console << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied "
414 << std::endl;
415 else if (converged == 2)
416 _console << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied "
417 << std::endl;
418 else if (converged == 3)
419 _console << "The Augmented Lagrangian contact frictional force enforcement is NOT satisfied "
420 << std::endl;
421 else
422 return true;
423
424 return false;
425}
426
427void
428MechanicalContactConstraint::updateContactStatefulData(bool beginning_of_step)
429{
430 for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
431 {
432 if (!pinfo)
433 continue;
434
435 const Node & node = _mesh.nodeRef(secondary_node_num);
436 if (node.n_comp(_sys.number(), _vars[_component]) < 1)
437 continue;
438
439 if (beginning_of_step)
440 {
441 if (_app.getExecutioner()->lastSolveConverged())
442 {
443 pinfo->_contact_force_old = pinfo->_contact_force;
444 pinfo->_accumulated_slip_old = pinfo->_accumulated_slip;
445 pinfo->_frictional_energy_old = pinfo->_frictional_energy;
446 pinfo->_mech_status_old = pinfo->_mech_status;
447 }
448 else if (pinfo->_mech_status_old == PenetrationInfo::MS_NO_CONTACT &&
449 pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
450 {
451 // The penetration info object could be based on a bad state so delete it
452 delete pinfo;
453 pinfo = NULL;
454 continue;
455 }
456
457 pinfo->_locked_this_step = 0;
458 pinfo->_stick_locked_this_step = 0;
459 pinfo->_starting_elem = pinfo->_elem;
460 pinfo->_starting_side_num = pinfo->_side_num;
461 pinfo->_starting_closest_point_ref = pinfo->_closest_point_ref;
462 }
463 pinfo->_incremental_slip_prev_iter = pinfo->_incremental_slip;
464 }
465}
466
467bool
468MechanicalContactConstraint::shouldApply()
469{
470 bool in_contact = false;
471
472 std::map<dof_id_type, PenetrationInfo *>::iterator found =
473 _penetration_locator._penetration_info.find(_current_node->id());
474 if (found != _penetration_locator._penetration_info.end())
475 {
476 PenetrationInfo * pinfo = found->second;
477 if (pinfo != NULL)
478 {
479 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 if (_component == 0)
484 computeContactForce(*_current_node, pinfo, is_nonlinear);
485
486 if (pinfo->isCaptured())
487 {
488 in_contact = true;
489 if (is_nonlinear)
490 {
491 Threads::spin_mutex::scoped_lock lock(_contact_set_mutex);
492 _current_contact_state.insert(_current_node->id());
493 }
494 }
495 }
496 }
497
498 return in_contact;
499}
500
501void
502MechanicalContactConstraint::computeContactForce(const Node & node,
503 PenetrationInfo * pinfo,
504 bool update_contact_set)
505{
506 // Build up residual vector
507 RealVectorValue res_vec;
508 for (unsigned int i = 0; i < _mesh_dimension; ++i)
509 {
510 dof_id_type dof_number = node.dof_number(0, _vars[i], 0);
511 res_vec(i) = _residual_copy(dof_number) / _var_objects[i]->scalingFactor();
512 }
513
514 RealVectorValue distance_vec(node - pinfo->_closest_point);
515 if (distance_vec.norm() != 0)
516 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 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 if (_formulation == ContactFormulation::KINEMATIC ||
534 _formulation == ContactFormulation::TANGENTIAL_PENALTY)
535 ++pinfo->_locked_this_step;
536 }
537
538 if (!pinfo->isCaptured())
539 return;
540
541 const Real penalty = getPenalty(node);
542 const Real penalty_slip = getTangentialPenalty(node);
543
544 RealVectorValue pen_force(penalty * distance_vec);
545
546 switch (_model)
547 {
548 case ContactModel::FRICTIONLESS:
549 switch (_formulation)
550 {
551 case ContactFormulation::KINEMATIC:
552 pinfo->_contact_force = -pinfo->_normal * (pinfo->_normal * res_vec);
553 break;
554
555 case ContactFormulation::PENALTY:
556 pinfo->_contact_force = pinfo->_normal * (pinfo->_normal * pen_force);
557 break;
558
559 case ContactFormulation::AUGMENTED_LAGRANGE:
560 pinfo->_contact_force =
561 (pinfo->_normal *
562 (pinfo->_normal * (pen_force + pinfo->_lagrange_multiplier * pinfo->_normal)));
563 break;
564
565 default:
566 mooseError("Invalid contact formulation");
567 break;
568 }
569 pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
570 break;
571 case ContactModel::COULOMB:
572 switch (_formulation)
573 {
574 case ContactFormulation::KINEMATIC:
575 {
576 // Frictional capacity
577 const Real capacity(_friction_coefficient *
578 (res_vec * pinfo->_normal > 0 ? res_vec * pinfo->_normal : 0));
579
580 // Normal and tangential components of predictor force
581 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 const Real tan_mag(contact_force_tangential.norm());
592 const Real tangential_inc_slip_mag = tangential_inc_slip.norm();
593 const Real slip_tol = capacity / penalty;
594 pinfo->_slip_tol = slip_tol;
595
596 if ((tangential_inc_slip_mag > slip_tol || tan_mag > capacity) &&
597 (pinfo->_stick_locked_this_step < _stick_lock_iterations ||
598 tan_mag > capacity * _stick_unlock_factor))
599 {
600 if (pinfo->_stick_locked_this_step >= _stick_lock_iterations)
601 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 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 if (slip_dot_tang_force < capacity)
612 slipped_too_far = true;
613 }
614
615 if (slipped_too_far) // slip back along slip increment
616 pinfo->_contact_force = contact_force_normal + capacity * slip_inc_direction;
617 else
618 {
619 if (tan_mag > 0) // slip along tangential force direction
620 pinfo->_contact_force =
621 contact_force_normal + capacity * contact_force_tangential / tan_mag;
622 else // treat as frictionless
623 pinfo->_contact_force = contact_force_normal;
624 }
625 if (capacity == 0)
626 pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
627 else
628 pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
629 }
630 else
631 {
632 if (pinfo->_mech_status != PenetrationInfo::MS_STICKING &&
633 pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
634 ++pinfo->_stick_locked_this_step;
635 pinfo->_mech_status = PenetrationInfo::MS_STICKING;
636 }
637 break;
638 }
639
640 case ContactFormulation::PENALTY:
641 {
642 distance_vec =
643 pinfo->_incremental_slip +
644 (pinfo->_normal * (node - pinfo->_closest_point) + gapOffset(node)) * pinfo->_normal;
645 pen_force = penalty * distance_vec;
646
647 // Frictional capacity
648 const Real capacity(_friction_coefficient *
649 (pen_force * pinfo->_normal < 0 ? -pen_force * pinfo->_normal : 0));
650
651 // Elastic predictor
652 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 const Real tan_mag(contact_force_tangential.norm());
661
662 if (tan_mag > capacity)
663 {
664 pinfo->_contact_force =
665 contact_force_normal + capacity * contact_force_tangential / tan_mag;
666 if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
667 pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
668 else
669 pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
670 }
671 else
672 pinfo->_mech_status = PenetrationInfo::MS_STICKING;
673 break;
674 }
675
676 case ContactFormulation::AUGMENTED_LAGRANGE:
677 {
678 distance_vec =
679 (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 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
696 pinfo->_contact_force =
697 contact_force_normal + contact_force_tangential + inc_pen_force_tangential;
698 else
699 pinfo->_contact_force = contact_force_normal + contact_force_tangential;
700
701 break;
702 }
703
704 case ContactFormulation::TANGENTIAL_PENALTY:
705 {
706 // Frictional capacity (kinematic formulation)
707 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 const Real tan_mag(contact_force_tangential.norm());
721
722 if (tan_mag > capacity)
723 {
724 pinfo->_contact_force =
725 contact_force_normal + capacity * contact_force_tangential / tan_mag;
726 if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
727 pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
728 else
729 pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
730 }
731 else
732 {
733 pinfo->_contact_force = contact_force_normal + contact_force_tangential;
734 pinfo->_mech_status = PenetrationInfo::MS_STICKING;
735 }
736 break;
737 }
738
739 default:
740 mooseError("Invalid contact formulation");
741 break;
742 }
743 break;
744
745 case ContactModel::GLUED:
746 switch (_formulation)
747 {
748 case ContactFormulation::KINEMATIC:
749 pinfo->_contact_force = -res_vec;
750 break;
751
752 case ContactFormulation::PENALTY:
753 pinfo->_contact_force = pen_force;
754 break;
755
756 case ContactFormulation::AUGMENTED_LAGRANGE:
757 pinfo->_contact_force =
758 pen_force + pinfo->_lagrange_multiplier * distance_vec / distance_vec.norm();
759 break;
760
761 default:
762 mooseError("Invalid contact formulation");
763 break;
764 }
765 pinfo->_mech_status = PenetrationInfo::MS_STICKING;
766 break;
767
768 default:
769 mooseError("Invalid or unavailable contact model");
770 break;
771 }
772
773 // Release
774 if (update_contact_set && _model != ContactModel::GLUED && pinfo->isCaptured() &&
775 !newly_captured && _tension_release >= 0.0 &&
776 (_contact_linesearch ? true : pinfo->_locked_this_step < 2))
777 {
778 const Real contact_pressure = -(pinfo->_normal * pinfo->_contact_force) / nodalArea(node);
779 if (-contact_pressure >= _tension_release)
780 {
781 pinfo->release();
782 pinfo->_contact_force.zero();
783 }
784 }
785}
786
787Real
792
793Real
794MechanicalContactConstraint::computeQpResidual(Moose::ConstraintType type)
795{
796 PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
797 Real resid = pinfo->_contact_force(_component);
798 switch (type)
799 {
800 case Moose::Secondary:
801 if (_formulation == ContactFormulation::KINEMATIC)
802 {
803 RealVectorValue distance_vec(*_current_node - pinfo->_closest_point);
804 if (distance_vec.norm() != 0)
805 distance_vec += gapOffset(*_current_node) * pinfo->_normal * distance_vec.unit() *
806 distance_vec.unit();
807
808 const Real penalty = getPenalty(*_current_node);
809 RealVectorValue pen_force(penalty * distance_vec);
810
811 if (_model == ContactModel::FRICTIONLESS)
812 resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
813 else if (_model == ContactModel::COULOMB)
814 {
815 distance_vec = distance_vec - pinfo->_incremental_slip;
816 pen_force = penalty * distance_vec;
817 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
818 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
819 resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
820 else
821 resid += pen_force(_component);
822 }
823 else if (_model == ContactModel::GLUED)
824 resid += pen_force(_component);
825 }
826 else if (_formulation == ContactFormulation::TANGENTIAL_PENALTY &&
827 _model == ContactModel::COULOMB)
828 {
829 RealVectorValue distance_vec = (pinfo->_normal * (*_current_node - pinfo->_closest_point) +
830 gapOffset(*_current_node)) *
831 pinfo->_normal;
832
833 const Real penalty = getPenalty(*_current_node);
834 RealVectorValue pen_force(penalty * distance_vec);
835 resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
836 }
837 return _test_secondary[_i][_qp] * resid;
838
839 case Moose::Primary:
840 return _test_primary[_i][_qp] * -resid;
841 }
842
843 return 0.0;
844}
845
846Real
847MechanicalContactConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
848{
849 PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
850
851 const Real penalty = getPenalty(*_current_node);
852 const Real penalty_slip = getTangentialPenalty(*_current_node);
853
854 switch (type)
855 {
856 default:
857 mooseError("Unhandled ConstraintJacobianType");
858
859 case Moose::SecondarySecondary:
860 switch (_model)
861 {
862 case ContactModel::FRICTIONLESS:
863 switch (_formulation)
864 {
865 case ContactFormulation::KINEMATIC:
866 {
867 RealVectorValue jac_vec;
868 for (unsigned int i = 0; i < _mesh_dimension; ++i)
869 {
870 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
871 jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
872 _var_objects[i]->scalingFactor();
873 }
874 return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
875 (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
876 pinfo->_normal(_component) * pinfo->_normal(_component);
877 }
878
879 case ContactFormulation::PENALTY:
880 case ContactFormulation::AUGMENTED_LAGRANGE:
881 return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
882 pinfo->_normal(_component) * pinfo->_normal(_component);
883
884 default:
885 mooseError("Invalid contact formulation");
886 }
887
888 case ContactModel::COULOMB:
889 switch (_formulation)
890 {
891 case ContactFormulation::KINEMATIC:
892 {
893 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
894 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
895 {
896 RealVectorValue jac_vec;
897 for (unsigned int i = 0; i < _mesh_dimension; ++i)
898 {
899 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
900 jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
901 _var_objects[i]->scalingFactor();
902 }
903 return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
904 (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
905 pinfo->_normal(_component) * pinfo->_normal(_component);
906 }
907 else
908 {
909 const Real curr_jac =
910 (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
911 _connected_dof_indices[_j]) /
912 _var.scalingFactor();
913 return (-curr_jac + _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
914 }
915 }
916
917 case ContactFormulation::PENALTY:
918 {
919 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
920 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
921 return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
922 pinfo->_normal(_component) * pinfo->_normal(_component);
923 else
924 return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp];
925 }
926 case ContactFormulation::AUGMENTED_LAGRANGE:
927 {
928 Real normal_comp = _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
929 pinfo->_normal(_component) * pinfo->_normal(_component);
930
931 Real tang_comp = 0.0;
932 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
933 tang_comp = _phi_secondary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
934 (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
935 return normal_comp + tang_comp;
936 }
937
938 case ContactFormulation::TANGENTIAL_PENALTY:
939 {
940 RealVectorValue jac_vec;
941 for (unsigned int i = 0; i < _mesh_dimension; ++i)
942 {
943 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
944 jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
945 _var_objects[i]->scalingFactor();
946 }
947 Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
948 (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
949 pinfo->_normal(_component) * pinfo->_normal(_component);
950
951 Real tang_comp = 0.0;
952 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
953 tang_comp = _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
954 (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
955
956 return normal_comp + tang_comp;
957 }
958
959 default:
960 mooseError("Invalid contact formulation");
961 }
962
963 case ContactModel::GLUED:
964 switch (_formulation)
965 {
966 case ContactFormulation::KINEMATIC:
967 {
968 const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
969 _connected_dof_indices[_j]) /
970 _var.scalingFactor();
971 return (-curr_jac + _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
972 }
973
974 case ContactFormulation::PENALTY:
975 case ContactFormulation::AUGMENTED_LAGRANGE:
976 return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp];
977
978 default:
979 mooseError("Invalid contact formulation");
980 }
981 default:
982 mooseError("Invalid or unavailable contact model");
983 }
984
985 case Moose::SecondaryPrimary:
986 switch (_model)
987 {
988 case ContactModel::FRICTIONLESS:
989 switch (_formulation)
990 {
991 case ContactFormulation::KINEMATIC:
992 {
993 const Node * curr_primary_node = _current_primary->node_ptr(_j);
994
995 RealVectorValue jac_vec;
996 for (unsigned int i = 0; i < _mesh_dimension; ++i)
997 {
998 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
999 jac_vec(i) = (*_jacobian)(dof_number,
1000 curr_primary_node->dof_number(0, _vars[_component], 0)) /
1001 _var_objects[i]->scalingFactor();
1002 }
1003 return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1004 (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1005 pinfo->_normal(_component) * pinfo->_normal(_component);
1006 }
1007
1008 case ContactFormulation::PENALTY:
1009 case ContactFormulation::AUGMENTED_LAGRANGE:
1010 return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1011 pinfo->_normal(_component) * pinfo->_normal(_component);
1012
1013 default:
1014 mooseError("Invalid contact formulation");
1015 }
1016
1017 case ContactModel::COULOMB:
1018 switch (_formulation)
1019 {
1020 case ContactFormulation::KINEMATIC:
1021 {
1022 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1023 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1024 {
1025 const Node * curr_primary_node = _current_primary->node_ptr(_j);
1026
1027 RealVectorValue jac_vec;
1028 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1029 {
1030 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1031 jac_vec(i) =
1032 (*_jacobian)(dof_number,
1033 curr_primary_node->dof_number(0, _vars[_component], 0)) /
1034 _var_objects[i]->scalingFactor();
1035 }
1036 return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1037 (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1038 pinfo->_normal(_component) * pinfo->_normal(_component);
1039 }
1040 else
1041 {
1042 const Node * curr_primary_node = _current_primary->node_ptr(_j);
1043 const Real curr_jac =
1044 (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1045 curr_primary_node->dof_number(0, _vars[_component], 0)) /
1046 _var.scalingFactor();
1047 return (-curr_jac - _phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
1048 }
1049 }
1050
1051 case ContactFormulation::PENALTY:
1052 {
1053 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1054 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1055 return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1056 pinfo->_normal(_component) * pinfo->_normal(_component);
1057 else
1058 return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp];
1059 }
1060 case ContactFormulation::AUGMENTED_LAGRANGE:
1061 {
1062 Real normal_comp = -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1063 pinfo->_normal(_component) * pinfo->_normal(_component);
1064
1065 Real tang_comp = 0.0;
1066 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1067 tang_comp = -_phi_primary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
1068 (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1069 return normal_comp + tang_comp;
1070 }
1071
1072 case ContactFormulation::TANGENTIAL_PENALTY:
1073 {
1074 const Node * curr_primary_node = _current_primary->node_ptr(_j);
1075
1076 RealVectorValue jac_vec;
1077 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1078 {
1079 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1080 jac_vec(i) = (*_jacobian)(dof_number,
1081 curr_primary_node->dof_number(0, _vars[_component], 0)) /
1082 _var_objects[i]->scalingFactor();
1083 }
1084 Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1085 (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1086 pinfo->_normal(_component) * pinfo->_normal(_component);
1087
1088 Real tang_comp = 0.0;
1089 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1090 tang_comp = -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1091 (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1092
1093 return normal_comp + tang_comp;
1094 }
1095
1096 default:
1097 mooseError("Invalid contact formulation");
1098 }
1099 case ContactModel::GLUED:
1100 switch (_formulation)
1101 {
1102 case ContactFormulation::KINEMATIC:
1103 {
1104 const Node * curr_primary_node = _current_primary->node_ptr(_j);
1105 const Real curr_jac =
1106 (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1107 curr_primary_node->dof_number(0, _vars[_component], 0)) /
1108 _var.scalingFactor();
1109 return (-curr_jac - _phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
1110 }
1111
1112 case ContactFormulation::PENALTY:
1113 case ContactFormulation::AUGMENTED_LAGRANGE:
1114 return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp];
1115
1116 default:
1117 mooseError("Invalid contact formulation");
1118 }
1119
1120 default:
1121 mooseError("Invalid or unavailable contact model");
1122 }
1123
1124 case Moose::PrimarySecondary:
1125 switch (_model)
1126 {
1127 case ContactModel::FRICTIONLESS:
1128 switch (_formulation)
1129 {
1130 case ContactFormulation::KINEMATIC:
1131 {
1132 RealVectorValue jac_vec;
1133 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1134 {
1135 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1136 jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1137 _var_objects[i]->scalingFactor();
1138 }
1139 return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1140 _test_primary[_i][_qp];
1141 }
1142
1143 case ContactFormulation::PENALTY:
1144 case ContactFormulation::AUGMENTED_LAGRANGE:
1145 return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1146 pinfo->_normal(_component) * pinfo->_normal(_component);
1147
1148 default:
1149 mooseError("Invalid contact formulation");
1150 }
1151 case ContactModel::COULOMB:
1152 switch (_formulation)
1153 {
1154 case ContactFormulation::KINEMATIC:
1155 {
1156 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1157 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1158 {
1159 RealVectorValue jac_vec;
1160 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1161 {
1162 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1163 jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1164 _var_objects[i]->scalingFactor();
1165 }
1166 return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1167 _test_primary[_i][_qp];
1168 }
1169 else
1170 {
1171 const Real secondary_jac =
1172 (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1173 _connected_dof_indices[_j]) /
1174 _var.scalingFactor();
1175 return secondary_jac * _test_primary[_i][_qp];
1176 }
1177 }
1178
1179 case ContactFormulation::PENALTY:
1180 {
1181 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1182 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1183 return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1184 pinfo->_normal(_component) * pinfo->_normal(_component);
1185 else
1186 return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp];
1187 }
1188 case ContactFormulation::AUGMENTED_LAGRANGE:
1189 {
1190 Real normal_comp = -_phi_secondary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1191 pinfo->_normal(_component) * pinfo->_normal(_component);
1192
1193 Real tang_comp = 0.0;
1194 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1195 tang_comp = -_phi_secondary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1196 (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1197 return normal_comp + tang_comp;
1198 }
1199
1200 case ContactFormulation::TANGENTIAL_PENALTY:
1201 {
1202 RealVectorValue jac_vec;
1203 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1204 {
1205 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1206 jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1207 _var_objects[i]->scalingFactor();
1208 }
1209 Real normal_comp =
1210 pinfo->_normal(_component) * (pinfo->_normal * jac_vec) * _test_primary[_i][_qp];
1211
1212 Real tang_comp = 0.0;
1213 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1214 tang_comp = -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1215 (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1216
1217 return normal_comp + tang_comp;
1218 }
1219
1220 default:
1221 mooseError("Invalid contact formulation");
1222 }
1223
1224 case ContactModel::GLUED:
1225 switch (_formulation)
1226 {
1227 case ContactFormulation::KINEMATIC:
1228 {
1229 const Real secondary_jac =
1230 (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1231 _connected_dof_indices[_j]) /
1232 _var.scalingFactor();
1233 return secondary_jac * _test_primary[_i][_qp];
1234 }
1235
1236 case ContactFormulation::PENALTY:
1237 case ContactFormulation::AUGMENTED_LAGRANGE:
1238 return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp];
1239
1240 default:
1241 mooseError("Invalid contact formulation");
1242 }
1243
1244 default:
1245 mooseError("Invalid or unavailable contact model");
1246 }
1247
1248 case Moose::PrimaryPrimary:
1249 switch (_model)
1250 {
1251 case ContactModel::FRICTIONLESS:
1252 switch (_formulation)
1253 {
1254 case ContactFormulation::KINEMATIC:
1255 return 0.0;
1256
1257 case ContactFormulation::PENALTY:
1258 case ContactFormulation::AUGMENTED_LAGRANGE:
1259 return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1260 pinfo->_normal(_component) * pinfo->_normal(_component);
1261
1262 default:
1263 mooseError("Invalid contact formulation");
1264 }
1265
1266 case ContactModel::COULOMB:
1267 case ContactModel::GLUED:
1268 switch (_formulation)
1269 {
1270 case ContactFormulation::KINEMATIC:
1271 return 0.0;
1272
1273 case ContactFormulation::PENALTY:
1274 {
1275 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1276 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1277 return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1278 pinfo->_normal(_component) * pinfo->_normal(_component);
1279 else
1280 return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp];
1281 }
1282
1283 case ContactFormulation::TANGENTIAL_PENALTY:
1284 {
1285 Real tang_comp = 0.0;
1286 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1287 tang_comp = _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1288 (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1289 return tang_comp; // normal component is zero
1290 }
1291
1292 case ContactFormulation::AUGMENTED_LAGRANGE:
1293 {
1294 Real normal_comp = _phi_primary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1295 pinfo->_normal(_component) * pinfo->_normal(_component);
1296
1297 Real tang_comp = 0.0;
1298 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1299 tang_comp = _phi_primary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1300 (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1301 return normal_comp + tang_comp;
1302 }
1303
1304 default:
1305 mooseError("Invalid contact formulation");
1306 }
1307
1308 default:
1309 mooseError("Invalid or unavailable contact model");
1310 }
1311 }
1312
1313 return 0.0;
1314}
1315
1316Real
1317MechanicalContactConstraint::computeQpOffDiagJacobian(Moose::ConstraintJacobianType type,
1318 unsigned int jvar)
1319{
1320 PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
1321
1322 const Real penalty = getPenalty(*_current_node);
1323 const Real penalty_slip = getTangentialPenalty(*_current_node);
1324
1325 unsigned int coupled_component;
1326 Real normal_component_in_coupled_var_dir = 1.0;
1327 if (getCoupledVarComponent(jvar, coupled_component))
1328 normal_component_in_coupled_var_dir = pinfo->_normal(coupled_component);
1329
1330 switch (type)
1331 {
1332 default:
1333 mooseError("Unhandled ConstraintJacobianType");
1334
1335 case Moose::SecondarySecondary:
1336 switch (_model)
1337 {
1338 case ContactModel::FRICTIONLESS:
1339 switch (_formulation)
1340 {
1341 case ContactFormulation::KINEMATIC:
1342 {
1343 RealVectorValue jac_vec;
1344 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1345 {
1346 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1347 jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1348 _var_objects[i]->scalingFactor();
1349 }
1350 return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1351 (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1352 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1353 }
1354
1355 case ContactFormulation::PENALTY:
1356 case ContactFormulation::AUGMENTED_LAGRANGE:
1357 return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1358 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1359
1360 default:
1361 mooseError("Invalid contact formulation");
1362 }
1363
1364 case ContactModel::COULOMB:
1365 {
1366 if ((_formulation == ContactFormulation::KINEMATIC ||
1367 _formulation == ContactFormulation::TANGENTIAL_PENALTY) &&
1368 (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1369 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1370 {
1371 RealVectorValue jac_vec;
1372 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1373 {
1374 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1375 jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1376 _var_objects[i]->scalingFactor();
1377 }
1378
1379 return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1380 (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1381 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1382 }
1383 else if ((_formulation == ContactFormulation::PENALTY) &&
1384 (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1385 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1386 return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1387 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1388 else if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
1389 {
1390 Real normal_comp = _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1391 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1392
1393 Real tang_comp = 0.0;
1394 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1395 tang_comp = _phi_secondary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
1396 (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1397 return normal_comp + tang_comp;
1398 }
1399 else
1400 {
1401 const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1402 _connected_dof_indices[_j]);
1403 return -curr_jac;
1404 }
1405 }
1406
1407 case ContactModel::GLUED:
1408 {
1409 const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1410 _connected_dof_indices[_j]);
1411 return -curr_jac;
1412 }
1413 default:
1414 mooseError("Invalid or unavailable contact model");
1415 }
1416
1417 case Moose::SecondaryPrimary:
1418 switch (_model)
1419 {
1420 case ContactModel::FRICTIONLESS:
1421 switch (_formulation)
1422 {
1423 case ContactFormulation::KINEMATIC:
1424 {
1425 const Node * curr_primary_node = _current_primary->node_ptr(_j);
1426
1427 RealVectorValue jac_vec;
1428 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1429 {
1430 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1431 jac_vec(i) = (*_jacobian)(dof_number,
1432 curr_primary_node->dof_number(0, _vars[_component], 0)) /
1433 _var_objects[i]->scalingFactor();
1434 }
1435 return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1436 (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1437 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1438 }
1439
1440 case ContactFormulation::PENALTY:
1441 case ContactFormulation::AUGMENTED_LAGRANGE:
1442 return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1443 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1444
1445 default:
1446 mooseError("Invalid contact formulation");
1447 }
1448
1449 case ContactModel::COULOMB:
1450 if ((_formulation == ContactFormulation::KINEMATIC ||
1451 _formulation == ContactFormulation::TANGENTIAL_PENALTY) &&
1452 (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1453 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1454 {
1455 const Node * curr_primary_node = _current_primary->node_ptr(_j);
1456
1457 RealVectorValue jac_vec;
1458 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1459 {
1460 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1461 jac_vec(i) =
1462 (*_jacobian)(dof_number, curr_primary_node->dof_number(0, _vars[_component], 0)) /
1463 _var_objects[i]->scalingFactor();
1464 }
1465
1466 return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1467 (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1468 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1469 }
1470 else if ((_formulation == ContactFormulation::PENALTY) &&
1471 (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1472 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1473 return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1474 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1475 else if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
1476 {
1477 Real normal_comp = -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1478 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1479
1480 Real tang_comp = 0.0;
1481 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1482 tang_comp = -_phi_primary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
1483 (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1484 return normal_comp + tang_comp;
1485 }
1486 else
1487 return 0.0;
1488
1489 case ContactModel::GLUED:
1490 return 0;
1491
1492 default:
1493 mooseError("Invalid or unavailable contact model");
1494 }
1495
1496 case Moose::PrimarySecondary:
1497 switch (_model)
1498 {
1499 case ContactModel::FRICTIONLESS:
1500 switch (_formulation)
1501 {
1502 case ContactFormulation::KINEMATIC:
1503 {
1504 RealVectorValue jac_vec;
1505 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1506 {
1507 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1508 jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1509 _var_objects[i]->scalingFactor();
1510 }
1511 return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1512 _test_primary[_i][_qp];
1513 }
1514
1515 case ContactFormulation::PENALTY:
1516 case ContactFormulation::AUGMENTED_LAGRANGE:
1517 return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1518 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1519
1520 default:
1521 mooseError("Invalid contact formulation");
1522 }
1523 case ContactModel::COULOMB:
1524 switch (_formulation)
1525 {
1526 case ContactFormulation::KINEMATIC:
1527 {
1528 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1529 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1530 {
1531 RealVectorValue jac_vec;
1532 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1533 {
1534 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1535 jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1536 _var_objects[i]->scalingFactor();
1537 }
1538 return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1539 _test_primary[_i][_qp];
1540 }
1541 else
1542 {
1543 const Real secondary_jac =
1544 (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1545 _connected_dof_indices[_j]) /
1546 _var.scalingFactor();
1547 return secondary_jac * _test_primary[_i][_qp];
1548 }
1549 }
1550
1551 case ContactFormulation::PENALTY:
1552 {
1553 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1554 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1555 return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1556 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1557 else
1558 return 0.0;
1559 }
1560 case ContactFormulation::AUGMENTED_LAGRANGE:
1561 {
1562 Real normal_comp = -_phi_secondary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1563 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1564
1565 Real tang_comp = 0.0;
1566 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1567 tang_comp = -_phi_secondary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1568 (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1569 return normal_comp + tang_comp;
1570 }
1571 case ContactFormulation::TANGENTIAL_PENALTY:
1572 {
1573 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1574 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1575 {
1576 RealVectorValue jac_vec;
1577 for (unsigned int i = 0; i < _mesh_dimension; ++i)
1578 {
1579 dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1580 jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1581 _var_objects[i]->scalingFactor();
1582 }
1583 return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1584 _test_primary[_i][_qp];
1585 }
1586 else
1587 return 0.0;
1588 }
1589
1590 default:
1591 mooseError("Invalid contact formulation");
1592 }
1593
1594 case ContactModel::GLUED:
1595 switch (_formulation)
1596 {
1597 case ContactFormulation::KINEMATIC:
1598 {
1599 const Real secondary_jac =
1600 (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1601 _connected_dof_indices[_j]) /
1602 _var.scalingFactor();
1603 return secondary_jac * _test_primary[_i][_qp];
1604 }
1605
1606 case ContactFormulation::PENALTY:
1607 case ContactFormulation::AUGMENTED_LAGRANGE:
1608 return 0.0;
1609
1610 default:
1611 mooseError("Invalid contact formulation");
1612 }
1613
1614 default:
1615 mooseError("Invalid or unavailable contact model");
1616 }
1617
1618 case Moose::PrimaryPrimary:
1619 switch (_model)
1620 {
1621 case ContactModel::FRICTIONLESS:
1622 switch (_formulation)
1623 {
1624 case ContactFormulation::KINEMATIC:
1625 return 0.0;
1626
1627 case ContactFormulation::PENALTY:
1628 case ContactFormulation::AUGMENTED_LAGRANGE:
1629 return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1630 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1631
1632 default:
1633 mooseError("Invalid contact formulation");
1634 }
1635
1636 case ContactModel::COULOMB:
1637 {
1638 if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
1639 {
1640 Real normal_comp = _phi_primary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1641 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1642
1643 if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1644 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1645 return normal_comp;
1646 else
1647 return 0.0;
1648 }
1649 else if (_formulation == ContactFormulation::PENALTY &&
1650 (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1651 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1652 return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1653 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1654 else
1655 return 0.0;
1656 }
1657
1658 case ContactModel::GLUED:
1659 if (_formulation == ContactFormulation::PENALTY &&
1660 (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1661 pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1662 return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1663 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1664 else if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
1665 {
1666 Real normal_comp = _phi_primary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1667 pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1668
1669 Real tang_comp = 0.0;
1670 if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1671 tang_comp = _phi_primary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1672 (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1673 return normal_comp + tang_comp;
1674 }
1675 else
1676 return 0.0;
1677
1678 default:
1679 mooseError("Invalid or unavailable contact model");
1680 }
1681 }
1682
1683 return 0.0;
1684}
1685
1686Real
1687MechanicalContactConstraint::gapOffset(const Node & node)
1688{
1689 Real val = 0;
1690
1691 if (_has_secondary_gap_offset)
1692 val += _secondary_gap_offset_var->getNodalValue(node);
1693
1694 if (_has_mapped_primary_gap_offset)
1695 val += _mapped_primary_gap_offset_var->getNodalValue(node);
1696
1697 return val;
1698}
1699
1700Real
1701MechanicalContactConstraint::nodalArea(const Node & node)
1702{
1703 dof_id_type dof = node.dof_number(_aux_system.number(), _nodal_area_var->number(), 0);
1704
1705 Real area = (*_aux_solution)(dof);
1706 if (area == 0.0)
1707 {
1708 if (_t_step > 1)
1709 mooseError("Zero nodal area found");
1710 else
1711 area = 1.0; // Avoid divide by zero during initialization
1712 }
1713
1714 return area;
1715}
1716
1717Real
1718MechanicalContactConstraint::getPenalty(const Node & node)
1719{
1720 Real penalty = _penalty;
1721 if (_normalize_penalty)
1722 penalty *= nodalArea(node);
1723 return penalty * MathUtils::pow(_penalty_multiplier, _lagrangian_iteration_number);
1724}
1725
1726Real
1727MechanicalContactConstraint::getTangentialPenalty(const Node & node)
1728{
1729 Real penalty = _penalty_tangential;
1730 if (_normalize_penalty)
1731 penalty *= nodalArea(node);
1732
1733 return penalty * MathUtils::pow(_penalty_multiplier, _lagrangian_iteration_number);
1734}
1735
1736void
1737MechanicalContactConstraint::computeJacobian()
1738{
1739 getConnectedDofIndices(_var.number());
1740
1741 prepareMatrixTagNeighbor(
1742 _assembly, _primary_var.number(), _var.number(), Moose::NeighborNeighbor, _Knn);
1743
1744 _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
1745
1746 for (_i = 0; _i < _test_secondary.size(); _i++)
1747 // Loop over the connected dof indices so we can get all the jacobian contributions
1748 for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1749 _Kee(_i, _j) += computeQpJacobian(Moose::SecondarySecondary);
1750
1751 if (_primary_secondary_jacobian)
1752 {
1753 prepareMatrixTagNeighbor(_assembly, _var.number(), _var.number(), Moose::ElementNeighbor, _Ken);
1754 if (_Ken.m() && _Ken.n())
1755 {
1756 for (_i = 0; _i < _test_secondary.size(); _i++)
1757 for (_j = 0; _j < _phi_primary.size(); _j++)
1758 _Ken(_i, _j) += computeQpJacobian(Moose::SecondaryPrimary);
1759 accumulateTaggedLocalMatrix(
1760 _assembly, _var.number(), _var.number(), Moose::ElementNeighbor, _Ken);
1761 }
1762
1763 _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
1764 for (_i = 0; _i < _test_primary.size(); _i++)
1765 // Loop over the connected dof indices so we can get all the jacobian contributions
1766 for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1767 _Kne(_i, _j) += computeQpJacobian(Moose::PrimarySecondary);
1768 }
1769
1770 if (_Knn.m() && _Knn.n())
1771 {
1772 for (_i = 0; _i < _test_primary.size(); _i++)
1773 for (_j = 0; _j < _phi_primary.size(); _j++)
1774 _Knn(_i, _j) += computeQpJacobian(Moose::PrimaryPrimary);
1775 accumulateTaggedLocalMatrix(
1776 _assembly, _primary_var.number(), _var.number(), Moose::NeighborNeighbor, _Knn);
1777 }
1778}
1779
1780void
1781MechanicalContactConstraint::computeOffDiagJacobian(const unsigned int jvar_num)
1782{
1783 getConnectedDofIndices(jvar_num);
1784
1785 _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
1786
1787 prepareMatrixTagNeighbor(
1788 _assembly, _primary_var.number(), jvar_num, Moose::NeighborNeighbor, _Knn);
1789
1790 for (_i = 0; _i < _test_secondary.size(); _i++)
1791 // Loop over the connected dof indices so we can get all the jacobian contributions
1792 for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1793 _Kee(_i, _j) += computeQpOffDiagJacobian(Moose::SecondarySecondary, jvar_num);
1794
1795 if (_primary_secondary_jacobian)
1796 {
1797 prepareMatrixTagNeighbor(_assembly, _var.number(), jvar_num, Moose::ElementNeighbor, _Ken);
1798 for (_i = 0; _i < _test_secondary.size(); _i++)
1799 for (_j = 0; _j < _phi_primary.size(); _j++)
1800 _Ken(_i, _j) += computeQpOffDiagJacobian(Moose::SecondaryPrimary, jvar_num);
1801 accumulateTaggedLocalMatrix(_assembly, _var.number(), jvar_num, Moose::ElementNeighbor, _Ken);
1802
1803 _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
1804 if (_Kne.m() && _Kne.n())
1805 for (_i = 0; _i < _test_primary.size(); _i++)
1806 // Loop over the connected dof indices so we can get all the jacobian contributions
1807 for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1808 _Kne(_i, _j) += computeQpOffDiagJacobian(Moose::PrimarySecondary, jvar_num);
1809 }
1810
1811 for (_i = 0; _i < _test_primary.size(); _i++)
1812 for (_j = 0; _j < _phi_primary.size(); _j++)
1813 _Knn(_i, _j) += computeQpOffDiagJacobian(Moose::PrimaryPrimary, jvar_num);
1814 accumulateTaggedLocalMatrix(
1815 _assembly, _primary_var.number(), jvar_num, Moose::NeighborNeighbor, _Knn);
1816}
1817
1818void
1819MechanicalContactConstraint::getConnectedDofIndices(unsigned int var_num)
1820{
1821 unsigned int component;
1822 if (getCoupledVarComponent(var_num, component) || _non_displacement_vars_jacobian)
1823 {
1824 if (_primary_secondary_jacobian && _connected_secondary_nodes_jacobian)
1825 NodeFaceConstraint::getConnectedDofIndices(var_num);
1826 else
1827 {
1828 _connected_dof_indices.clear();
1829 MooseVariableFEBase & var = _sys.getVariable(0, var_num);
1830 _connected_dof_indices.push_back(var.nodalDofIndex());
1831 }
1832 }
1833
1834 _phi_secondary.resize(_connected_dof_indices.size());
1835
1836 dof_id_type current_node_var_dof_index = _sys.getVariable(0, var_num).nodalDofIndex();
1837 _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 for (unsigned int j = 0; j < _connected_dof_indices.size(); j++)
1843 {
1844 _phi_secondary[j].resize(1);
1845
1846 if (_connected_dof_indices[j] == current_node_var_dof_index)
1847 _phi_secondary[j][_qp] = 1.0;
1848 else
1849 _phi_secondary[j][_qp] = 0.0;
1850 }
1851}
1852
1853bool
1854MechanicalContactConstraint::getCoupledVarComponent(unsigned int var_num, unsigned int & component)
1855{
1856 component = std::numeric_limits<unsigned int>::max();
1857 bool coupled_var_is_disp_var = false;
1858 for (const auto i : make_range(Moose::dim))
1859 {
1860 if (var_num == _vars[i])
1861 {
1862 coupled_var_is_disp_var = true;
1863 component = i;
1864 break;
1865 }
1866 }
1867
1868 return coupled_var_is_disp_var;
1869}
1870
1871void
1872MechanicalContactConstraint::residualEnd()
1873{
1874 if (_component == 0 && (_print_contact_nodes || _contact_linesearch))
1875 {
1876 _communicator.set_union(_current_contact_state);
1877 if (_print_contact_nodes)
1878 {
1879 if (_current_contact_state == _old_contact_state)
1880 _console << "Unchanged contact state. " << _current_contact_state.size()
1881 << " nodes in contact.\n";
1882 else
1883 _console << "Changed contact state!!! " << _current_contact_state.size()
1884 << " nodes in contact.\n";
1885 }
1886 if (_contact_linesearch)
1887 _contact_linesearch->insertSet(_current_contact_state);
1888 _old_contact_state.swap(_current_contact_state);
1889 _current_contact_state.clear();
1890 }
1891}
registerMooseObject
registerMooseObject("ContactApp", MechanicalContactConstraint)
int
void ErrorVector unsigned int
AugmentedLagrangianContactProblemInterface
Class to provide an interface for parameters and routines required to check convergence for the augme...
Definition AugmentedLagrangianContactProblemInterface.h:19
NodeFaceConstraint::_console
const ConsoleStream _console
NodeFaceConstraint::_qp
unsigned int _qp
NodeFaceConstraint::_j
unsigned int _j
NodeFaceConstraint::_i
unsigned int _i
ContactAction::commonParameters
static InputParameters commonParameters()
Define parameters used by multiple contact objects.
Definition ContactAction.C:1829
ContactLineSearchBase
This class implements a custom line search for use with mechanical contact.
Definition ContactLineSearchBase.h:35
ContactLineSearchBase::reset
virtual void reset()
Reset the line search data.
Definition ContactLineSearchBase.C:63
ContactLineSearchBase::insertSet
void insertSet(const std::set< dof_id_type > &mech_set)
Unionize sets from different constraints.
Definition ContactLineSearchBase.C:53
NodeFaceConstraint::getVar
MooseVariable * getVar(const std::string &var_name, unsigned int comp)
NodeFaceConstraint::coupledComponents
unsigned int coupledComponents(const std::string &var_name) const
NodeFaceConstraint::coupled
virtual unsigned int coupled(const std::string &var_name, unsigned int comp=0) const
Executioner::lastSolveConverged
virtual bool lastSolveConverged() const=0
InputParameters::addRequiredCoupledVar
void addRequiredCoupledVar(const std::string &name, const std::string &doc_string)
InputParameters::addRequiredParam
void addRequiredParam(const std::string &name, const std::string &doc_string)
InputParameters::addParam
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
InputParameters::get
std::vector< std::pair< R1, R2 > > get(const std::string &param1, const std::string &param2) const
InputParameters::addClassDescription
void addClassDescription(const std::string &doc_string)
InputParameters::set
T & set(const std::string &name, bool quiet_mode=false)
InputParameters::addCoupledVar
void addCoupledVar(const std::string &name, const std::string &doc_string)
InputParameters::isParamValid
bool isParamValid(const std::string &name) const
MechanicalContactConstraint
A MechanicalContactConstraint forces the value of a variable to be the same on both sides of an inter...
Definition MechanicalContactConstraint.h:28
MechanicalContactConstraint::_connected_secondary_nodes_jacobian
const bool _connected_secondary_nodes_jacobian
Whether to include coupling terms with the nodes connected to the secondary nodes in the Jacobian.
Definition MechanicalContactConstraint.h:132
MechanicalContactConstraint::computeQpResidual
virtual Real computeQpResidual(Moose::ConstraintType type) override
Definition MechanicalContactConstraint.C:794
MechanicalContactConstraint::_al_frictional_force_tolerance
Real _al_frictional_force_tolerance
The tolerance of the frictional force for augmented Lagrangian method.
Definition MechanicalContactConstraint.h:141
MechanicalContactConstraint::updateContactStatefulData
virtual void updateContactStatefulData(bool beginning_of_step)
Definition MechanicalContactConstraint.C:428
MechanicalContactConstraint::_contact_set_mutex
static Threads::spin_mutex _contact_set_mutex
Definition MechanicalContactConstraint.h:148
MechanicalContactConstraint::_augmented_lagrange_problem
AugmentedLagrangianContactProblemInterface *const _augmented_lagrange_problem
Definition MechanicalContactConstraint.h:150
MechanicalContactConstraint::_no_iterations
static const unsigned int _no_iterations
Definition MechanicalContactConstraint.h:151
MechanicalContactConstraint::getCoupledVarComponent
bool getCoupledVarComponent(unsigned int var_num, unsigned int &component)
Determine whether the coupled variable is one of the displacement variables, and find its component.
Definition MechanicalContactConstraint.C:1854
MechanicalContactConstraint::_mapped_primary_gap_offset_var
const MooseVariable *const _mapped_primary_gap_offset_var
Definition MechanicalContactConstraint.h:123
MechanicalContactConstraint::MechanicalContactConstraint
MechanicalContactConstraint(const InputParameters &parameters)
Definition MechanicalContactConstraint.C:123
MechanicalContactConstraint::_contact_linesearch
ContactLineSearchBase * _contact_linesearch
Definition MechanicalContactConstraint.h:143
MechanicalContactConstraint::updateAugmentedLagrangianMultiplier
virtual void updateAugmentedLagrangianMultiplier(bool beginning_of_step)
Definition MechanicalContactConstraint.C:256
MechanicalContactConstraint::_non_displacement_vars_jacobian
const bool _non_displacement_vars_jacobian
Whether to include coupling terms with non-displacement variables in the Jacobian.
Definition MechanicalContactConstraint.h:134
MechanicalContactConstraint::_var_objects
std::vector< MooseVariable * > _var_objects
Definition MechanicalContactConstraint.h:117
MechanicalContactConstraint::_primary_secondary_jacobian
const bool _primary_secondary_jacobian
Whether to include coupling between the primary and secondary nodes in the Jacobian.
Definition MechanicalContactConstraint.h:130
MechanicalContactConstraint::_al_incremental_slip_tolerance
Real _al_incremental_slip_tolerance
The tolerance of the incremental slip for augmented Lagrangian method.
Definition MechanicalContactConstraint.h:139
MechanicalContactConstraint::computeQpSecondaryValue
virtual Real computeQpSecondaryValue() override
Definition MechanicalContactConstraint.C:788
MechanicalContactConstraint::computeContactForce
void computeContactForce(const Node &node, PenetrationInfo *pinfo, bool update_contact_set)
Definition MechanicalContactConstraint.C:502
MechanicalContactConstraint::computeJacobian
virtual void computeJacobian() override
Computes the jacobian for the current element.
Definition MechanicalContactConstraint.C:1737
MechanicalContactConstraint::getConnectedDofIndices
virtual void getConnectedDofIndices(unsigned int var_num) override
Get the dof indices of the nodes connected to the secondary node for a specific variable.
Definition MechanicalContactConstraint.C:1819
MechanicalContactConstraint::_secondary_gap_offset_var
const MooseVariable *const _secondary_gap_offset_var
Definition MechanicalContactConstraint.h:121
MechanicalContactConstraint::computeQpJacobian
virtual Real computeQpJacobian(Moose::ConstraintJacobianType type) override
Definition MechanicalContactConstraint.C:847
MechanicalContactConstraint::_lagrangian_iteration_number
const unsigned int & _lagrangian_iteration_number
Definition MechanicalContactConstraint.h:152
MechanicalContactConstraint::_residual_copy
NumericVector< Number > & _residual_copy
Definition MechanicalContactConstraint.h:111
MechanicalContactConstraint::_has_secondary_gap_offset
const bool _has_secondary_gap_offset
gap offset from either secondary, primary or both
Definition MechanicalContactConstraint.h:120
MechanicalContactConstraint::computeOffDiagJacobian
virtual void computeOffDiagJacobian(unsigned int jvar) override
Compute off-diagonal Jacobian entries.
Definition MechanicalContactConstraint.C:1781
MechanicalContactConstraint::computeQpOffDiagJacobian
virtual Real computeQpOffDiagJacobian(Moose::ConstraintJacobianType type, unsigned int jvar) override
Compute off-diagonal Jacobian entries.
Definition MechanicalContactConstraint.C:1317
MechanicalContactConstraint::_al_penetration_tolerance
Real _al_penetration_tolerance
The tolerance of the penetration for augmented Lagrangian method.
Definition MechanicalContactConstraint.h:137
MechanicalContactConstraint::_current_contact_state
std::set< dof_id_type > _current_contact_state
Definition MechanicalContactConstraint.h:144
MooseApp::getExecutioner
Executioner * getExecutioner() const
NodeFaceConstraint::parameters
const InputParameters & parameters() const
NodeFaceConstraint::type
const std::string & type() const
NodeFaceConstraint::mooseError
void mooseError(Args &&... args) const
NodeFaceConstraint::isParamValid
bool isParamValid(const std::string &name) const
MooseMesh::nodeRef
virtual const Node & nodeRef(const dof_id_type i) const
MooseVariableFE::scalingFactor
void scalingFactor(const std::vector< Real > &factor)
MooseVariableFE::number
unsigned int number() const
MooseVariableFE::getNodalValue
DofValue getNodalValue(const Node &node) const
MooseVariableFieldBase::nodalDofIndex
virtual const dof_id_type & nodalDofIndex() const=0
NodeFaceConstraint::_penetration_locator
PenetrationLocator & _penetration_locator
NodeFaceConstraint::_current_node
const Node *const & _current_node
NodeFaceConstraint::_Kee
DenseMatrix< Number > _Kee
NodeFaceConstraint::_test_secondary
VariableTestValue _test_secondary
NodeFaceConstraint::_current_primary
const Elem *const & _current_primary
NodeFaceConstraint::_connected_dof_indices
std::vector< dof_id_type > _connected_dof_indices
NodeFaceConstraint::_test_primary
const VariableTestValue & _test_primary
NodeFaceConstraint::_overwrite_secondary_residual
bool _overwrite_secondary_residual
NodeFaceConstraint::_Kne
DenseMatrix< Number > _Kne
NodeFaceConstraint::validParams
static InputParameters validParams()
NodeFaceConstraint::_var
MooseVariable & _var
NodeFaceConstraint::_u_secondary
const VariableValue & _u_secondary
NodeFaceConstraint::_phi_primary
const VariablePhiValue & _phi_primary
NodeFaceConstraint::_phi_secondary
VariablePhiValue _phi_secondary
NodeFaceConstraint::_primary_var
MooseVariable & _primary_var
NodeFaceConstraint::getConnectedDofIndices
virtual void getConnectedDofIndices(unsigned int var_num)
PenetrationInfo::_mech_status
MECH_STATUS_ENUM _mech_status
PenetrationInfo::_contact_force_old
RealVectorValue _contact_force_old
PenetrationInfo::_locked_this_step
unsigned int _locked_this_step
PenetrationInfo::_contact_force
RealVectorValue _contact_force
PenetrationInfo::_incremental_slip
Point _incremental_slip
PenetrationInfo::isCaptured
bool isCaptured() const
PenetrationInfo::_lagrange_multiplier
Real _lagrange_multiplier
PenetrationInfo::_lagrange_multiplier_slip
RealVectorValue _lagrange_multiplier_slip
PenetrationInfo::_normal
RealVectorValue _normal
PenetrationInfo::_stick_locked_this_step
unsigned int _stick_locked_this_step
PenetrationLocator::setNormalSmoothingMethod
void setNormalSmoothingMethod(std::string nsmString)
PenetrationLocator::setUpdate
void setUpdate(bool update)
PenetrationLocator::setTangentialTolerance
void setTangentialTolerance(Real tangential_tolerance)
PenetrationLocator::setNormalSmoothingDistance
void setNormalSmoothingDistance(Real normal_smoothing_distance)
PenetrationLocator::_penetration_info
std::map< dof_id_type, PenetrationInfo * > & _penetration_info
NodeFaceConstraint::_mesh
MooseMesh & _mesh
NodeFaceConstraint::_assembly
Assembly & _assembly
NodeFaceConstraint::_subproblem
SubProblem & _subproblem
NodeFaceConstraint::_sys
SystemBase & _sys
SubProblem::computingNonlinearResid
bool computingNonlinearResid() const
SystemBase::getVariable
MooseVariableFieldBase & getVariable(THREAD_ID tid, const std::string &var_name) const
SystemBase::number
unsigned int number() const
libMesh::Parallel::Communicator::max
void max(const T &r, T &o, Request &req) const
libMesh::Parallel::Communicator::set_union
void set_union(T &data, const unsigned int root_id) const
NodeFaceConstraint::accumulateTaggedLocalMatrix
void accumulateTaggedLocalMatrix()
NodeFaceConstraint::prepareMatrixTagNeighbor
void prepareMatrixTagNeighbor(Assembly &assembly, unsigned int ivar, unsigned int jvar, Moose::DGJacobianType type)
libMesh::DenseMatrix::n
unsigned int n() const
libMesh::DenseMatrix::m
unsigned int m() const
libMesh::DenseMatrix::resize
void resize(const unsigned int new_m, const unsigned int new_n)
libMesh::ParallelObject::_communicator
const Parallel::Communicator & _communicator
MathUtils::pow
T pow(T x, int e)
Moose::ConstraintType
ConstraintType
Moose::ConstraintJacobianType
ConstraintJacobianType
Moose::SecondarySecondary
SecondarySecondary
Moose::SecondaryPrimary
SecondaryPrimary
Moose::PrimarySecondary
PrimarySecondary
Moose::PrimaryPrimary
PrimaryPrimary
Moose::dim
static constexpr std::size_t dim
Moose::NeighborNeighbor
NeighborNeighbor
Moose::ElementNeighbor
ElementNeighbor
libMesh
The following methods are specializations for using the Parallel::packed_range_* routines for a vecto...
distance
Real distance(const Point &p)
Generated on Fri Aug 21 2026 13:40:58 for https://mooseframework.inl.gov by doxygen 1.9.8