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