www.mooseframework.org
MechanicalContactConstraint.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 // MOOSE includes
10 #include "FEProblem.h"
11 #include "DisplacedProblem.h"
12 #include "AuxiliarySystem.h"
13 #include "PenetrationLocator.h"
14 #include "NearestNodeLocator.h"
15 
16 #include "SystemBase.h"
17 #include "Assembly.h"
18 #include "MooseMesh.h"
20 #include "Executioner.h"
21 #include "AddVariableAction.h"
22 
23 #include "libmesh/string_to_enum.h"
24 #include "libmesh/sparse_matrix.h"
25 
26 template <>
27 InputParameters
29 {
30  MooseEnum orders(AddVariableAction::getNonlinearVariableOrders());
31 
32  InputParameters params = validParams<NodeFaceConstraint>();
33  params.addRequiredParam<BoundaryName>("boundary", "The master boundary");
34  params.addRequiredParam<BoundaryName>("slave", "The slave boundary");
35  params.addRequiredParam<unsigned int>("component",
36  "An integer corresponding to the direction "
37  "the variable this kernel acts in. (0 for x, "
38  "1 for y, 2 for z)");
39 
40  params.addCoupledVar("disp_x", "The x displacement");
41  params.addCoupledVar("disp_y", "The y displacement");
42  params.addCoupledVar("disp_z", "The z displacement");
43 
44  params.addCoupledVar(
45  "displacements",
46  "The displacements appropriate for the simulation geometry and coordinate system");
47 
48  params.addRequiredCoupledVar("nodal_area", "The nodal area");
49  params.addParam<std::string>("model", "frictionless", "The contact model to use");
50 
51  params.set<bool>("use_displaced_mesh") = true;
52  params.addParam<Real>(
53  "penalty",
54  1e8,
55  "The penalty to apply. This can vary depending on the stiffness of your materials");
56  params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
57  params.addParam<Real>("tangential_tolerance",
58  "Tangential distance to extend edges of contact surfaces");
59  params.addParam<Real>(
60  "capture_tolerance", 0, "Normal distance from surface within which nodes are captured");
61  params.addParam<Real>(
62  "normal_smoothing_distance",
63  "Distance from edge in parametric coordinates over which to smooth contact normal");
64  params.addParam<std::string>("normal_smoothing_method",
65  "Method to use to smooth normals (edge_based|nodal_normal_based)");
66  params.addParam<MooseEnum>("order", orders, "The finite element order");
67 
68  params.addParam<Real>("tension_release",
69  0.0,
70  "Tension release threshold. A node in contact "
71  "will not be released if its tensile load is below "
72  "this value. No tension release if negative.");
73 
74  params.addParam<std::string>("formulation", "default", "The contact formulation");
75  params.addParam<bool>(
76  "normalize_penalty",
77  false,
78  "Whether to normalize the penalty parameter with the nodal area for penalty contact.");
79  params.addParam<bool>("master_slave_jacobian",
80  true,
81  "Whether to include jacobian entries coupling master and slave nodes.");
82  params.addParam<bool>(
83  "connected_slave_nodes_jacobian",
84  true,
85  "Whether to include jacobian entries coupling nodes connected to slave nodes.");
86  params.addParam<bool>("non_displacement_variables_jacobian",
87  true,
88  "Whether to include jacobian "
89  "entries coupling with "
90  "variables that are not "
91  "displacement variables.");
92  params.addParam<unsigned int>("stick_lock_iterations",
93  std::numeric_limits<unsigned int>::max(),
94  "Number of times permitted to switch between sticking and slipping "
95  "in a solution before locking node in a sticked state.");
96  params.addParam<Real>("stick_unlock_factor",
97  1.5,
98  "Factor by which frictional capacity must be "
99  "exceeded to permit stick-locked node to slip "
100  "again.");
101  params.addParam<Real>("al_penetration_tolerance",
102  "The tolerance of the penetration for augmented Lagrangian method.");
103  params.addParam<Real>("al_incremental_slip_tolerance",
104  "The tolerance of the incremental slip for augmented Lagrangian method.");
105 
106  params.addParam<Real>("al_frictional_force_tolerance",
107  "The tolerance of the frictional force for augmented Lagrangian method.");
108  return params;
109 }
110 
112  : NodeFaceConstraint(parameters),
113  _displaced_problem(parameters.get<FEProblemBase *>("_fe_problem_base")->getDisplacedProblem()),
114  _component(getParam<unsigned int>("component")),
115  _model(ContactMaster::contactModel(getParam<std::string>("model"))),
116  _formulation(ContactMaster::contactFormulation(getParam<std::string>("formulation"))),
117  _normalize_penalty(getParam<bool>("normalize_penalty")),
118  _penalty(getParam<Real>("penalty")),
119  _friction_coefficient(getParam<Real>("friction_coefficient")),
120  _tension_release(getParam<Real>("tension_release")),
121  _capture_tolerance(getParam<Real>("capture_tolerance")),
122  _stick_lock_iterations(getParam<unsigned int>("stick_lock_iterations")),
123  _stick_unlock_factor(getParam<Real>("stick_unlock_factor")),
124  _update_contact_set(true),
125  _residual_copy(_sys.residualGhosted()),
126  _mesh_dimension(_mesh.dimension()),
127  _vars(3, libMesh::invalid_uint),
128  _nodal_area_var(getVar("nodal_area", 0)),
129  _aux_system(_nodal_area_var->sys()),
130  _aux_solution(_aux_system.currentSolution()),
131  _master_slave_jacobian(getParam<bool>("master_slave_jacobian")),
132  _connected_slave_nodes_jacobian(getParam<bool>("connected_slave_nodes_jacobian")),
133  _non_displacement_vars_jacobian(getParam<bool>("non_displacement_variables_jacobian"))
134 {
135  _overwrite_slave_residual = false;
136 
137  if (isParamValid("displacements"))
138  {
139  // modern parameter scheme for displacements
140  for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
141  _vars[i] = coupled("displacements", i);
142  }
143  else
144  {
145  // Legacy parameter scheme for displacements
146  if (isParamValid("disp_x"))
147  _vars[0] = coupled("disp_x");
148  if (isParamValid("disp_y"))
149  _vars[1] = coupled("disp_y");
150  if (isParamValid("disp_z"))
151  _vars[2] = coupled("disp_z");
152 
153  mooseDeprecated("use the `displacements` parameter rather than the `disp_*` parameters (those "
154  "will go away with the deprecation of the Solid Mechanics module).");
155  }
156 
157  if (parameters.isParamValid("tangential_tolerance"))
158  _penetration_locator.setTangentialTolerance(getParam<Real>("tangential_tolerance"));
159 
160  if (parameters.isParamValid("normal_smoothing_distance"))
161  _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
162 
163  if (parameters.isParamValid("normal_smoothing_method"))
164  _penetration_locator.setNormalSmoothingMethod(
165  parameters.get<std::string>("normal_smoothing_method"));
166 
168  mooseError("The 'tangential_penalty' formulation can only be used with the 'coulomb' model");
169 
170  if (_model == CM_GLUED)
171  _penetration_locator.setUpdate(false);
172 
173  if (_friction_coefficient < 0)
174  mooseError("The friction coefficient must be nonnegative");
175 
176  // set _penalty_tangential to the value of _penalty for now
178 
180  {
181  if (_model == CM_GLUED)
182  mooseError("The Augmented Lagrangian contact formulation does not support GLUED case.");
183 
184  FEProblemBase * fe_problem = getParam<FEProblemBase *>("_fe_problem_base");
185  if (dynamic_cast<AugmentedLagrangianContactProblem *>(fe_problem) == NULL)
186  mooseError("The Augmented Lagrangian contact formulation must use "
187  "AugmentedLagrangianContactProblem.");
188 
189  if (!parameters.isParamValid("al_penetration_tolerance"))
190  mooseError("For Augmented Lagrangian contact, al_penetration_tolerance must be provided.");
191  else
192  _al_penetration_tolerance = parameters.get<Real>("al_penetration_tolerance");
193 
194  if (_model != CM_FRICTIONLESS)
195  {
196  if (!parameters.isParamValid("al_incremental_slip_tolerance") ||
197  !parameters.isParamValid("al_frictional_force_tolerance"))
198  {
199  mooseError("For the Augmented Lagrangian frictional contact formualton, "
200  "al_incremental_slip_tolerance and "
201  "al_frictional_force_tolerance must be provided.");
202  }
203  else
204  {
205  _al_incremental_slip_tolerance = parameters.get<Real>("al_incremental_slip_tolerance");
206  _al_frictional_force_tolerance = parameters.get<Real>("al_frictional_force_tolerance");
207  }
208  }
209  }
210 }
211 
212 void
214 {
215  if (_component == 0)
216  {
217  updateContactSet(true);
220 
221  _update_contact_set = false;
222  }
223 }
224 
225 void
227 {
228  if (_component == 0)
229  {
232  _update_contact_set = true;
233  }
234 }
235 
236 void
238 {
239  for (auto & pinfo_pair : _penetration_locator._penetration_info)
240  {
241  const dof_id_type slave_node_num = pinfo_pair.first;
242  PenetrationInfo * pinfo = pinfo_pair.second;
243 
244  if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
245  continue;
246 
247  const Real distance = pinfo->_normal * (pinfo->_closest_point - _mesh.nodeRef(slave_node_num));
248 
249  if (beginning_of_step && _model == CM_COULOMB)
250  {
251  pinfo->_lagrange_multiplier_slip.zero();
252  if (pinfo->isCaptured())
253  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
254  }
255 
256  if (pinfo->isCaptured())
257  {
258  if (_model == CM_FRICTIONLESS)
259  pinfo->_lagrange_multiplier -= getPenalty(*pinfo) * distance;
260 
261  if (_model == CM_COULOMB)
262  {
263  if (!beginning_of_step)
264  {
265  Real penalty = getPenalty(*pinfo);
266  RealVectorValue pen_force_normal =
267  penalty * (-distance) * pinfo->_normal + pinfo->_lagrange_multiplier * pinfo->_normal;
268 
269  // update normal lagrangian multiplier
270  pinfo->_lagrange_multiplier += penalty * (-distance);
271 
272  // Frictional capacity
273  const Real capacity(_friction_coefficient * (pen_force_normal * pinfo->_normal < 0
274  ? -pen_force_normal * pinfo->_normal
275  : 0));
276 
277  RealVectorValue tangential_inc_slip =
278  pinfo->_incremental_slip -
279  (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
280 
281  Real penalty_slip = getTangentialPenalty(*pinfo);
282 
283  RealVectorValue inc_pen_force_tangential =
284  pinfo->_lagrange_multiplier_slip + penalty_slip * tangential_inc_slip;
285 
286  RealVectorValue tau_old = pinfo->_contact_force_old -
287  pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old);
288 
289  RealVectorValue contact_force_tangential = inc_pen_force_tangential + tau_old;
290  const Real tan_mag(contact_force_tangential.norm());
291 
292  if (tan_mag > capacity * (_al_frictional_force_tolerance + 1.0))
293  {
294  pinfo->_lagrange_multiplier_slip =
295  -tau_old + capacity * contact_force_tangential / tan_mag;
296  if (MooseUtils::absoluteFuzzyEqual(capacity, 0.0))
297  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
298  else
299  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
300  }
301  else
302  {
303  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
304  pinfo->_lagrange_multiplier_slip += penalty_slip * tangential_inc_slip;
305  }
306  }
307  }
308  }
309  }
310 }
311 
312 bool
314 {
315  Real contactResidual = 0.0;
316  unsigned int converged = 0;
317 
318  for (auto & pinfo_pair : _penetration_locator._penetration_info)
319  {
320  const dof_id_type slave_node_num = pinfo_pair.first;
321  PenetrationInfo * pinfo = pinfo_pair.second;
322 
323  // Skip this pinfo if there are no DOFs on this node.
324  if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
325  continue;
326 
327  const Real distance = pinfo->_normal * (pinfo->_closest_point - _mesh.nodeRef(slave_node_num));
328 
329  if (pinfo->isCaptured())
330  {
331  if (contactResidual < std::abs(distance))
332  contactResidual = std::abs(distance);
333 
334  // penetration < tol
335  if (contactResidual > _al_penetration_tolerance)
336  {
337  converged = 1;
338  break;
339  }
340 
341  if (_model == CM_COULOMB)
342  {
343  RealVectorValue contact_force_normal((pinfo->_contact_force * pinfo->_normal) *
344  pinfo->_normal);
345  RealVectorValue contact_force_tangential(pinfo->_contact_force - contact_force_normal);
346 
347  RealVectorValue tangential_inc_slip =
348  pinfo->_incremental_slip - (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
349 
350  const Real tan_mag(contact_force_tangential.norm());
351  const Real tangential_inc_slip_mag = tangential_inc_slip.norm();
352 
353  RealVectorValue distance_vec =
354  (pinfo->_normal * (_mesh.nodeRef(slave_node_num) - pinfo->_closest_point)) *
355  pinfo->_normal;
356 
357  Real penalty = getPenalty(*pinfo);
358  RealVectorValue pen_force_normal =
359  penalty * distance_vec + pinfo->_lagrange_multiplier * pinfo->_normal;
360 
361  // Frictional capacity
362  Real capacity(_friction_coefficient * (pen_force_normal * pinfo->_normal < 0
363  ? -pen_force_normal * pinfo->_normal
364  : 0.0));
365 
366  // incremental slip <= tol for all pinfo_pair such that tan_mag < capacity
367  if (MooseUtils::absoluteFuzzyLessThan(tan_mag, capacity) &&
368  pinfo->_mech_status == PenetrationInfo::MS_STICKING)
369  {
370  if (MooseUtils::absoluteFuzzyGreaterThan(tangential_inc_slip_mag,
372  {
373  converged = 2;
374  break;
375  }
376  }
377 
378  // for all pinfo_pair, tag_mag should be less than (1 + tol) * (capacity + tol)
379  if (tan_mag >
381  {
382  converged = 3;
383  break;
384  }
385  }
386  }
387  }
388 
389  _communicator.max(converged);
390 
391  if (converged == 1)
392  _console
393  << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied \n";
394  else if (converged == 2)
395  _console
396  << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied \n";
397  else if (converged == 3)
398  _console << "The Augmented Lagrangian contact frictional force enforcement is NOT satisfied \n";
399  else
400  return true;
401 
402  return false;
403 }
404 
405 void
407 {
408  for (auto & pinfo_pair : _penetration_locator._penetration_info)
409  {
410  const dof_id_type slave_node_num = pinfo_pair.first;
411  PenetrationInfo * pinfo = pinfo_pair.second;
412 
413  // Skip this pinfo if there are no DOFs on this node.
414  if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
415  continue;
416 
417  if (beginning_of_step)
418  {
419  if (_app.getExecutioner()->lastSolveConverged())
420  {
421  pinfo->_contact_force_old = pinfo->_contact_force;
422  pinfo->_accumulated_slip_old = pinfo->_accumulated_slip;
423  pinfo->_frictional_energy_old = pinfo->_frictional_energy;
424  pinfo->_mech_status_old = pinfo->_mech_status;
425  }
426  else if (pinfo->_mech_status_old == PenetrationInfo::MS_NO_CONTACT &&
427  pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
428  {
429  // The penetration info object could be based on a bad state so delete it
430  delete pinfo_pair.second;
431  pinfo_pair.second = NULL;
432  continue;
433  }
434 
435  pinfo->_locked_this_step = 0;
436  pinfo->_stick_locked_this_step = 0;
437  pinfo->_starting_elem = pinfo->_elem;
438  pinfo->_starting_side_num = pinfo->_side_num;
439  pinfo->_starting_closest_point_ref = pinfo->_closest_point_ref;
440  }
441  pinfo->_incremental_slip_prev_iter = pinfo->_incremental_slip;
442 
443  const Real contact_pressure = -(pinfo->_normal * pinfo->_contact_force) / nodalArea(*pinfo);
444  const Real distance = pinfo->_normal * (pinfo->_closest_point - _mesh.nodeRef(slave_node_num));
445 
446  // Capture
447  if (!pinfo->isCaptured() &&
448  MooseUtils::absoluteFuzzyGreaterEqual(distance, 0.0, _capture_tolerance))
449  {
450  pinfo->capture();
451 
452  // Increment the lock count every time the node comes back into contact from not being in
453  // contact.
455  ++pinfo->_locked_this_step;
456  }
457  // Release
458  else if (_model != CM_GLUED && pinfo->isCaptured() && _tension_release >= 0.0 &&
459  -contact_pressure >= _tension_release && pinfo->_locked_this_step < 2)
460  {
461  pinfo->release();
462  pinfo->_contact_force.zero();
463  }
464  }
465 }
466 
467 bool
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 && pinfo->isCaptured())
478  {
479  in_contact = true;
480 
481  // This computes the contact force once per constraint, rather than once per quad point
482  // and for both master and slave cases.
483  if (_component == 0)
484  computeContactForce(pinfo);
485  }
486  }
487 
488  return in_contact;
489 }
490 
491 void
493 {
494  const Node * node = pinfo->_node;
495 
496  // Build up residual vector
497  RealVectorValue res_vec;
498  for (unsigned int i = 0; i < _mesh_dimension; ++i)
499  {
500  dof_id_type dof_number = node->dof_number(0, _vars[i], 0);
501  res_vec(i) = _residual_copy(dof_number);
502  }
503 
504  RealVectorValue distance_vec(_mesh.nodeRef(node->id()) - pinfo->_closest_point);
505  const Real penalty = getPenalty(*pinfo);
506  const Real penalty_slip = getTangentialPenalty(*pinfo);
507 
508  RealVectorValue pen_force(penalty * distance_vec);
509 
510  switch (_model)
511  {
512  case CM_FRICTIONLESS:
513  switch (_formulation)
514  {
515  case CF_KINEMATIC:
516  pinfo->_contact_force = -pinfo->_normal * (pinfo->_normal * res_vec);
517  break;
518 
519  case CF_PENALTY:
520  pinfo->_contact_force = pinfo->_normal * (pinfo->_normal * pen_force);
521  break;
522 
524  pinfo->_contact_force =
525  (pinfo->_normal *
526  (pinfo->_normal * (pen_force + pinfo->_lagrange_multiplier * pinfo->_normal)));
527  break;
528 
529  default:
530  mooseError("Invalid contact formulation");
531  break;
532  }
533  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
534  break;
535  case CM_COULOMB:
536  switch (_formulation)
537  {
538  case CF_KINEMATIC:
539  {
540  // Frictional capacity
541  const Real capacity(_friction_coefficient *
542  (res_vec * pinfo->_normal > 0 ? res_vec * pinfo->_normal : 0));
543 
544  // Normal and tangential components of predictor force
545  pinfo->_contact_force = -res_vec;
546  RealVectorValue contact_force_normal((pinfo->_contact_force * pinfo->_normal) *
547  pinfo->_normal);
548  RealVectorValue contact_force_tangential(pinfo->_contact_force - contact_force_normal);
549 
550  RealVectorValue tangential_inc_slip =
551  pinfo->_incremental_slip -
552  (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
553 
554  // Magnitude of tangential predictor force
555  const Real tan_mag(contact_force_tangential.size());
556  const Real tangential_inc_slip_mag = tangential_inc_slip.size();
557  const Real slip_tol = capacity / penalty;
558  pinfo->_slip_tol = slip_tol;
559 
560  if ((tangential_inc_slip_mag > slip_tol || tan_mag > capacity) &&
561  (pinfo->_stick_locked_this_step < _stick_lock_iterations ||
562  tan_mag > capacity * _stick_unlock_factor))
563  {
564  if (pinfo->_stick_locked_this_step >= _stick_lock_iterations)
565  pinfo->_stick_locked_this_step = 0;
566 
567  // Check for scenario where node has slipped too far in a previous iteration
568  // for special treatment (only done if the incremental slip is non-trivial)
569  bool slipped_too_far = false;
570  RealVectorValue slip_inc_direction;
571  if (tangential_inc_slip_mag > slip_tol)
572  {
573  slip_inc_direction = tangential_inc_slip / tangential_inc_slip_mag;
574  Real slip_dot_tang_force = slip_inc_direction * contact_force_tangential;
575  if (slip_dot_tang_force < capacity)
576  slipped_too_far = true;
577  }
578 
579  if (slipped_too_far) // slip back along slip increment
580  pinfo->_contact_force = contact_force_normal + capacity * slip_inc_direction;
581  else
582  {
583  if (tan_mag > 0) // slip along tangential force direction
584  pinfo->_contact_force =
585  contact_force_normal + capacity * contact_force_tangential / tan_mag;
586  else // treat as frictionless
587  pinfo->_contact_force = contact_force_normal;
588  }
589  if (capacity == 0)
590  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
591  else
592  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
593  }
594  else
595  {
596  if (pinfo->_mech_status != PenetrationInfo::MS_STICKING &&
597  pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
598  ++pinfo->_stick_locked_this_step;
599  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
600  }
601  break;
602  }
603 
604  case CF_PENALTY:
605  {
606  distance_vec = pinfo->_incremental_slip +
607  (pinfo->_normal * (_mesh.nodeRef(node->id()) - pinfo->_closest_point)) *
608  pinfo->_normal;
609  pen_force = penalty * distance_vec;
610 
611  // Frictional capacity
612  // const Real capacity( _friction_coefficient * (pen_force * pinfo->_normal < 0 ?
613  // -pen_force * pinfo->_normal : 0) );
614  const Real capacity(_friction_coefficient *
615  (res_vec * pinfo->_normal > 0 ? res_vec * pinfo->_normal : 0));
616 
617  // Elastic predictor
618  pinfo->_contact_force =
619  pen_force + (pinfo->_contact_force_old -
620  pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old));
621  RealVectorValue contact_force_normal((pinfo->_contact_force * pinfo->_normal) *
622  pinfo->_normal);
623  RealVectorValue contact_force_tangential(pinfo->_contact_force - contact_force_normal);
624 
625  // Tangential magnitude of elastic predictor
626  const Real tan_mag(contact_force_tangential.norm());
627 
628  if (tan_mag > capacity)
629  {
630  pinfo->_contact_force =
631  contact_force_normal + capacity * contact_force_tangential / tan_mag;
632  if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
633  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
634  else
635  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
636  }
637  else
638  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
639  break;
640  }
641 
643  {
644  distance_vec = (pinfo->_normal * (_mesh.nodeRef(node->id()) - pinfo->_closest_point)) *
645  pinfo->_normal;
646 
647  RealVectorValue contact_force_normal =
648  penalty * distance_vec + pinfo->_lagrange_multiplier * pinfo->_normal;
649 
650  RealVectorValue tangential_inc_slip =
651  pinfo->_incremental_slip -
652  (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
653 
654  RealVectorValue contact_force_tangential =
655  pinfo->_lagrange_multiplier_slip +
656  (pinfo->_contact_force_old -
657  pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old));
658 
659  RealVectorValue inc_pen_force_tangential = penalty_slip * tangential_inc_slip;
660 
661  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
662  pinfo->_contact_force =
663  contact_force_normal + contact_force_tangential + inc_pen_force_tangential;
664  else
665  pinfo->_contact_force = contact_force_normal + contact_force_tangential;
666 
667  break;
668  }
669 
671  {
672  // Frictional capacity (kinematic formulation)
673  const Real capacity = _friction_coefficient * std::max(res_vec * pinfo->_normal, 0.0);
674 
675  // Normal component of contact force (kinematic formulation)
676  RealVectorValue contact_force_normal((-res_vec * pinfo->_normal) * pinfo->_normal);
677 
678  // Predictor tangential component of contact force (penalty formulation)
679  RealVectorValue inc_pen_force_tangential = penalty * pinfo->_incremental_slip;
680  RealVectorValue contact_force_tangential =
681  inc_pen_force_tangential +
682  (pinfo->_contact_force_old -
683  pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old));
684 
685  // Magnitude of tangential predictor
686  const Real tan_mag(contact_force_tangential.norm());
687 
688  if (tan_mag > capacity)
689  {
690  pinfo->_contact_force =
691  contact_force_normal + capacity * contact_force_tangential / tan_mag;
692  if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
693  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
694  else
695  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
696  }
697  else
698  {
699  pinfo->_contact_force = contact_force_normal + contact_force_tangential;
700  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
701  }
702  break;
703  }
704 
705  default:
706  mooseError("Invalid contact formulation");
707  break;
708  }
709  break;
710 
711  case CM_GLUED:
712  switch (_formulation)
713  {
714  case CF_KINEMATIC:
715  pinfo->_contact_force = -res_vec;
716  break;
717 
718  case CF_PENALTY:
719  pinfo->_contact_force = pen_force;
720  break;
721 
723  pinfo->_contact_force =
724  pen_force + pinfo->_lagrange_multiplier * distance_vec / distance_vec.norm();
725  break;
726 
727  default:
728  mooseError("Invalid contact formulation");
729  break;
730  }
731  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
732  break;
733 
734  default:
735  mooseError("Invalid or unavailable contact model");
736  break;
737  }
738 }
739 
740 Real
742 {
743  return _u_slave[_qp];
744 }
745 
746 Real
748 {
749  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
750  Real resid = pinfo->_contact_force(_component);
751  switch (type)
752  {
753  case Moose::Slave:
754  if (_formulation == CF_KINEMATIC)
755  {
756  RealVectorValue distance_vec(*_current_node - pinfo->_closest_point);
757  const Real penalty = getPenalty(*pinfo);
758  RealVectorValue pen_force(penalty * distance_vec);
759 
760  if (_model == CM_FRICTIONLESS)
761  resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
762 
763  else if (_model == CM_COULOMB)
764  {
765  distance_vec = distance_vec - pinfo->_incremental_slip;
766  pen_force = penalty * distance_vec;
767  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
768  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
769  resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
770  else
771  resid += pen_force(_component);
772  }
773  else if (_model == CM_GLUED)
774  resid += pen_force(_component);
775  }
777  {
778  RealVectorValue distance_vec(*_current_node - pinfo->_closest_point);
779  const Real penalty = getPenalty(*pinfo);
780  RealVectorValue pen_force(penalty * distance_vec);
781  resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
782  }
783  return _test_slave[_i][_qp] * resid;
784 
785  case Moose::Master:
786  return _test_master[_i][_qp] * -resid;
787  }
788 
789  return 0.0;
790 }
791 
792 Real
793 MechanicalContactConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
794 {
795  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
796 
797  const Real penalty = getPenalty(*pinfo);
798  const Real penalty_slip = getTangentialPenalty(*pinfo);
799 
800  switch (type)
801  {
802  case Moose::SlaveSlave:
803  switch (_model)
804  {
805  case CM_FRICTIONLESS:
806  switch (_formulation)
807  {
808  case CF_KINEMATIC:
809  {
810  RealVectorValue jac_vec;
811  for (unsigned int i = 0; i < _mesh_dimension; ++i)
812  {
813  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
814  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
815  }
816  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
817  (_phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp]) *
818  pinfo->_normal(_component) * pinfo->_normal(_component);
819  }
820 
821  case CF_PENALTY:
823  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
824  pinfo->_normal(_component) * pinfo->_normal(_component);
825 
826  default:
827  mooseError("Invalid contact formulation");
828  }
829 
830  case CM_COULOMB:
831  switch (_formulation)
832  {
833  case CF_KINEMATIC:
834  {
835  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
836  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
837  {
838  RealVectorValue jac_vec;
839  for (unsigned int i = 0; i < _mesh_dimension; ++i)
840  {
841  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
842  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
843  }
844  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
845  (_phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp]) *
846  pinfo->_normal(_component) * pinfo->_normal(_component);
847  }
848  else
849  {
850  const Real curr_jac = (*_jacobian)(
851  _current_node->dof_number(0, _vars[_component], 0), _connected_dof_indices[_j]);
852  return -curr_jac + _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp];
853  }
854  }
855 
856  case CF_PENALTY:
857  {
858  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
859  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
860  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
861  pinfo->_normal(_component) * pinfo->_normal(_component);
862  else
863  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp];
864  }
866  {
867  Real normal_comp = _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
868  pinfo->_normal(_component) * pinfo->_normal(_component);
869 
870  Real tang_comp = 0.0;
871  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
872  tang_comp = _phi_slave[_j][_qp] * penalty_slip * _test_slave[_i][_qp] *
873  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
874  return normal_comp + tang_comp;
875  }
876 
878  {
879  RealVectorValue jac_vec;
880  for (unsigned int i = 0; i < _mesh_dimension; ++i)
881  {
882  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
883  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
884  }
885  Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
886  (_phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp]) *
887  pinfo->_normal(_component) * pinfo->_normal(_component);
888 
889  Real tang_comp = 0.0;
890  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
891  tang_comp = _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
892  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
893 
894  return normal_comp + tang_comp;
895  }
896 
897  default:
898  mooseError("Invalid contact formulation");
899  }
900 
901  case CM_GLUED:
902  switch (_formulation)
903  {
904  case CF_KINEMATIC:
905  {
906  const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
907  _connected_dof_indices[_j]);
908  return -curr_jac + _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp];
909  }
910 
911  case CF_PENALTY:
913  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp];
914 
915  default:
916  mooseError("Invalid contact formulation");
917  }
918  default:
919  mooseError("Invalid or unavailable contact model");
920  }
921 
922  case Moose::SlaveMaster:
923  switch (_model)
924  {
925  case CM_FRICTIONLESS:
926  switch (_formulation)
927  {
928  case CF_KINEMATIC:
929  {
930  Node * curr_master_node = _current_master->get_node(_j);
931 
932  RealVectorValue jac_vec;
933  for (unsigned int i = 0; i < _mesh_dimension; ++i)
934  {
935  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
936  jac_vec(i) =
937  (*_jacobian)(dof_number, curr_master_node->dof_number(0, _vars[_component], 0));
938  }
939  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
940  (_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp]) *
941  pinfo->_normal(_component) * pinfo->_normal(_component);
942  }
943 
944  case CF_PENALTY:
946  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
947  pinfo->_normal(_component) * pinfo->_normal(_component);
948 
949  default:
950  mooseError("Invalid contact formulation");
951  }
952 
953  case CM_COULOMB:
954  switch (_formulation)
955  {
956  case CF_KINEMATIC:
957  {
958  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
959  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
960  {
961  Node * curr_master_node = _current_master->get_node(_j);
962 
963  RealVectorValue jac_vec;
964  for (unsigned int i = 0; i < _mesh_dimension; ++i)
965  {
966  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
967  jac_vec(i) = (*_jacobian)(dof_number,
968  curr_master_node->dof_number(0, _vars[_component], 0));
969  }
970  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
971  (_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp]) *
972  pinfo->_normal(_component) * pinfo->_normal(_component);
973  }
974  else
975  {
976  Node * curr_master_node = _current_master->get_node(_j);
977  const Real curr_jac =
978  (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
979  curr_master_node->dof_number(0, _vars[_component], 0));
980  return -curr_jac - _phi_master[_j][_qp] * penalty * _test_slave[_i][_qp];
981  }
982  }
983 
984  case CF_PENALTY:
985  {
986  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
987  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
988  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
989  pinfo->_normal(_component) * pinfo->_normal(_component);
990  else
991  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp];
992  }
994  {
995  Real normal_comp = -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
996  pinfo->_normal(_component) * pinfo->_normal(_component);
997 
998  Real tang_comp = 0.0;
999  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1000  tang_comp = -_phi_master[_j][_qp] * penalty_slip * _test_slave[_i][_qp] *
1001  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1002  return normal_comp + tang_comp;
1003  }
1004 
1005  case CF_TANGENTIAL_PENALTY:
1006  {
1007  Node * curr_master_node = _current_master->get_node(_j);
1008 
1009  RealVectorValue jac_vec;
1010  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1011  {
1012  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1013  jac_vec(i) =
1014  (*_jacobian)(dof_number, curr_master_node->dof_number(0, _vars[_component], 0));
1015  }
1016  Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1017  (_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp]) *
1018  pinfo->_normal(_component) * pinfo->_normal(_component);
1019 
1020  Real tang_comp = 0.0;
1021  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1022  tang_comp = -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
1023  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1024 
1025  return normal_comp + tang_comp;
1026  }
1027 
1028  default:
1029  mooseError("Invalid contact formulation");
1030  }
1031  case CM_GLUED:
1032  switch (_formulation)
1033  {
1034  case CF_KINEMATIC:
1035  {
1036  Node * curr_master_node = _current_master->get_node(_j);
1037  const Real curr_jac =
1038  (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1039  curr_master_node->dof_number(0, _vars[_component], 0));
1040  return -curr_jac - _phi_master[_j][_qp] * penalty * _test_slave[_i][_qp];
1041  }
1042 
1043  case CF_PENALTY:
1044  case CF_AUGMENTED_LAGRANGE:
1045  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp];
1046 
1047  default:
1048  mooseError("Invalid contact formulation");
1049  }
1050 
1051  default:
1052  mooseError("Invalid or unavailable contact model");
1053  }
1054 
1055  case Moose::MasterSlave:
1056  switch (_model)
1057  {
1058  case CM_FRICTIONLESS:
1059  switch (_formulation)
1060  {
1061  case CF_KINEMATIC:
1062  {
1063  RealVectorValue jac_vec;
1064  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1065  {
1066  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1067  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1068  }
1069  return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1070  _test_master[_i][_qp];
1071  }
1072 
1073  case CF_PENALTY:
1074  case CF_AUGMENTED_LAGRANGE:
1075  return -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp] *
1076  pinfo->_normal(_component) * pinfo->_normal(_component);
1077 
1078  default:
1079  mooseError("Invalid contact formulation");
1080  }
1081  case CM_COULOMB:
1082  switch (_formulation)
1083  {
1084  case CF_KINEMATIC:
1085  {
1086  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1087  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1088  {
1089  RealVectorValue jac_vec;
1090  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1091  {
1092  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1093  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1094  }
1095  return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1096  _test_master[_i][_qp];
1097  }
1098  else
1099  {
1100  const Real slave_jac = (*_jacobian)(
1101  _current_node->dof_number(0, _vars[_component], 0), _connected_dof_indices[_j]);
1102  return slave_jac * _test_master[_i][_qp];
1103  }
1104  }
1105 
1106  case CF_PENALTY:
1107  {
1108  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1109  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1110  return -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp] *
1111  pinfo->_normal(_component) * pinfo->_normal(_component);
1112  else
1113  return -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp];
1114  }
1115  case CF_AUGMENTED_LAGRANGE:
1116  {
1117  Real normal_comp = -_phi_slave[_j][_qp] * penalty * _test_master[_i][_qp] *
1118  pinfo->_normal(_component) * pinfo->_normal(_component);
1119 
1120  Real tang_comp = 0.0;
1121  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1122  tang_comp = -_phi_slave[_j][_qp] * penalty_slip * _test_master[_i][_qp] *
1123  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1124  return normal_comp + tang_comp;
1125  }
1126 
1127  case CF_TANGENTIAL_PENALTY:
1128  {
1129  RealVectorValue jac_vec;
1130  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1131  {
1132  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1133  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1134  }
1135  Real normal_comp =
1136  pinfo->_normal(_component) * (pinfo->_normal * jac_vec) * _test_master[_i][_qp];
1137 
1138  Real tang_comp = 0.0;
1139  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1140  tang_comp = -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp] *
1141  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1142 
1143  return normal_comp + tang_comp;
1144  }
1145 
1146  default:
1147  mooseError("Invalid contact formulation");
1148  }
1149 
1150  case CM_GLUED:
1151  switch (_formulation)
1152  {
1153  case CF_KINEMATIC:
1154  {
1155  const Real slave_jac = (*_jacobian)(
1156  _current_node->dof_number(0, _vars[_component], 0), _connected_dof_indices[_j]);
1157  return slave_jac * _test_master[_i][_qp];
1158  }
1159 
1160  case CF_PENALTY:
1161  case CF_AUGMENTED_LAGRANGE:
1162  return -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp];
1163 
1164  default:
1165  mooseError("Invalid contact formulation");
1166  }
1167 
1168  default:
1169  mooseError("Invalid or unavailable contact model");
1170  }
1171 
1172  case Moose::MasterMaster:
1173  switch (_model)
1174  {
1175  case CM_FRICTIONLESS:
1176  switch (_formulation)
1177  {
1178  case CF_KINEMATIC:
1179  return 0.0;
1180 
1181  case CF_PENALTY:
1182  case CF_AUGMENTED_LAGRANGE:
1183  return _test_master[_i][_qp] * penalty * _phi_master[_j][_qp] *
1184  pinfo->_normal(_component) * pinfo->_normal(_component);
1185 
1186  default:
1187  mooseError("Invalid contact formulation");
1188  }
1189 
1190  case CM_COULOMB:
1191  case CM_GLUED:
1192  switch (_formulation)
1193  {
1194  case CF_KINEMATIC:
1195  return 0.0;
1196 
1197  case CF_PENALTY:
1198  {
1199  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1200  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1201  return _test_master[_i][_qp] * penalty * _phi_master[_j][_qp] *
1202  pinfo->_normal(_component) * pinfo->_normal(_component);
1203  else
1204  return _test_master[_i][_qp] * penalty * _phi_master[_j][_qp];
1205  }
1206 
1207  case CF_TANGENTIAL_PENALTY:
1208  {
1209  Real tang_comp = 0.0;
1210  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1211  tang_comp = _test_master[_i][_qp] * penalty * _phi_master[_j][_qp] *
1212  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1213  return tang_comp; // normal component is zero
1214  }
1215 
1216  case CF_AUGMENTED_LAGRANGE:
1217  {
1218  Real normal_comp = _phi_master[_j][_qp] * penalty * _test_master[_i][_qp] *
1219  pinfo->_normal(_component) * pinfo->_normal(_component);
1220 
1221  Real tang_comp = 0.0;
1222  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1223  tang_comp = _phi_master[_j][_qp] * penalty_slip * _test_master[_i][_qp] *
1224  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1225  return normal_comp + tang_comp;
1226  }
1227 
1228  default:
1229  mooseError("Invalid contact formulation");
1230  }
1231 
1232  default:
1233  mooseError("Invalid or unavailable contact model");
1234  }
1235  }
1236 
1237  return 0.0;
1238 }
1239 
1240 Real
1242  unsigned int jvar)
1243 {
1244  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
1245 
1246  const Real penalty = getPenalty(*pinfo);
1247  const Real penalty_slip = getTangentialPenalty(*pinfo);
1248 
1249  unsigned int coupled_component;
1250  Real normal_component_in_coupled_var_dir = 1.0;
1251  if (getCoupledVarComponent(jvar, coupled_component))
1252  normal_component_in_coupled_var_dir = pinfo->_normal(coupled_component);
1253 
1254  switch (type)
1255  {
1256  case Moose::SlaveSlave:
1257  switch (_model)
1258  {
1259  case CM_FRICTIONLESS:
1260  switch (_formulation)
1261  {
1262  case CF_KINEMATIC:
1263  {
1264  RealVectorValue jac_vec;
1265  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1266  {
1267  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1268  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1269  }
1270  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1271  (_phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp]) *
1272  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1273  }
1274 
1275  case CF_PENALTY:
1276  case CF_AUGMENTED_LAGRANGE:
1277  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
1278  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1279 
1280  default:
1281  mooseError("Invalid contact formulation");
1282  }
1283 
1284  case CM_COULOMB:
1285  {
1287  (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1288  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1289  {
1290  RealVectorValue jac_vec;
1291  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1292  {
1293  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1294  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1295  }
1296  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1297  (_phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp]) *
1298  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1299  }
1300  else if ((_formulation == CF_PENALTY) &&
1301  (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1302  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1303  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
1304  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1305 
1306  else if (_formulation == CF_AUGMENTED_LAGRANGE)
1307  {
1308  Real normal_comp = _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
1309  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1310 
1311  Real tang_comp = 0.0;
1312  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1313  tang_comp = _phi_slave[_j][_qp] * penalty_slip * _test_slave[_i][_qp] *
1314  (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1315  return normal_comp + tang_comp;
1316  }
1317  else
1318  {
1319  const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1320  _connected_dof_indices[_j]);
1321  return -curr_jac;
1322  }
1323  }
1324 
1325  case CM_GLUED:
1326  {
1327  const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1328  _connected_dof_indices[_j]);
1329  return -curr_jac;
1330  }
1331  default:
1332  mooseError("Invalid or unavailable contact model");
1333  }
1334 
1335  case Moose::SlaveMaster:
1336  switch (_model)
1337  {
1338  case CM_FRICTIONLESS:
1339  switch (_formulation)
1340  {
1341  case CF_KINEMATIC:
1342  {
1343  Node * curr_master_node = _current_master->get_node(_j);
1344 
1345  RealVectorValue jac_vec;
1346  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1347  {
1348  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1349  jac_vec(i) =
1350  (*_jacobian)(dof_number, curr_master_node->dof_number(0, _vars[_component], 0));
1351  }
1352  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1353  (_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp]) *
1354  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1355  }
1356 
1357  case CF_PENALTY:
1358  case CF_AUGMENTED_LAGRANGE:
1359  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
1360  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1361 
1362  default:
1363  mooseError("Invalid contact formulation");
1364  }
1365 
1366  case CM_COULOMB:
1368  (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1369  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1370  {
1371  Node * curr_master_node = _current_master->get_node(_j);
1372 
1373  RealVectorValue jac_vec;
1374  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1375  {
1376  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1377  jac_vec(i) =
1378  (*_jacobian)(dof_number, curr_master_node->dof_number(0, _vars[_component], 0));
1379  }
1380  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1381  (_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp]) *
1382  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1383  }
1384  else if ((_formulation == CF_PENALTY) &&
1385  (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1386  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1387  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
1388  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1389 
1390  else if (_formulation == CF_AUGMENTED_LAGRANGE)
1391  {
1392  Real normal_comp = -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
1393  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1394 
1395  Real tang_comp = 0.0;
1396  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1397  tang_comp = -_phi_master[_j][_qp] * penalty_slip * _test_slave[_i][_qp] *
1398  (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1399  return normal_comp + tang_comp;
1400  }
1401  else
1402  return 0.0;
1403 
1404  case CM_GLUED:
1405  return 0;
1406 
1407  default:
1408  mooseError("Invalid or unavailable contact model");
1409  }
1410 
1411  case Moose::MasterSlave:
1412  switch (_model)
1413  {
1414  case CM_FRICTIONLESS:
1415  switch (_formulation)
1416  {
1417  case CF_KINEMATIC:
1418  {
1419  RealVectorValue jac_vec;
1420  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1421  {
1422  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1423  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1424  }
1425  return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1426  _test_master[_i][_qp];
1427  }
1428 
1429  case CF_PENALTY:
1430  case CF_AUGMENTED_LAGRANGE:
1431  return -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp] *
1432  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1433 
1434  default:
1435  mooseError("Invalid contact formulation");
1436  }
1437  case CM_COULOMB:
1438  switch (_formulation)
1439  {
1440  case CF_KINEMATIC:
1441  {
1442  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1443  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1444  {
1445  RealVectorValue jac_vec;
1446  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1447  {
1448  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1449  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1450  }
1451  return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1452  _test_master[_i][_qp];
1453  }
1454  else
1455  {
1456  const Real slave_jac = (*_jacobian)(
1457  _current_node->dof_number(0, _vars[_component], 0), _connected_dof_indices[_j]);
1458  return slave_jac * _test_master[_i][_qp];
1459  }
1460  }
1461 
1462  case CF_PENALTY:
1463  {
1464  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1465  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1466  return -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp] *
1467  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1468  else
1469  return 0.0;
1470  }
1471  case CF_AUGMENTED_LAGRANGE:
1472  {
1473  Real normal_comp = -_phi_slave[_j][_qp] * penalty * _test_master[_i][_qp] *
1474  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1475 
1476  Real tang_comp = 0.0;
1477  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1478  tang_comp = -_phi_slave[_j][_qp] * penalty_slip * _test_master[_i][_qp] *
1479  (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1480  return normal_comp + tang_comp;
1481  }
1482  case CF_TANGENTIAL_PENALTY:
1483  {
1484  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1485  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1486  {
1487  RealVectorValue jac_vec;
1488  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1489  {
1490  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1491  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1492  }
1493  return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1494  _test_master[_i][_qp];
1495  }
1496  else
1497  return 0.0;
1498  }
1499 
1500  default:
1501  mooseError("Invalid contact formulation");
1502  }
1503 
1504  case CM_GLUED:
1505  switch (_formulation)
1506  {
1507  case CF_KINEMATIC:
1508  {
1509  const Real slave_jac = (*_jacobian)(
1510  _current_node->dof_number(0, _vars[_component], 0), _connected_dof_indices[_j]);
1511  return slave_jac * _test_master[_i][_qp];
1512  }
1513 
1514  case CF_PENALTY:
1515  case CF_AUGMENTED_LAGRANGE:
1516  return 0.0;
1517 
1518  default:
1519  mooseError("Invalid contact formulation");
1520  }
1521 
1522  default:
1523  mooseError("Invalid or unavailable contact model");
1524  }
1525 
1526  case Moose::MasterMaster:
1527  switch (_model)
1528  {
1529  case CM_FRICTIONLESS:
1530  switch (_formulation)
1531  {
1532  case CF_KINEMATIC:
1533  return 0.0;
1534 
1535  case CF_PENALTY:
1536  case CF_AUGMENTED_LAGRANGE:
1537  return _test_master[_i][_qp] * penalty * _phi_master[_j][_qp] *
1538  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1539 
1540  default:
1541  mooseError("Invalid contact formulation");
1542  }
1543 
1544  case CM_COULOMB:
1545  {
1547  {
1548 
1549  Real normal_comp = _phi_master[_j][_qp] * penalty * _test_master[_i][_qp] *
1550  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1551 
1552  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1553  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1554  return normal_comp;
1555  }
1556  else if (_formulation == CF_PENALTY &&
1557  (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1558  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1559  return _test_master[_i][_qp] * penalty * _phi_master[_j][_qp] *
1560  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1561  else
1562  return 0.0;
1563  }
1564 
1565  case CM_GLUED:
1566  if (_formulation == CF_PENALTY &&
1567  (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1568  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1569  return _test_master[_i][_qp] * penalty * _phi_master[_j][_qp] *
1570  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1571  else if (_formulation == CF_AUGMENTED_LAGRANGE)
1572  {
1573  Real normal_comp = _phi_master[_j][_qp] * penalty * _test_master[_i][_qp] *
1574  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1575 
1576  Real tang_comp = 0.0;
1577  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1578  tang_comp = _phi_master[_j][_qp] * penalty_slip * _test_master[_i][_qp] *
1579  (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1580  return normal_comp + tang_comp;
1581  }
1582  else
1583  return 0.0;
1584 
1585  default:
1586  mooseError("Invalid or unavailable contact model");
1587  }
1588  }
1589 
1590  return 0.0;
1591 }
1592 
1593 Real
1595 {
1596  const Node * node = pinfo._node;
1597 
1598  dof_id_type dof = node->dof_number(_aux_system.number(), _nodal_area_var->number(), 0);
1599 
1600  Real area = (*_aux_solution)(dof);
1601  if (area == 0.0)
1602  {
1603  if (_t_step > 1)
1604  mooseError("Zero nodal area found");
1605  else
1606  area = 1.0; // Avoid divide by zero during initialization
1607  }
1608  return area;
1609 }
1610 
1611 Real
1613 {
1614  Real penalty = _penalty;
1615  if (_normalize_penalty)
1616  penalty *= nodalArea(pinfo);
1617 
1618  return penalty;
1619 }
1620 
1621 Real
1623 {
1624  Real penalty = _penalty_tangential;
1625  if (_normalize_penalty)
1626  penalty *= nodalArea(pinfo);
1627 
1628  return penalty;
1629 }
1630 
1631 void
1633 {
1634  getConnectedDofIndices(_var.number());
1635 
1636  DenseMatrix<Number> & Knn =
1637  _assembly.jacobianBlockNeighbor(Moose::NeighborNeighbor, _master_var.number(), _var.number());
1638 
1639  _Kee.resize(_test_slave.size(), _connected_dof_indices.size());
1640 
1641  for (_i = 0; _i < _test_slave.size(); _i++)
1642  // Loop over the connected dof indices so we can get all the jacobian contributions
1643  for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1644  _Kee(_i, _j) += computeQpJacobian(Moose::SlaveSlave);
1645 
1647  {
1648  DenseMatrix<Number> & Ken =
1649  _assembly.jacobianBlockNeighbor(Moose::ElementNeighbor, _var.number(), _var.number());
1650  if (Ken.m() && Ken.n())
1651  for (_i = 0; _i < _test_slave.size(); _i++)
1652  for (_j = 0; _j < _phi_master.size(); _j++)
1653  Ken(_i, _j) += computeQpJacobian(Moose::SlaveMaster);
1654 
1655  _Kne.resize(_test_master.size(), _connected_dof_indices.size());
1656  for (_i = 0; _i < _test_master.size(); _i++)
1657  // Loop over the connected dof indices so we can get all the jacobian contributions
1658  for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1659  _Kne(_i, _j) += computeQpJacobian(Moose::MasterSlave);
1660  }
1661 
1662  if (Knn.m() && Knn.n())
1663  for (_i = 0; _i < _test_master.size(); _i++)
1664  for (_j = 0; _j < _phi_master.size(); _j++)
1665  Knn(_i, _j) += computeQpJacobian(Moose::MasterMaster);
1666 }
1667 
1668 void
1670 {
1671  getConnectedDofIndices(jvar);
1672 
1673  _Kee.resize(_test_slave.size(), _connected_dof_indices.size());
1674 
1675  DenseMatrix<Number> & Knn =
1676  _assembly.jacobianBlockNeighbor(Moose::NeighborNeighbor, _master_var.number(), jvar);
1677 
1678  for (_i = 0; _i < _test_slave.size(); _i++)
1679  // Loop over the connected dof indices so we can get all the jacobian contributions
1680  for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1681  _Kee(_i, _j) += computeQpOffDiagJacobian(Moose::SlaveSlave, jvar);
1682 
1684  {
1685  DenseMatrix<Number> & Ken =
1686  _assembly.jacobianBlockNeighbor(Moose::ElementNeighbor, _var.number(), jvar);
1687  for (_i = 0; _i < _test_slave.size(); _i++)
1688  for (_j = 0; _j < _phi_master.size(); _j++)
1689  Ken(_i, _j) += computeQpOffDiagJacobian(Moose::SlaveMaster, jvar);
1690 
1691  _Kne.resize(_test_master.size(), _connected_dof_indices.size());
1692  if (_Kne.m() && _Kne.n())
1693  for (_i = 0; _i < _test_master.size(); _i++)
1694  // Loop over the connected dof indices so we can get all the jacobian contributions
1695  for (_j = 0; _j < _connected_dof_indices.size(); _j++)
1696  _Kne(_i, _j) += computeQpOffDiagJacobian(Moose::MasterSlave, jvar);
1697  }
1698 
1699  for (_i = 0; _i < _test_master.size(); _i++)
1700  for (_j = 0; _j < _phi_master.size(); _j++)
1701  Knn(_i, _j) += computeQpOffDiagJacobian(Moose::MasterMaster, jvar);
1702 }
1703 
1704 void
1706 {
1707  unsigned int component;
1708  if (getCoupledVarComponent(var_num, component) || _non_displacement_vars_jacobian)
1709  {
1711  NodeFaceConstraint::getConnectedDofIndices(var_num);
1712  else
1713  {
1714  _connected_dof_indices.clear();
1715  MooseVariable & var = _sys.getVariable(0, var_num);
1716  _connected_dof_indices.push_back(var.nodalDofIndex());
1717  }
1718  }
1719 
1720  _phi_slave.resize(_connected_dof_indices.size());
1721 
1722  dof_id_type current_node_var_dof_index = _sys.getVariable(0, var_num).nodalDofIndex();
1723  _qp = 0;
1724 
1725  // Fill up _phi_slave so that it is 1 when j corresponds to the dof associated with this node
1726  // and 0 for every other dof
1727  // This corresponds to evaluating all of the connected shape functions at _this_ node
1728  for (unsigned int j = 0; j < _connected_dof_indices.size(); j++)
1729  {
1730  _phi_slave[j].resize(1);
1731 
1732  if (_connected_dof_indices[j] == current_node_var_dof_index)
1733  _phi_slave[j][_qp] = 1.0;
1734  else
1735  _phi_slave[j][_qp] = 0.0;
1736  }
1737 }
1738 
1739 bool
1741 {
1742  component = std::numeric_limits<unsigned int>::max();
1743  bool coupled_var_is_disp_var = false;
1744  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1745  {
1746  if (var_num == _vars[i])
1747  {
1748  coupled_var_is_disp_var = true;
1749  component = i;
1750  break;
1751  }
1752  }
1753  return coupled_var_is_disp_var;
1754 }
NumericVector< Number > & _residual_copy
Real nodalArea(PenetrationInfo &pinfo)
ContactFormulation contactFormulation(const std::string &the_name)
bool getCoupledVarComponent(unsigned int var_num, unsigned int &component)
Determine whether the coupled variable is one of the displacement variables, and find its component...
const ContactFormulation _formulation
virtual void getConnectedDofIndices(unsigned int var_num) override
Get the dof indices of the nodes connected to the slave node for a specific variable.
virtual Real computeQpOffDiagJacobian(Moose::ConstraintJacobianType type, unsigned int jvar) override
Compute off-diagonal Jacobian entries.
virtual void updateAugmentedLagrangianMultiplier(bool beginning_of_step=false)
Real component(const SymmTensor &symm_tensor, unsigned int index)
Real _al_penetration_tolerance
The tolerance of the penetration for augmented Lagrangian method.
virtual Real computeQpResidual(Moose::ConstraintType type) override
void computeContactForce(PenetrationInfo *pinfo)
Real getPenalty(PenetrationInfo &pinfo)
virtual void updateContactSet(bool beginning_of_step=false)
virtual Real computeQpSlaveValue() override
virtual Real computeQpJacobian(Moose::ConstraintJacobianType type) override
virtual void computeOffDiagJacobian(unsigned int jvar) override
Compute off-diagonal Jacobian entries.
const bool _non_displacement_vars_jacobian
Whether to include coupling terms with non-displacement variables in the Jacobian.
const bool _connected_slave_nodes_jacobian
Whether to include coupling terms with the nodes connected to the slave nodes in the Jacobian...
Real getTangentialPenalty(PenetrationInfo &pinfo)
std::vector< unsigned int > _vars
ContactModel contactModel(const std::string &the_name)
virtual void computeJacobian() override
Computes the jacobian for the current element.
Real _al_incremental_slip_tolerance
The tolerance of the incremental slip for augmented Lagrangian method.
MechanicalContactConstraint(const InputParameters &parameters)
Real _al_frictional_force_tolerance
The tolerance of the frictional force for augmented Lagrangian method.
const bool _master_slave_jacobian
Whether to include coupling between the master and slave nodes in the Jacobian.
InputParameters validParams< MechanicalContactConstraint >()