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 3948 : MechanicalContactConstraint::validParams()
36 : {
37 3948 : InputParameters params = NodeFaceConstraint::validParams();
38 3948 : params += ContactAction::commonParameters();
39 :
40 7896 : params.addRequiredParam<BoundaryName>("boundary", "The primary boundary");
41 7896 : params.addParam<BoundaryName>("secondary", "The secondary boundary");
42 7896 : 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 7896 : params.addCoupledVar(
48 : "displacements",
49 : "The displacements appropriate for the simulation geometry and coordinate system");
50 :
51 7896 : params.addCoupledVar("secondary_gap_offset", "offset to the gap distance from secondary side");
52 7896 : params.addCoupledVar("mapped_primary_gap_offset",
53 : "offset to the gap distance mapped from primary side");
54 7896 : params.addRequiredCoupledVar("nodal_area", "The nodal area");
55 :
56 3948 : params.set<bool>("use_displaced_mesh") = true;
57 7896 : params.addParam<Real>(
58 : "penalty",
59 7896 : 1e8,
60 : "The penalty to apply. This can vary depending on the stiffness of your materials");
61 7896 : params.addParam<Real>("penalty_multiplier",
62 7896 : 1.0,
63 : "The growth factor for the penalty applied at the end of each augmented "
64 : "Lagrange update iteration");
65 7896 : params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
66 7896 : params.addParam<Real>("tangential_tolerance",
67 : "Tangential distance to extend edges of contact surfaces");
68 7896 : params.addParam<Real>(
69 7896 : "capture_tolerance", 0, "Normal distance from surface within which nodes are captured");
70 :
71 7896 : params.addParam<Real>("tension_release",
72 7896 : 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 7896 : params.addParam<bool>(
78 : "normalize_penalty",
79 7896 : false,
80 : "Whether to normalize the penalty parameter with the nodal area for penalty contact.");
81 7896 : params.addParam<bool>(
82 : "primary_secondary_jacobian",
83 7896 : true,
84 : "Whether to include jacobian entries coupling primary and secondary nodes.");
85 7896 : params.addParam<bool>(
86 : "connected_secondary_nodes_jacobian",
87 7896 : true,
88 : "Whether to include jacobian entries coupling nodes connected to secondary nodes.");
89 7896 : params.addParam<bool>("non_displacement_variables_jacobian",
90 7896 : true,
91 : "Whether to include jacobian entries coupling with variables that are not "
92 : "displacement variables.");
93 7896 : params.addParam<unsigned int>("stick_lock_iterations",
94 7896 : 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 7896 : params.addParam<Real>("stick_unlock_factor",
98 7896 : 1.5,
99 : "Factor by which frictional capacity must be "
100 : "exceeded to permit stick-locked node to slip "
101 : "again.");
102 7896 : params.addParam<Real>("al_penetration_tolerance",
103 : "The tolerance of the penetration for augmented Lagrangian method.");
104 7896 : params.addParam<Real>("al_incremental_slip_tolerance",
105 : "The tolerance of the incremental slip for augmented Lagrangian method.");
106 :
107 7896 : params.addParam<Real>("al_frictional_force_tolerance",
108 : "The tolerance of the frictional force for augmented Lagrangian method.");
109 7896 : params.addParam<bool>(
110 7896 : "print_contact_nodes", false, "Whether to print the number of nodes in contact.");
111 :
112 3948 : 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 3948 : return params;
119 0 : }
120 :
121 : Threads::spin_mutex MechanicalContactConstraint::_contact_set_mutex;
122 :
123 2791 : MechanicalContactConstraint::MechanicalContactConstraint(const InputParameters & parameters)
124 : : NodeFaceConstraint(parameters),
125 2791 : _displaced_problem(parameters.get<FEProblemBase *>("_fe_problem_base")->getDisplacedProblem()),
126 5582 : _component(getParam<unsigned int>("component")),
127 5582 : _model(getParam<MooseEnum>("model").getEnum<ContactModel>()),
128 5582 : _formulation(getParam<MooseEnum>("formulation").getEnum<ContactFormulation>()),
129 5582 : _normalize_penalty(getParam<bool>("normalize_penalty")),
130 5582 : _penalty(getParam<Real>("penalty")),
131 5582 : _penalty_multiplier(getParam<Real>("penalty_multiplier")),
132 5582 : _friction_coefficient(getParam<Real>("friction_coefficient")),
133 5582 : _tension_release(getParam<Real>("tension_release")),
134 5582 : _capture_tolerance(getParam<Real>("capture_tolerance")),
135 5582 : _stick_lock_iterations(getParam<unsigned int>("stick_lock_iterations")),
136 5582 : _stick_unlock_factor(getParam<Real>("stick_unlock_factor")),
137 2791 : _update_stateful_data(true),
138 2791 : _residual_copy(_sys.residualGhosted()),
139 2791 : _mesh_dimension(_mesh.dimension()),
140 2791 : _vars(3, libMesh::invalid_uint),
141 2791 : _var_objects(3, nullptr),
142 2791 : _has_secondary_gap_offset(isCoupled("secondary_gap_offset")),
143 2791 : _secondary_gap_offset_var(_has_secondary_gap_offset ? getVar("secondary_gap_offset", 0)
144 : : nullptr),
145 2791 : _has_mapped_primary_gap_offset(isCoupled("mapped_primary_gap_offset")),
146 2791 : _mapped_primary_gap_offset_var(
147 2791 : _has_mapped_primary_gap_offset ? getVar("mapped_primary_gap_offset", 0) : nullptr),
148 2791 : _nodal_area_var(getVar("nodal_area", 0)),
149 2791 : _aux_system(_nodal_area_var->sys()),
150 2791 : _aux_solution(_aux_system.currentSolution()),
151 5582 : _primary_secondary_jacobian(getParam<bool>("primary_secondary_jacobian")),
152 5582 : _connected_secondary_nodes_jacobian(getParam<bool>("connected_secondary_nodes_jacobian")),
153 5582 : _non_displacement_vars_jacobian(getParam<bool>("non_displacement_variables_jacobian")),
154 2791 : _contact_linesearch(dynamic_cast<ContactLineSearchBase *>(_subproblem.getLineSearch())),
155 5582 : _print_contact_nodes(getParam<bool>("print_contact_nodes")),
156 2791 : _augmented_lagrange_problem(
157 2791 : dynamic_cast<AugmentedLagrangianContactProblemInterface *>(&_fe_problem)),
158 2791 : _lagrangian_iteration_number(_augmented_lagrange_problem
159 2791 : ? _augmented_lagrange_problem->getLagrangianIterationNumber()
160 5582 : : _no_iterations)
161 : {
162 2791 : _overwrite_secondary_residual = false;
163 :
164 5582 : if (isParamValid("displacements"))
165 : {
166 : // modern parameter scheme for displacements
167 18348 : for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
168 : {
169 6383 : _vars[i] = coupled("displacements", i);
170 6383 : _var_objects[i] = getVar("displacements", i);
171 : }
172 : }
173 :
174 5582 : if (parameters.isParamValid("tangential_tolerance"))
175 5892 : _penetration_locator.setTangentialTolerance(getParam<Real>("tangential_tolerance"));
176 :
177 5582 : if (parameters.isParamValid("normal_smoothing_distance"))
178 669 : _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
179 :
180 5582 : if (parameters.isParamValid("normal_smoothing_method"))
181 0 : _penetration_locator.setNormalSmoothingMethod(
182 : parameters.get<std::string>("normal_smoothing_method"));
183 :
184 2791 : 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 2791 : if (_model == ContactModel::GLUED)
188 392 : _penetration_locator.setUpdate(false);
189 :
190 2791 : 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 2791 : _penalty_tangential = _penalty;
195 :
196 2791 : if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
197 : {
198 348 : if (_model == ContactModel::GLUED)
199 0 : mooseError("The Augmented Lagrangian contact formulation does not support GLUED case.");
200 :
201 348 : if (!_augmented_lagrange_problem)
202 0 : mooseError("The Augmented Lagrangian contact formulation must use "
203 : "AugmentedLagrangianContactProblem.");
204 :
205 696 : if (!parameters.isParamValid("al_penetration_tolerance"))
206 0 : mooseError("For Augmented Lagrangian contact, al_penetration_tolerance must be provided.");
207 : else
208 348 : _al_penetration_tolerance = parameters.get<Real>("al_penetration_tolerance");
209 :
210 348 : if (_model != ContactModel::FRICTIONLESS)
211 : {
212 450 : if (!parameters.isParamValid("al_incremental_slip_tolerance") ||
213 450 : !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 150 : _al_incremental_slip_tolerance = parameters.get<Real>("al_incremental_slip_tolerance");
222 150 : _al_frictional_force_tolerance = parameters.get<Real>("al_frictional_force_tolerance");
223 : }
224 : }
225 : }
226 2791 : }
227 :
228 : void
229 12816 : MechanicalContactConstraint::timestepSetup()
230 : {
231 12816 : if (_component == 0)
232 : {
233 6150 : updateContactStatefulData(/* beginning_of_step = */ true);
234 6150 : if (_formulation == ContactFormulation::AUGMENTED_LAGRANGE)
235 151 : updateAugmentedLagrangianMultiplier(/* beginning_of_step = */ true);
236 :
237 6150 : _update_stateful_data = false;
238 :
239 6150 : if (_contact_linesearch)
240 13 : _contact_linesearch->reset();
241 : }
242 12816 : }
243 :
244 : void
245 37090 : MechanicalContactConstraint::jacobianSetup()
246 : {
247 37090 : if (_component == 0)
248 : {
249 17927 : if (_update_stateful_data)
250 12224 : updateContactStatefulData(/* beginning_of_step = */ false);
251 17927 : _update_stateful_data = true;
252 : }
253 37090 : }
254 :
255 : void
256 303 : MechanicalContactConstraint::updateAugmentedLagrangianMultiplier(bool beginning_of_step)
257 : {
258 2230 : for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
259 : {
260 1927 : if (!pinfo)
261 175 : continue;
262 :
263 1752 : const Node & node = _mesh.nodeRef(secondary_node_num);
264 1752 : if (node.n_comp(_sys.number(), _vars[_component]) < 1)
265 0 : continue;
266 :
267 1752 : const Real distance = pinfo->_normal * (pinfo->_closest_point - node) - gapOffset(node);
268 :
269 1752 : if (beginning_of_step && _model == ContactModel::COULOMB)
270 : {
271 280 : pinfo->_lagrange_multiplier_slip.zero();
272 280 : if (pinfo->isCaptured())
273 0 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
274 : }
275 :
276 1752 : if (pinfo->isCaptured())
277 : {
278 1028 : if (_model == ContactModel::FRICTIONLESS)
279 678 : pinfo->_lagrange_multiplier -= getPenalty(node) * distance;
280 :
281 1028 : 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 55 : pinfo->_lagrange_multiplier_slip =
315 : -tau_old + capacity * contact_force_tangential / tan_mag;
316 55 : if (MooseUtils::absoluteFuzzyEqual(capacity, 0.0))
317 55 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
318 : else
319 0 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
320 : }
321 : else
322 : {
323 295 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
324 : pinfo->_lagrange_multiplier_slip += penalty_slip * tangential_inc_slip;
325 : }
326 : }
327 : }
328 : }
329 : }
330 303 : }
331 :
332 : bool
333 500 : MechanicalContactConstraint::AugmentedLagrangianContactConverged()
334 : {
335 : Real contactResidual = 0.0;
336 500 : unsigned int converged = 0;
337 :
338 2862 : for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
339 : {
340 2514 : if (!pinfo)
341 0 : continue;
342 :
343 2514 : const auto & node = _mesh.nodeRef(secondary_node_num);
344 :
345 : // Skip this pinfo if there are no DOFs on this node.
346 2514 : if (node.n_comp(_sys.number(), _vars[_component]) < 1)
347 0 : continue;
348 :
349 2514 : const Real distance = pinfo->_normal * (pinfo->_closest_point - node) - gapOffset(node);
350 :
351 2514 : if (pinfo->isCaptured())
352 : {
353 2514 : if (contactResidual < std::abs(distance))
354 : contactResidual = std::abs(distance);
355 :
356 : // penetration < tol
357 2514 : if (contactResidual > _al_penetration_tolerance)
358 : {
359 152 : converged = 1;
360 152 : break;
361 : }
362 :
363 2362 : 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 500 : _communicator.max(converged);
411 :
412 500 : if (converged == 1)
413 152 : _console << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied "
414 152 : << std::endl;
415 348 : else if (converged == 2)
416 0 : _console << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied "
417 0 : << std::endl;
418 348 : 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 18374 : MechanicalContactConstraint::updateContactStatefulData(bool beginning_of_step)
429 : {
430 160504 : for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
431 : {
432 142130 : if (!pinfo)
433 5447 : continue;
434 :
435 136683 : const Node & node = _mesh.nodeRef(secondary_node_num);
436 136683 : if (node.n_comp(_sys.number(), _vars[_component]) < 1)
437 0 : continue;
438 :
439 136683 : if (beginning_of_step)
440 : {
441 43153 : if (_app.getExecutioner()->lastSolveConverged())
442 : {
443 41291 : pinfo->_contact_force_old = pinfo->_contact_force;
444 41291 : pinfo->_accumulated_slip_old = pinfo->_accumulated_slip;
445 41291 : pinfo->_frictional_energy_old = pinfo->_frictional_energy;
446 41291 : pinfo->_mech_status_old = pinfo->_mech_status;
447 : }
448 1862 : else if (pinfo->_mech_status_old == PenetrationInfo::MS_NO_CONTACT &&
449 1686 : pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
450 : {
451 : // The penetration info object could be based on a bad state so delete it
452 1226 : delete pinfo;
453 1226 : pinfo = NULL;
454 1226 : continue;
455 : }
456 :
457 41927 : pinfo->_locked_this_step = 0;
458 41927 : pinfo->_stick_locked_this_step = 0;
459 41927 : pinfo->_starting_elem = pinfo->_elem;
460 41927 : pinfo->_starting_side_num = pinfo->_side_num;
461 41927 : pinfo->_starting_closest_point_ref = pinfo->_closest_point_ref;
462 : }
463 135457 : pinfo->_incremental_slip_prev_iter = pinfo->_incremental_slip;
464 : }
465 18374 : }
466 :
467 : bool
468 2166427 : MechanicalContactConstraint::shouldApply()
469 : {
470 : bool in_contact = false;
471 :
472 : std::map<dof_id_type, PenetrationInfo *>::iterator found =
473 2166427 : _penetration_locator._penetration_info.find(_current_node->id());
474 2166427 : if (found != _penetration_locator._penetration_info.end())
475 : {
476 2166427 : PenetrationInfo * pinfo = found->second;
477 2166427 : if (pinfo != NULL)
478 : {
479 2166427 : 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 2166427 : if (_component == 0)
484 924473 : computeContactForce(*_current_node, pinfo, is_nonlinear);
485 :
486 2166427 : if (pinfo->isCaptured())
487 : {
488 : in_contact = true;
489 1408388 : if (is_nonlinear)
490 : {
491 : Threads::spin_mutex::scoped_lock lock(_contact_set_mutex);
492 168507 : _current_contact_state.insert(_current_node->id());
493 : }
494 : }
495 : }
496 : }
497 :
498 2166427 : return in_contact;
499 : }
500 :
501 : void
502 924473 : MechanicalContactConstraint::computeContactForce(const Node & node,
503 : PenetrationInfo * pinfo,
504 : bool update_contact_set)
505 : {
506 : // Build up residual vector
507 : RealVectorValue res_vec;
508 3090900 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
509 : {
510 2166427 : dof_id_type dof_number = node.dof_number(0, _vars[i], 0);
511 2166427 : res_vec(i) = _residual_copy(dof_number) / _var_objects[i]->scalingFactor();
512 : }
513 :
514 : RealVectorValue distance_vec(node - pinfo->_closest_point);
515 924473 : if (distance_vec.norm() != 0)
516 888177 : 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 924473 : 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 9504 : if (_formulation == ContactFormulation::KINEMATIC ||
534 : _formulation == ContactFormulation::TANGENTIAL_PENALTY)
535 5134 : ++pinfo->_locked_this_step;
536 : }
537 :
538 924473 : if (!pinfo->isCaptured())
539 340958 : return;
540 :
541 583515 : const Real penalty = getPenalty(node);
542 583515 : const Real penalty_slip = getTangentialPenalty(node);
543 :
544 : RealVectorValue pen_force(penalty * distance_vec);
545 :
546 583515 : switch (_model)
547 : {
548 404344 : case ContactModel::FRICTIONLESS:
549 404344 : switch (_formulation)
550 : {
551 223133 : case ContactFormulation::KINEMATIC:
552 223133 : pinfo->_contact_force = -pinfo->_normal * (pinfo->_normal * res_vec);
553 223133 : break;
554 :
555 162617 : case ContactFormulation::PENALTY:
556 162617 : pinfo->_contact_force = pinfo->_normal * (pinfo->_normal * pen_force);
557 162617 : break;
558 :
559 18594 : case ContactFormulation::AUGMENTED_LAGRANGE:
560 18594 : pinfo->_contact_force =
561 : (pinfo->_normal *
562 : (pinfo->_normal * (pen_force + pinfo->_lagrange_multiplier * pinfo->_normal)));
563 18594 : break;
564 :
565 0 : default:
566 0 : mooseError("Invalid contact formulation");
567 : break;
568 : }
569 404344 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
570 404344 : break;
571 81255 : case ContactModel::COULOMB:
572 81255 : switch (_formulation)
573 : {
574 25819 : case ContactFormulation::KINEMATIC:
575 : {
576 : // Frictional capacity
577 49055 : const Real capacity(_friction_coefficient *
578 25819 : (res_vec * pinfo->_normal > 0 ? res_vec * pinfo->_normal : 0));
579 :
580 : // Normal and tangential components of predictor force
581 25819 : 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 25819 : const Real tan_mag(contact_force_tangential.norm());
592 25819 : const Real tangential_inc_slip_mag = tangential_inc_slip.norm();
593 25819 : const Real slip_tol = capacity / penalty;
594 25819 : pinfo->_slip_tol = slip_tol;
595 :
596 25819 : if ((tangential_inc_slip_mag > slip_tol || tan_mag > capacity) &&
597 24655 : (pinfo->_stick_locked_this_step < _stick_lock_iterations ||
598 0 : tan_mag > capacity * _stick_unlock_factor))
599 : {
600 24655 : 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 24655 : 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 24157 : if (slip_dot_tang_force < capacity)
612 : slipped_too_far = true;
613 : }
614 :
615 : if (slipped_too_far) // slip back along slip increment
616 13557 : pinfo->_contact_force = contact_force_normal + capacity * slip_inc_direction;
617 : else
618 : {
619 11098 : if (tan_mag > 0) // slip along tangential force direction
620 11090 : pinfo->_contact_force =
621 : contact_force_normal + capacity * contact_force_tangential / tan_mag;
622 : else // treat as frictionless
623 8 : pinfo->_contact_force = contact_force_normal;
624 : }
625 24655 : if (capacity == 0)
626 2602 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
627 : else
628 22053 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
629 : }
630 : else
631 : {
632 1164 : if (pinfo->_mech_status != PenetrationInfo::MS_STICKING &&
633 : pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
634 304 : ++pinfo->_stick_locked_this_step;
635 1164 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
636 : }
637 : break;
638 : }
639 :
640 12475 : case ContactFormulation::PENALTY:
641 : {
642 12475 : distance_vec =
643 : pinfo->_incremental_slip +
644 12475 : (pinfo->_normal * (node - pinfo->_closest_point) + gapOffset(node)) * pinfo->_normal;
645 : pen_force = penalty * distance_vec;
646 :
647 : // Frictional capacity
648 21045 : const Real capacity(_friction_coefficient *
649 12475 : (pen_force * pinfo->_normal < 0 ? -pen_force * pinfo->_normal : 0));
650 :
651 : // Elastic predictor
652 12475 : 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 12475 : const Real tan_mag(contact_force_tangential.norm());
661 :
662 12475 : if (tan_mag > capacity)
663 : {
664 2251 : pinfo->_contact_force =
665 : contact_force_normal + capacity * contact_force_tangential / tan_mag;
666 2251 : if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
667 1744 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
668 : else
669 507 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
670 : }
671 : else
672 10224 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
673 : break;
674 : }
675 :
676 12774 : case ContactFormulation::AUGMENTED_LAGRANGE:
677 : {
678 12774 : distance_vec =
679 12774 : (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 12774 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
696 5472 : pinfo->_contact_force =
697 : contact_force_normal + contact_force_tangential + inc_pen_force_tangential;
698 : else
699 7302 : pinfo->_contact_force = contact_force_normal + contact_force_tangential;
700 :
701 : break;
702 : }
703 :
704 30187 : case ContactFormulation::TANGENTIAL_PENALTY:
705 : {
706 : // Frictional capacity (kinematic formulation)
707 32846 : 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 30187 : const Real tan_mag(contact_force_tangential.norm());
721 :
722 30187 : if (tan_mag > capacity)
723 : {
724 24226 : pinfo->_contact_force =
725 : contact_force_normal + capacity * contact_force_tangential / tan_mag;
726 24226 : if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
727 3597 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
728 : else
729 20629 : pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
730 : }
731 : else
732 : {
733 5961 : pinfo->_contact_force = contact_force_normal + contact_force_tangential;
734 5961 : 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 97916 : case ContactModel::GLUED:
746 97916 : switch (_formulation)
747 : {
748 : case ContactFormulation::KINEMATIC:
749 66529 : pinfo->_contact_force = -res_vec;
750 66529 : break;
751 :
752 31387 : case ContactFormulation::PENALTY:
753 31387 : pinfo->_contact_force = pen_force;
754 31387 : 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 97916 : pinfo->_mech_status = PenetrationInfo::MS_STICKING;
766 97916 : break;
767 :
768 0 : default:
769 0 : mooseError("Invalid or unavailable contact model");
770 : break;
771 : }
772 :
773 : // Release
774 72367 : if (update_contact_set && _model != ContactModel::GLUED && pinfo->isCaptured() &&
775 637566 : !newly_captured && _tension_release >= 0.0 &&
776 52285 : (_contact_linesearch ? true : pinfo->_locked_this_step < 2))
777 : {
778 52279 : const Real contact_pressure = -(pinfo->_normal * pinfo->_contact_force) / nodalArea(node);
779 52279 : if (-contact_pressure >= _tension_release)
780 : {
781 : pinfo->release();
782 : pinfo->_contact_force.zero();
783 : }
784 : }
785 : }
786 :
787 : Real
788 20023 : MechanicalContactConstraint::computeQpSecondaryValue()
789 : {
790 20023 : return _u_secondary[_qp];
791 : }
792 :
793 : Real
794 13725740 : MechanicalContactConstraint::computeQpResidual(Moose::ConstraintType type)
795 : {
796 13725740 : PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
797 13725740 : Real resid = pinfo->_contact_force(_component);
798 13725740 : switch (type)
799 : {
800 1261451 : case Moose::Secondary:
801 1261451 : if (_formulation == ContactFormulation::KINEMATIC)
802 : {
803 717582 : RealVectorValue distance_vec(*_current_node - pinfo->_closest_point);
804 717582 : if (distance_vec.norm() != 0)
805 680486 : distance_vec += gapOffset(*_current_node) * pinfo->_normal * distance_vec.unit() *
806 1360972 : distance_vec.unit();
807 :
808 717582 : const Real penalty = getPenalty(*_current_node);
809 : RealVectorValue pen_force(penalty * distance_vec);
810 :
811 717582 : if (_model == ContactModel::FRICTIONLESS)
812 511289 : resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
813 206293 : else if (_model == ContactModel::COULOMB)
814 : {
815 : distance_vec = distance_vec - pinfo->_incremental_slip;
816 42754 : pen_force = penalty * distance_vec;
817 42754 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
818 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
819 40502 : resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
820 : else
821 2252 : resid += pen_force(_component);
822 : }
823 163539 : else if (_model == ContactModel::GLUED)
824 163539 : resid += pen_force(_component);
825 : }
826 543869 : else if (_formulation == ContactFormulation::TANGENTIAL_PENALTY &&
827 51424 : _model == ContactModel::COULOMB)
828 : {
829 51424 : RealVectorValue distance_vec = (pinfo->_normal * (*_current_node - pinfo->_closest_point) +
830 51424 : gapOffset(*_current_node)) *
831 : pinfo->_normal;
832 :
833 51424 : const Real penalty = getPenalty(*_current_node);
834 : RealVectorValue pen_force(penalty * distance_vec);
835 51424 : resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
836 : }
837 1261451 : return _test_secondary[_i][_qp] * resid;
838 :
839 12464289 : case Moose::Primary:
840 12464289 : return _test_primary[_i][_qp] * -resid;
841 : }
842 :
843 : return 0.0;
844 : }
845 :
846 : Real
847 50271509 : MechanicalContactConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
848 : {
849 50271509 : PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
850 :
851 50271509 : const Real penalty = getPenalty(*_current_node);
852 50271509 : const Real penalty_slip = getTangentialPenalty(*_current_node);
853 :
854 50271509 : switch (type)
855 : {
856 0 : default:
857 0 : mooseError("Unhandled ConstraintJacobianType");
858 :
859 1735693 : case Moose::SecondarySecondary:
860 1735693 : switch (_model)
861 : {
862 1358576 : case ContactModel::FRICTIONLESS:
863 1358576 : switch (_formulation)
864 : {
865 : case ContactFormulation::KINEMATIC:
866 : {
867 : RealVectorValue jac_vec;
868 3621180 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
869 : {
870 2689818 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
871 2689818 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
872 2689818 : _var_objects[i]->scalingFactor();
873 : }
874 931362 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
875 931362 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
876 931362 : pinfo->_normal(_component) * pinfo->_normal(_component);
877 : }
878 :
879 427214 : case ContactFormulation::PENALTY:
880 : case ContactFormulation::AUGMENTED_LAGRANGE:
881 427214 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
882 427214 : pinfo->_normal(_component) * pinfo->_normal(_component);
883 :
884 0 : default:
885 0 : mooseError("Invalid contact formulation");
886 : }
887 :
888 205863 : case ContactModel::COULOMB:
889 205863 : switch (_formulation)
890 : {
891 43920 : case ContactFormulation::KINEMATIC:
892 : {
893 43920 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
894 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
895 : {
896 : RealVectorValue jac_vec;
897 130536 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
898 : {
899 87024 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
900 87024 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
901 87024 : _var_objects[i]->scalingFactor();
902 : }
903 43512 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
904 43512 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
905 43512 : pinfo->_normal(_component) * pinfo->_normal(_component);
906 : }
907 : else
908 : {
909 : const Real curr_jac =
910 408 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
911 408 : _connected_dof_indices[_j]) /
912 408 : _var.scalingFactor();
913 408 : return (-curr_jac + _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
914 : }
915 : }
916 :
917 48793 : case ContactFormulation::PENALTY:
918 : {
919 48793 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
920 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
921 10174 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
922 10174 : pinfo->_normal(_component) * pinfo->_normal(_component);
923 : else
924 38619 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp];
925 : }
926 72714 : case ContactFormulation::AUGMENTED_LAGRANGE:
927 : {
928 72714 : Real normal_comp = _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
929 72714 : pinfo->_normal(_component) * pinfo->_normal(_component);
930 :
931 : Real tang_comp = 0.0;
932 72714 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
933 25662 : tang_comp = _phi_secondary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
934 25662 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
935 72714 : return normal_comp + tang_comp;
936 : }
937 :
938 : case ContactFormulation::TANGENTIAL_PENALTY:
939 : {
940 : RealVectorValue jac_vec;
941 121308 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
942 : {
943 80872 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
944 80872 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
945 80872 : _var_objects[i]->scalingFactor();
946 : }
947 40436 : Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
948 40436 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
949 40436 : pinfo->_normal(_component) * pinfo->_normal(_component);
950 :
951 : Real tang_comp = 0.0;
952 40436 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
953 408 : tang_comp = _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
954 408 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
955 :
956 40436 : return normal_comp + tang_comp;
957 : }
958 :
959 0 : default:
960 0 : mooseError("Invalid contact formulation");
961 : }
962 :
963 171254 : case ContactModel::GLUED:
964 171254 : switch (_formulation)
965 : {
966 104281 : case ContactFormulation::KINEMATIC:
967 : {
968 104281 : const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
969 104281 : _connected_dof_indices[_j]) /
970 104281 : _var.scalingFactor();
971 104281 : return (-curr_jac + _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
972 : }
973 :
974 66973 : case ContactFormulation::PENALTY:
975 : case ContactFormulation::AUGMENTED_LAGRANGE:
976 66973 : 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 1238432 : case Moose::SecondaryPrimary:
986 1238432 : switch (_model)
987 : {
988 946188 : case ContactModel::FRICTIONLESS:
989 946188 : switch (_formulation)
990 : {
991 635616 : case ContactFormulation::KINEMATIC:
992 : {
993 635616 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
994 :
995 : RealVectorValue jac_vec;
996 2461632 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
997 : {
998 1826016 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
999 3652032 : jac_vec(i) = (*_jacobian)(dof_number,
1000 1826016 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1001 1826016 : _var_objects[i]->scalingFactor();
1002 : }
1003 635616 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1004 635616 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1005 635616 : pinfo->_normal(_component) * pinfo->_normal(_component);
1006 : }
1007 :
1008 310572 : case ContactFormulation::PENALTY:
1009 : case ContactFormulation::AUGMENTED_LAGRANGE:
1010 310572 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1011 310572 : pinfo->_normal(_component) * pinfo->_normal(_component);
1012 :
1013 0 : default:
1014 0 : mooseError("Invalid contact formulation");
1015 : }
1016 :
1017 164940 : case ContactModel::COULOMB:
1018 164940 : switch (_formulation)
1019 : {
1020 33528 : case ContactFormulation::KINEMATIC:
1021 : {
1022 33528 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1023 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1024 : {
1025 33240 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1026 :
1027 : RealVectorValue jac_vec;
1028 99720 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1029 : {
1030 66480 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1031 66480 : jac_vec(i) =
1032 66480 : (*_jacobian)(dof_number,
1033 66480 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1034 66480 : _var_objects[i]->scalingFactor();
1035 : }
1036 33240 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1037 33240 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1038 33240 : pinfo->_normal(_component) * pinfo->_normal(_component);
1039 : }
1040 : else
1041 : {
1042 288 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1043 : const Real curr_jac =
1044 288 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1045 288 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1046 288 : _var.scalingFactor();
1047 288 : return (-curr_jac - _phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
1048 : }
1049 : }
1050 :
1051 39604 : case ContactFormulation::PENALTY:
1052 : {
1053 39604 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1054 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1055 7620 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1056 7620 : pinfo->_normal(_component) * pinfo->_normal(_component);
1057 : else
1058 31984 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp];
1059 : }
1060 60040 : case ContactFormulation::AUGMENTED_LAGRANGE:
1061 : {
1062 60040 : Real normal_comp = -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1063 60040 : pinfo->_normal(_component) * pinfo->_normal(_component);
1064 :
1065 : Real tang_comp = 0.0;
1066 60040 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1067 22460 : tang_comp = -_phi_primary[_j][_qp] * penalty_slip * _test_secondary[_i][_qp] *
1068 22460 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1069 60040 : return normal_comp + tang_comp;
1070 : }
1071 :
1072 31768 : case ContactFormulation::TANGENTIAL_PENALTY:
1073 : {
1074 31768 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1075 :
1076 : RealVectorValue jac_vec;
1077 95304 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1078 : {
1079 63536 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1080 127072 : jac_vec(i) = (*_jacobian)(dof_number,
1081 63536 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1082 63536 : _var_objects[i]->scalingFactor();
1083 : }
1084 31768 : Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1085 31768 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1086 31768 : pinfo->_normal(_component) * pinfo->_normal(_component);
1087 :
1088 : Real tang_comp = 0.0;
1089 31768 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1090 304 : tang_comp = -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1091 304 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1092 :
1093 31768 : return normal_comp + tang_comp;
1094 : }
1095 :
1096 0 : default:
1097 0 : mooseError("Invalid contact formulation");
1098 : }
1099 127304 : case ContactModel::GLUED:
1100 127304 : switch (_formulation)
1101 : {
1102 73764 : case ContactFormulation::KINEMATIC:
1103 : {
1104 73764 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1105 : const Real curr_jac =
1106 73764 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1107 73764 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1108 73764 : _var.scalingFactor();
1109 73764 : return (-curr_jac - _phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]);
1110 : }
1111 :
1112 53540 : case ContactFormulation::PENALTY:
1113 : case ContactFormulation::AUGMENTED_LAGRANGE:
1114 53540 : 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 27758412 : case Moose::PrimarySecondary:
1125 27758412 : switch (_model)
1126 : {
1127 23553112 : case ContactModel::FRICTIONLESS:
1128 23553112 : switch (_formulation)
1129 : {
1130 : case ContactFormulation::KINEMATIC:
1131 : {
1132 : RealVectorValue jac_vec;
1133 71905536 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1134 : {
1135 53795948 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1136 53795948 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1137 53795948 : _var_objects[i]->scalingFactor();
1138 : }
1139 18109588 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1140 18109588 : _test_primary[_i][_qp];
1141 : }
1142 :
1143 5443524 : case ContactFormulation::PENALTY:
1144 : case ContactFormulation::AUGMENTED_LAGRANGE:
1145 5443524 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1146 5443524 : pinfo->_normal(_component) * pinfo->_normal(_component);
1147 :
1148 0 : default:
1149 0 : mooseError("Invalid contact formulation");
1150 : }
1151 2275740 : case ContactModel::COULOMB:
1152 2275740 : switch (_formulation)
1153 : {
1154 175680 : case ContactFormulation::KINEMATIC:
1155 : {
1156 175680 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1157 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1158 : {
1159 : RealVectorValue jac_vec;
1160 522144 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1161 : {
1162 348096 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1163 348096 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1164 348096 : _var_objects[i]->scalingFactor();
1165 : }
1166 174048 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1167 174048 : _test_primary[_i][_qp];
1168 : }
1169 : else
1170 : {
1171 : const Real secondary_jac =
1172 1632 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1173 1632 : _connected_dof_indices[_j]) /
1174 1632 : _var.scalingFactor();
1175 1632 : return secondary_jac * _test_primary[_i][_qp];
1176 : }
1177 : }
1178 :
1179 739972 : case ContactFormulation::PENALTY:
1180 : {
1181 739972 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1182 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1183 189336 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1184 189336 : pinfo->_normal(_component) * pinfo->_normal(_component);
1185 : else
1186 550636 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp];
1187 : }
1188 1198344 : case ContactFormulation::AUGMENTED_LAGRANGE:
1189 : {
1190 1198344 : Real normal_comp = -_phi_secondary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1191 1198344 : pinfo->_normal(_component) * pinfo->_normal(_component);
1192 :
1193 : Real tang_comp = 0.0;
1194 1198344 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1195 357984 : tang_comp = -_phi_secondary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1196 357984 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1197 1198344 : return normal_comp + tang_comp;
1198 : }
1199 :
1200 : case ContactFormulation::TANGENTIAL_PENALTY:
1201 : {
1202 : RealVectorValue jac_vec;
1203 485232 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1204 : {
1205 323488 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1206 323488 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1207 323488 : _var_objects[i]->scalingFactor();
1208 : }
1209 : Real normal_comp =
1210 161744 : pinfo->_normal(_component) * (pinfo->_normal * jac_vec) * _test_primary[_i][_qp];
1211 :
1212 : Real tang_comp = 0.0;
1213 161744 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1214 1632 : tang_comp = -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1215 1632 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1216 :
1217 161744 : return normal_comp + tang_comp;
1218 : }
1219 :
1220 0 : default:
1221 0 : mooseError("Invalid contact formulation");
1222 : }
1223 :
1224 1929560 : case ContactModel::GLUED:
1225 1929560 : switch (_formulation)
1226 : {
1227 1074844 : case ContactFormulation::KINEMATIC:
1228 : {
1229 : const Real secondary_jac =
1230 1074844 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1231 1074844 : _connected_dof_indices[_j]) /
1232 1074844 : _var.scalingFactor();
1233 1074844 : return secondary_jac * _test_primary[_i][_qp];
1234 : }
1235 :
1236 854716 : case ContactFormulation::PENALTY:
1237 : case ContactFormulation::AUGMENTED_LAGRANGE:
1238 854716 : 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 19538972 : case Moose::PrimaryPrimary:
1249 19538972 : switch (_model)
1250 : {
1251 16240780 : case ContactModel::FRICTIONLESS:
1252 16240780 : switch (_formulation)
1253 : {
1254 : case ContactFormulation::KINEMATIC:
1255 : return 0.0;
1256 :
1257 4034800 : case ContactFormulation::PENALTY:
1258 : case ContactFormulation::AUGMENTED_LAGRANGE:
1259 4034800 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1260 4034800 : pinfo->_normal(_component) * pinfo->_normal(_component);
1261 :
1262 0 : default:
1263 0 : mooseError("Invalid contact formulation");
1264 : }
1265 :
1266 3298192 : case ContactModel::COULOMB:
1267 : case ContactModel::GLUED:
1268 3298192 : switch (_formulation)
1269 : {
1270 : case ContactFormulation::KINEMATIC:
1271 : return 0.0;
1272 :
1273 1281696 : case ContactFormulation::PENALTY:
1274 : {
1275 1281696 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1276 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1277 141200 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1278 141200 : pinfo->_normal(_component) * pinfo->_normal(_component);
1279 : else
1280 1140496 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp];
1281 : }
1282 :
1283 127072 : case ContactFormulation::TANGENTIAL_PENALTY:
1284 : {
1285 : Real tang_comp = 0.0;
1286 127072 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1287 1216 : tang_comp = _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1288 1216 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1289 : return tang_comp; // normal component is zero
1290 : }
1291 :
1292 979040 : case ContactFormulation::AUGMENTED_LAGRANGE:
1293 : {
1294 979040 : Real normal_comp = _phi_primary[_j][_qp] * penalty * _test_primary[_i][_qp] *
1295 979040 : pinfo->_normal(_component) * pinfo->_normal(_component);
1296 :
1297 : Real tang_comp = 0.0;
1298 979040 : if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1299 313840 : tang_comp = _phi_primary[_j][_qp] * penalty_slip * _test_primary[_i][_qp] *
1300 313840 : (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1301 979040 : 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 7939266 : MechanicalContactConstraint::computeQpOffDiagJacobian(Moose::ConstraintJacobianType type,
1318 : unsigned int jvar)
1319 : {
1320 7939266 : PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
1321 :
1322 7939266 : const Real penalty = getPenalty(*_current_node);
1323 7939266 : const Real penalty_slip = getTangentialPenalty(*_current_node);
1324 :
1325 : unsigned int coupled_component;
1326 : Real normal_component_in_coupled_var_dir = 1.0;
1327 7939266 : if (getCoupledVarComponent(jvar, coupled_component))
1328 7939266 : normal_component_in_coupled_var_dir = pinfo->_normal(coupled_component);
1329 :
1330 7939266 : switch (type)
1331 : {
1332 0 : default:
1333 0 : mooseError("Unhandled ConstraintJacobianType");
1334 :
1335 399354 : case Moose::SecondarySecondary:
1336 399354 : switch (_model)
1337 : {
1338 317532 : case ContactModel::FRICTIONLESS:
1339 317532 : switch (_formulation)
1340 : {
1341 : case ContactFormulation::KINEMATIC:
1342 : {
1343 : RealVectorValue jac_vec;
1344 661056 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1345 : {
1346 491028 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1347 491028 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1348 491028 : _var_objects[i]->scalingFactor();
1349 : }
1350 170028 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1351 170028 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1352 170028 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1353 : }
1354 :
1355 147504 : case ContactFormulation::PENALTY:
1356 : case ContactFormulation::AUGMENTED_LAGRANGE:
1357 147504 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1358 147504 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1359 :
1360 0 : default:
1361 0 : mooseError("Invalid contact formulation");
1362 : }
1363 :
1364 81822 : case ContactModel::COULOMB:
1365 : {
1366 81822 : if ((_formulation == ContactFormulation::KINEMATIC ||
1367 75520 : _formulation == ContactFormulation::TANGENTIAL_PENALTY) &&
1368 75520 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1369 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1370 : {
1371 : RealVectorValue jac_vec;
1372 224832 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1373 : {
1374 149888 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1375 149888 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1376 149888 : _var_objects[i]->scalingFactor();
1377 : }
1378 :
1379 74944 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1380 74944 : (_phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1381 74944 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1382 : }
1383 6878 : else if ((_formulation == ContactFormulation::PENALTY) &&
1384 6302 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1385 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1386 552 : return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1387 552 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1388 6326 : 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 6326 : const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1402 6326 : _connected_dof_indices[_j]);
1403 6326 : 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 295760 : case Moose::SecondaryPrimary:
1418 295760 : switch (_model)
1419 : {
1420 232432 : case ContactModel::FRICTIONLESS:
1421 232432 : switch (_formulation)
1422 : {
1423 128368 : case ContactFormulation::KINEMATIC:
1424 : {
1425 128368 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1426 :
1427 : RealVectorValue jac_vec;
1428 498096 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1429 : {
1430 369728 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1431 739456 : jac_vec(i) = (*_jacobian)(dof_number,
1432 369728 : curr_primary_node->dof_number(0, _vars[_component], 0)) /
1433 369728 : _var_objects[i]->scalingFactor();
1434 : }
1435 128368 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1436 128368 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1437 128368 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1438 : }
1439 :
1440 104064 : case ContactFormulation::PENALTY:
1441 : case ContactFormulation::AUGMENTED_LAGRANGE:
1442 104064 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1443 104064 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1444 :
1445 0 : default:
1446 0 : mooseError("Invalid contact formulation");
1447 : }
1448 :
1449 63328 : case ContactModel::COULOMB:
1450 63328 : if ((_formulation == ContactFormulation::KINEMATIC ||
1451 59200 : _formulation == ContactFormulation::TANGENTIAL_PENALTY) &&
1452 59200 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1453 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1454 : {
1455 58768 : const Node * curr_primary_node = _current_primary->node_ptr(_j);
1456 :
1457 : RealVectorValue jac_vec;
1458 176304 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1459 : {
1460 117536 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1461 117536 : jac_vec(i) =
1462 117536 : (*_jacobian)(dof_number, curr_primary_node->dof_number(0, _vars[_component], 0)) /
1463 117536 : _var_objects[i]->scalingFactor();
1464 : }
1465 :
1466 58768 : return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1467 58768 : (_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp]) *
1468 58768 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1469 : }
1470 4560 : else if ((_formulation == ContactFormulation::PENALTY) &&
1471 4128 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1472 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1473 368 : return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp] *
1474 368 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1475 4192 : 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 4172920 : case Moose::PrimarySecondary:
1497 4172920 : switch (_model)
1498 : {
1499 3845632 : case ContactModel::FRICTIONLESS:
1500 3845632 : switch (_formulation)
1501 : {
1502 : case ContactFormulation::KINEMATIC:
1503 : {
1504 : RealVectorValue jac_vec;
1505 11808672 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1506 : {
1507 8829536 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1508 8829536 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1509 8829536 : _var_objects[i]->scalingFactor();
1510 : }
1511 2979136 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1512 2979136 : _test_primary[_i][_qp];
1513 : }
1514 :
1515 866496 : case ContactFormulation::PENALTY:
1516 : case ContactFormulation::AUGMENTED_LAGRANGE:
1517 866496 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1518 866496 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1519 :
1520 0 : default:
1521 0 : mooseError("Invalid contact formulation");
1522 : }
1523 327288 : case ContactModel::COULOMB:
1524 327288 : switch (_formulation)
1525 : {
1526 157664 : case ContactFormulation::KINEMATIC:
1527 : {
1528 157664 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1529 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1530 : {
1531 : RealVectorValue jac_vec;
1532 469824 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1533 : {
1534 313216 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1535 313216 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1536 313216 : _var_objects[i]->scalingFactor();
1537 : }
1538 156608 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1539 156608 : _test_primary[_i][_qp];
1540 : }
1541 : else
1542 : {
1543 : const Real secondary_jac =
1544 1056 : (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1545 1056 : _connected_dof_indices[_j]) /
1546 1056 : _var.scalingFactor();
1547 1056 : return secondary_jac * _test_primary[_i][_qp];
1548 : }
1549 : }
1550 :
1551 25208 : case ContactFormulation::PENALTY:
1552 : {
1553 25208 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1554 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1555 2208 : return -_test_primary[_i][_qp] * penalty * _phi_secondary[_j][_qp] *
1556 2208 : 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 144416 : case ContactFormulation::TANGENTIAL_PENALTY:
1572 : {
1573 144416 : if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1574 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1575 : {
1576 : RealVectorValue jac_vec;
1577 429504 : for (unsigned int i = 0; i < _mesh_dimension; ++i)
1578 : {
1579 286336 : dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1580 286336 : jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]) /
1581 286336 : _var_objects[i]->scalingFactor();
1582 : }
1583 143168 : return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1584 143168 : _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 3071232 : case Moose::PrimaryPrimary:
1619 3071232 : switch (_model)
1620 : {
1621 2817920 : case ContactModel::FRICTIONLESS:
1622 2817920 : switch (_formulation)
1623 : {
1624 : case ContactFormulation::KINEMATIC:
1625 : return 0.0;
1626 :
1627 582144 : case ContactFormulation::PENALTY:
1628 : case ContactFormulation::AUGMENTED_LAGRANGE:
1629 582144 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1630 582144 : pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1631 :
1632 0 : default:
1633 0 : mooseError("Invalid contact formulation");
1634 : }
1635 :
1636 253312 : case ContactModel::COULOMB:
1637 : {
1638 253312 : 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 253312 : else if (_formulation == ContactFormulation::PENALTY &&
1650 16512 : (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1651 : pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1652 1472 : return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp] *
1653 1472 : 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 1650512 : MechanicalContactConstraint::gapOffset(const Node & node)
1688 : {
1689 : Real val = 0;
1690 :
1691 1650512 : if (_has_secondary_gap_offset)
1692 53298 : val += _secondary_gap_offset_var->getNodalValue(node);
1693 :
1694 1650512 : if (_has_mapped_primary_gap_offset)
1695 53298 : val += _mapped_primary_gap_offset_var->getNodalValue(node);
1696 :
1697 1650512 : return val;
1698 : }
1699 :
1700 : Real
1701 26650188 : MechanicalContactConstraint::nodalArea(const Node & node)
1702 : {
1703 26650188 : dof_id_type dof = node.dof_number(_aux_system.number(), _nodal_area_var->number(), 0);
1704 :
1705 26650188 : Real area = (*_aux_solution)(dof);
1706 26650188 : 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 26650188 : return area;
1715 : }
1716 :
1717 : Real
1718 59565234 : MechanicalContactConstraint::getPenalty(const Node & node)
1719 : {
1720 59565234 : Real penalty = _penalty;
1721 59565234 : if (_normalize_penalty)
1722 13325483 : penalty *= nodalArea(node);
1723 59565234 : return penalty * MathUtils::pow(_penalty_multiplier, _lagrangian_iteration_number);
1724 : }
1725 :
1726 : Real
1727 58794640 : MechanicalContactConstraint::getTangentialPenalty(const Node & node)
1728 : {
1729 58794640 : Real penalty = _penalty_tangential;
1730 58794640 : if (_normalize_penalty)
1731 13272426 : penalty *= nodalArea(node);
1732 :
1733 58794640 : return penalty * MathUtils::pow(_penalty_multiplier, _lagrangian_iteration_number);
1734 : }
1735 :
1736 : void
1737 126914 : MechanicalContactConstraint::computeJacobian()
1738 : {
1739 126914 : getConnectedDofIndices(_var.number());
1740 :
1741 126914 : prepareMatrixTagNeighbor(
1742 126914 : _assembly, _primary_var.number(), _var.number(), Moose::NeighborNeighbor, _Knn);
1743 :
1744 126914 : _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
1745 :
1746 253828 : for (_i = 0; _i < _test_secondary.size(); _i++)
1747 : // Loop over the connected dof indices so we can get all the jacobian contributions
1748 1862607 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1749 1735693 : _Kee(_i, _j) += computeQpJacobian(Moose::SecondarySecondary);
1750 :
1751 126914 : if (_primary_secondary_jacobian)
1752 : {
1753 124793 : prepareMatrixTagNeighbor(_assembly, _var.number(), _var.number(), Moose::ElementNeighbor, _Ken);
1754 124793 : if (_Ken.m() && _Ken.n())
1755 : {
1756 249586 : for (_i = 0; _i < _test_secondary.size(); _i++)
1757 1363225 : for (_j = 0; _j < _phi_primary.size(); _j++)
1758 1238432 : _Ken(_i, _j) += computeQpJacobian(Moose::SecondaryPrimary);
1759 124793 : accumulateTaggedLocalMatrix(
1760 124793 : _assembly, _var.number(), _var.number(), Moose::ElementNeighbor, _Ken);
1761 : }
1762 :
1763 124793 : _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
1764 1363225 : for (_i = 0; _i < _test_primary.size(); _i++)
1765 : // Loop over the connected dof indices so we can get all the jacobian contributions
1766 28996844 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1767 27758412 : _Kne(_i, _j) += computeQpJacobian(Moose::PrimarySecondary);
1768 : }
1769 :
1770 126914 : if (_Knn.m() && _Knn.n())
1771 : {
1772 1382314 : for (_i = 0; _i < _test_primary.size(); _i++)
1773 20794372 : for (_j = 0; _j < _phi_primary.size(); _j++)
1774 19538972 : _Knn(_i, _j) += computeQpJacobian(Moose::PrimaryPrimary);
1775 126914 : accumulateTaggedLocalMatrix(
1776 126914 : _assembly, _primary_var.number(), _var.number(), Moose::NeighborNeighbor, _Knn);
1777 : }
1778 126914 : }
1779 :
1780 : void
1781 46110 : MechanicalContactConstraint::computeOffDiagJacobian(const unsigned int jvar_num)
1782 : {
1783 46110 : getConnectedDofIndices(jvar_num);
1784 :
1785 46110 : _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
1786 :
1787 46110 : prepareMatrixTagNeighbor(
1788 46110 : _assembly, _primary_var.number(), jvar_num, Moose::NeighborNeighbor, _Knn);
1789 :
1790 92220 : for (_i = 0; _i < _test_secondary.size(); _i++)
1791 : // Loop over the connected dof indices so we can get all the jacobian contributions
1792 445464 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1793 399354 : _Kee(_i, _j) += computeQpOffDiagJacobian(Moose::SecondarySecondary, jvar_num);
1794 :
1795 46110 : if (_primary_secondary_jacobian)
1796 : {
1797 46110 : prepareMatrixTagNeighbor(_assembly, _var.number(), jvar_num, Moose::ElementNeighbor, _Ken);
1798 92220 : for (_i = 0; _i < _test_secondary.size(); _i++)
1799 341870 : for (_j = 0; _j < _phi_primary.size(); _j++)
1800 295760 : _Ken(_i, _j) += computeQpOffDiagJacobian(Moose::SecondaryPrimary, jvar_num);
1801 46110 : accumulateTaggedLocalMatrix(_assembly, _var.number(), jvar_num, Moose::ElementNeighbor, _Ken);
1802 :
1803 46110 : _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
1804 46110 : if (_Kne.m() && _Kne.n())
1805 341870 : for (_i = 0; _i < _test_primary.size(); _i++)
1806 : // Loop over the connected dof indices so we can get all the jacobian contributions
1807 4468680 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1808 4172920 : _Kne(_i, _j) += computeQpOffDiagJacobian(Moose::PrimarySecondary, jvar_num);
1809 : }
1810 :
1811 341870 : for (_i = 0; _i < _test_primary.size(); _i++)
1812 3366992 : for (_j = 0; _j < _phi_primary.size(); _j++)
1813 3071232 : _Knn(_i, _j) += computeQpOffDiagJacobian(Moose::PrimaryPrimary, jvar_num);
1814 46110 : accumulateTaggedLocalMatrix(
1815 46110 : _assembly, _primary_var.number(), jvar_num, Moose::NeighborNeighbor, _Knn);
1816 46110 : }
1817 :
1818 : void
1819 173024 : MechanicalContactConstraint::getConnectedDofIndices(unsigned int var_num)
1820 : {
1821 : unsigned int component;
1822 173024 : if (getCoupledVarComponent(var_num, component) || _non_displacement_vars_jacobian)
1823 : {
1824 173024 : if (_primary_secondary_jacobian && _connected_secondary_nodes_jacobian)
1825 170903 : NodeFaceConstraint::getConnectedDofIndices(var_num);
1826 : else
1827 : {
1828 2121 : _connected_dof_indices.clear();
1829 2121 : MooseVariableFEBase & var = _sys.getVariable(0, var_num);
1830 2121 : _connected_dof_indices.push_back(var.nodalDofIndex());
1831 : }
1832 : }
1833 :
1834 173024 : _phi_secondary.resize(_connected_dof_indices.size());
1835 :
1836 173024 : dof_id_type current_node_var_dof_index = _sys.getVariable(0, var_num).nodalDofIndex();
1837 173024 : _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 2308071 : for (unsigned int j = 0; j < _connected_dof_indices.size(); j++)
1843 : {
1844 2135047 : _phi_secondary[j].resize(1);
1845 :
1846 2135047 : if (_connected_dof_indices[j] == current_node_var_dof_index)
1847 173024 : _phi_secondary[j][_qp] = 1.0;
1848 : else
1849 1962023 : _phi_secondary[j][_qp] = 0.0;
1850 : }
1851 173024 : }
1852 :
1853 : bool
1854 8112290 : MechanicalContactConstraint::getCoupledVarComponent(unsigned int var_num, unsigned int & component)
1855 : {
1856 8112290 : component = std::numeric_limits<unsigned int>::max();
1857 : bool coupled_var_is_disp_var = false;
1858 15343530 : for (const auto i : make_range(Moose::dim))
1859 : {
1860 15343530 : if (var_num == _vars[i])
1861 : {
1862 : coupled_var_is_disp_var = true;
1863 8112290 : component = i;
1864 8112290 : break;
1865 : }
1866 : }
1867 :
1868 8112290 : return coupled_var_is_disp_var;
1869 : }
1870 :
1871 : void
1872 49888 : MechanicalContactConstraint::residualEnd()
1873 : {
1874 49888 : if (_component == 0 && (_print_contact_nodes || _contact_linesearch))
1875 : {
1876 66 : _communicator.set_union(_current_contact_state);
1877 66 : 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 66 : if (_contact_linesearch)
1887 66 : _contact_linesearch->insertSet(_current_contact_state);
1888 : _old_contact_state.swap(_current_contact_state);
1889 : _current_contact_state.clear();
1890 : }
1891 49888 : }
|
