www.mooseframework.org
Public Member Functions | Protected Member Functions | Protected Attributes | Static Protected Attributes | List of all members
MechanicalContactConstraint Class Reference

A MechanicalContactConstraint forces the value of a variable to be the same on both sides of an interface. More...

#include <MechanicalContactConstraint.h>

Inheritance diagram for MechanicalContactConstraint:
[legend]

Public Member Functions

 MechanicalContactConstraint (const InputParameters &parameters)
 
virtual void timestepSetup () override
 
virtual void jacobianSetup () override
 
virtual void residualEnd () override
 
virtual bool AugmentedLagrangianContactConverged ()
 
virtual void updateAugmentedLagrangianMultiplier (bool beginning_of_step=false)
 
virtual void updateContactStatefulData (bool beginning_of_step=false)
 
virtual Real computeQpSlaveValue () override
 
virtual Real computeQpResidual (Moose::ConstraintType type) override
 
virtual void computeJacobian () override
 Computes the jacobian for the current element. More...
 
virtual void computeOffDiagJacobian (unsigned int jvar) override
 Compute off-diagonal Jacobian entries. More...
 
virtual Real computeQpJacobian (Moose::ConstraintJacobianType type) override
 
virtual Real computeQpOffDiagJacobian (Moose::ConstraintJacobianType type, unsigned int jvar) override
 Compute off-diagonal Jacobian entries. More...
 
virtual void getConnectedDofIndices (unsigned int var_num) override
 Get the dof indices of the nodes connected to the slave node for a specific variable. More...
 
bool getCoupledVarComponent (unsigned int var_num, unsigned int &component)
 Determine whether the coupled variable is one of the displacement variables, and find its component. More...
 
virtual bool addCouplingEntriesToJacobian () override
 
bool shouldApply () override
 
void computeContactForce (PenetrationInfo *pinfo, bool update_contact_set)
 

Protected Member Functions

Real nodalArea (PenetrationInfo &pinfo)
 
Real getPenalty (PenetrationInfo &pinfo)
 
Real getTangentialPenalty (PenetrationInfo &pinfo)
 

Protected Attributes

MooseSharedPointer< DisplacedProblem > _displaced_problem
 
const unsigned int _component
 
const ContactModel _model
 
const ContactFormulation _formulation
 
const bool _normalize_penalty
 
const Real _penalty
 
Real _penalty_tangential
 
const Real _friction_coefficient
 
const Real _tension_release
 
const Real _capture_tolerance
 
const unsigned int _stick_lock_iterations
 
const Real _stick_unlock_factor
 
bool _update_stateful_data
 
NumericVector< Number > & _residual_copy
 
const unsigned int _mesh_dimension
 
std::vector< unsigned int > _vars
 
std::vector< MooseVariable * > _var_objects
 
MooseVariable * _nodal_area_var
 
SystemBase & _aux_system
 
const NumericVector< Number > * _aux_solution
 
const bool _master_slave_jacobian
 Whether to include coupling between the master and slave nodes in the Jacobian. More...
 
const bool _connected_slave_nodes_jacobian
 Whether to include coupling terms with the nodes connected to the slave nodes in the Jacobian. More...
 
const bool _non_displacement_vars_jacobian
 Whether to include coupling terms with non-displacement variables in the Jacobian. More...
 
Real _al_penetration_tolerance
 The tolerance of the penetration for augmented Lagrangian method. More...
 
Real _al_incremental_slip_tolerance
 The tolerance of the incremental slip for augmented Lagrangian method. More...
 
Real _al_frictional_force_tolerance
 The tolerance of the frictional force for augmented Lagrangian method. More...
 
ContactLineSearchBase * _contact_linesearch
 
std::set< dof_id_type > _current_contact_state
 
std::set< dof_id_type > _old_contact_state
 
const bool _print_contact_nodes
 

Static Protected Attributes

static Threads::spin_mutex _contact_set_mutex
 

Detailed Description

A MechanicalContactConstraint forces the value of a variable to be the same on both sides of an interface.

Definition at line 29 of file MechanicalContactConstraint.h.

Constructor & Destructor Documentation

◆ MechanicalContactConstraint()

MechanicalContactConstraint::MechanicalContactConstraint ( const InputParameters &  parameters)

Definition at line 110 of file MechanicalContactConstraint.C.

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 }

Member Function Documentation

◆ addCouplingEntriesToJacobian()

virtual bool MechanicalContactConstraint::addCouplingEntriesToJacobian ( )
inlineoverridevirtual

Definition at line 85 of file MechanicalContactConstraint.h.

85 { return _master_slave_jacobian; }

◆ AugmentedLagrangianContactConverged()

bool MechanicalContactConstraint::AugmentedLagrangianContactConverged ( )
virtual

Definition at line 330 of file MechanicalContactConstraint.C.

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 }

◆ computeContactForce()

void MechanicalContactConstraint::computeContactForce ( PenetrationInfo *  pinfo,
bool  update_contact_set 
)

Definition at line 496 of file MechanicalContactConstraint.C.

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 }

Referenced by shouldApply().

◆ computeJacobian()

void MechanicalContactConstraint::computeJacobian ( )
overridevirtual

Computes the jacobian for the current element.

Definition at line 1684 of file MechanicalContactConstraint.C.

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 }

◆ computeOffDiagJacobian()

void MechanicalContactConstraint::computeOffDiagJacobian ( unsigned int  jvar)
overridevirtual

Compute off-diagonal Jacobian entries.

Parameters
jvarThe index of the coupled variable

Definition at line 1721 of file MechanicalContactConstraint.C.

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 }

◆ computeQpJacobian()

Real MechanicalContactConstraint::computeQpJacobian ( Moose::ConstraintJacobianType  type)
overridevirtual

Definition at line 835 of file MechanicalContactConstraint.C.

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 }

Referenced by computeJacobian().

◆ computeQpOffDiagJacobian()

Real MechanicalContactConstraint::computeQpOffDiagJacobian ( Moose::ConstraintJacobianType  type,
unsigned int  jvar 
)
overridevirtual

Compute off-diagonal Jacobian entries.

Parameters
typeThe type of coupling
jvarThe index of the coupled variable

Definition at line 1286 of file MechanicalContactConstraint.C.

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 }

Referenced by computeOffDiagJacobian().

◆ computeQpResidual()

Real MechanicalContactConstraint::computeQpResidual ( Moose::ConstraintType  type)
overridevirtual

Definition at line 788 of file MechanicalContactConstraint.C.

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 }

◆ computeQpSlaveValue()

Real MechanicalContactConstraint::computeQpSlaveValue ( )
overridevirtual

Definition at line 782 of file MechanicalContactConstraint.C.

783 {
784  return _u_slave[_qp];
785 }

◆ getConnectedDofIndices()

void MechanicalContactConstraint::getConnectedDofIndices ( unsigned int  var_num)
overridevirtual

Get the dof indices of the nodes connected to the slave node for a specific variable.

Parameters
var_numThe number of the variable for which dof indices are gathered
Returns
bool indicating whether the coupled variable is one of the displacement variables

Definition at line 1757 of file MechanicalContactConstraint.C.

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 }

Referenced by computeJacobian(), and computeOffDiagJacobian().

◆ getCoupledVarComponent()

bool MechanicalContactConstraint::getCoupledVarComponent ( unsigned int  var_num,
unsigned int &  component 
)

Determine whether the coupled variable is one of the displacement variables, and find its component.

Parameters
var_numThe number of the variable to be checked
componentThe component index computed in this routine
Returns
bool indicating whether the coupled variable is one of the displacement variables

Definition at line 1792 of file MechanicalContactConstraint.C.

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 }

Referenced by computeQpOffDiagJacobian(), and getConnectedDofIndices().

◆ getPenalty()

Real MechanicalContactConstraint::getPenalty ( PenetrationInfo &  pinfo)
protected

Definition at line 1664 of file MechanicalContactConstraint.C.

1665 {
1666  Real penalty = _penalty;
1667  if (_normalize_penalty)
1668  penalty *= nodalArea(pinfo);
1669 
1670  return penalty;
1671 }

Referenced by AugmentedLagrangianContactConverged(), computeContactForce(), computeQpJacobian(), computeQpOffDiagJacobian(), computeQpResidual(), and updateAugmentedLagrangianMultiplier().

◆ getTangentialPenalty()

Real MechanicalContactConstraint::getTangentialPenalty ( PenetrationInfo &  pinfo)
protected

Definition at line 1674 of file MechanicalContactConstraint.C.

1675 {
1676  Real penalty = _penalty_tangential;
1677  if (_normalize_penalty)
1678  penalty *= nodalArea(pinfo);
1679 
1680  return penalty;
1681 }

Referenced by computeContactForce(), computeQpJacobian(), computeQpOffDiagJacobian(), and updateAugmentedLagrangianMultiplier().

◆ jacobianSetup()

void MechanicalContactConstraint::jacobianSetup ( )
overridevirtual

Definition at line 243 of file MechanicalContactConstraint.C.

244 {
245  if (_component == 0)
246  {
247  if (_update_stateful_data)
248  updateContactStatefulData();
249  _update_stateful_data = true;
250  }
251 }

◆ nodalArea()

Real MechanicalContactConstraint::nodalArea ( PenetrationInfo &  pinfo)
protected

Definition at line 1646 of file MechanicalContactConstraint.C.

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 }

Referenced by computeContactForce(), getPenalty(), and getTangentialPenalty().

◆ residualEnd()

void MechanicalContactConstraint::residualEnd ( )
overridevirtual

Definition at line 1809 of file MechanicalContactConstraint.C.

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 }

◆ shouldApply()

bool MechanicalContactConstraint::shouldApply ( )
override

Definition at line 462 of file MechanicalContactConstraint.C.

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 }

◆ timestepSetup()

void MechanicalContactConstraint::timestepSetup ( )
overridevirtual

Definition at line 227 of file MechanicalContactConstraint.C.

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 }

◆ updateAugmentedLagrangianMultiplier()

void MechanicalContactConstraint::updateAugmentedLagrangianMultiplier ( bool  beginning_of_step = false)
virtual

Definition at line 254 of file MechanicalContactConstraint.C.

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 }

Referenced by timestepSetup().

◆ updateContactStatefulData()

void MechanicalContactConstraint::updateContactStatefulData ( bool  beginning_of_step = false)
virtual

Definition at line 423 of file MechanicalContactConstraint.C.

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 }

Referenced by jacobianSetup(), and timestepSetup().

Member Data Documentation

◆ _al_frictional_force_tolerance

Real MechanicalContactConstraint::_al_frictional_force_tolerance
protected

The tolerance of the frictional force for augmented Lagrangian method.

Definition at line 134 of file MechanicalContactConstraint.h.

Referenced by AugmentedLagrangianContactConverged(), MechanicalContactConstraint(), and updateAugmentedLagrangianMultiplier().

◆ _al_incremental_slip_tolerance

Real MechanicalContactConstraint::_al_incremental_slip_tolerance
protected

The tolerance of the incremental slip for augmented Lagrangian method.

Definition at line 132 of file MechanicalContactConstraint.h.

Referenced by AugmentedLagrangianContactConverged(), and MechanicalContactConstraint().

◆ _al_penetration_tolerance

Real MechanicalContactConstraint::_al_penetration_tolerance
protected

The tolerance of the penetration for augmented Lagrangian method.

Definition at line 130 of file MechanicalContactConstraint.h.

Referenced by AugmentedLagrangianContactConverged(), and MechanicalContactConstraint().

◆ _aux_solution

const NumericVector<Number>* MechanicalContactConstraint::_aux_solution
protected

Definition at line 120 of file MechanicalContactConstraint.h.

◆ _aux_system

SystemBase& MechanicalContactConstraint::_aux_system
protected

Definition at line 119 of file MechanicalContactConstraint.h.

Referenced by nodalArea().

◆ _capture_tolerance

const Real MechanicalContactConstraint::_capture_tolerance
protected

Definition at line 105 of file MechanicalContactConstraint.h.

Referenced by computeContactForce().

◆ _component

const unsigned int MechanicalContactConstraint::_component
protected

Definition at line 96 of file MechanicalContactConstraint.h.

Referenced by AugmentedLagrangianContactConverged(), computeQpJacobian(), computeQpOffDiagJacobian(), computeQpResidual(), jacobianSetup(), residualEnd(), shouldApply(), timestepSetup(), updateAugmentedLagrangianMultiplier(), and updateContactStatefulData().

◆ _connected_slave_nodes_jacobian

const bool MechanicalContactConstraint::_connected_slave_nodes_jacobian
protected

Whether to include coupling terms with the nodes connected to the slave nodes in the Jacobian.

Definition at line 125 of file MechanicalContactConstraint.h.

Referenced by getConnectedDofIndices().

◆ _contact_linesearch

ContactLineSearchBase* MechanicalContactConstraint::_contact_linesearch
protected

Definition at line 136 of file MechanicalContactConstraint.h.

Referenced by computeContactForce(), residualEnd(), and timestepSetup().

◆ _contact_set_mutex

Threads::spin_mutex MechanicalContactConstraint::_contact_set_mutex
staticprotected

Definition at line 141 of file MechanicalContactConstraint.h.

Referenced by shouldApply().

◆ _current_contact_state

std::set<dof_id_type> MechanicalContactConstraint::_current_contact_state
protected

Definition at line 137 of file MechanicalContactConstraint.h.

Referenced by residualEnd(), and shouldApply().

◆ _displaced_problem

MooseSharedPointer<DisplacedProblem> MechanicalContactConstraint::_displaced_problem
protected

Definition at line 91 of file MechanicalContactConstraint.h.

◆ _formulation

const ContactFormulation MechanicalContactConstraint::_formulation
protected

Definition at line 98 of file MechanicalContactConstraint.h.

Referenced by computeContactForce(), computeQpJacobian(), computeQpOffDiagJacobian(), computeQpResidual(), MechanicalContactConstraint(), and timestepSetup().

◆ _friction_coefficient

const Real MechanicalContactConstraint::_friction_coefficient
protected

Definition at line 103 of file MechanicalContactConstraint.h.

Referenced by AugmentedLagrangianContactConverged(), computeContactForce(), MechanicalContactConstraint(), and updateAugmentedLagrangianMultiplier().

◆ _master_slave_jacobian

const bool MechanicalContactConstraint::_master_slave_jacobian
protected

Whether to include coupling between the master and slave nodes in the Jacobian.

Definition at line 123 of file MechanicalContactConstraint.h.

Referenced by addCouplingEntriesToJacobian(), computeJacobian(), computeOffDiagJacobian(), and getConnectedDofIndices().

◆ _mesh_dimension

const unsigned int MechanicalContactConstraint::_mesh_dimension
protected

Definition at line 113 of file MechanicalContactConstraint.h.

Referenced by computeContactForce(), computeQpJacobian(), and computeQpOffDiagJacobian().

◆ _model

const ContactModel MechanicalContactConstraint::_model
protected

Definition at line 97 of file MechanicalContactConstraint.h.

Referenced by AugmentedLagrangianContactConverged(), computeContactForce(), computeQpJacobian(), computeQpOffDiagJacobian(), computeQpResidual(), MechanicalContactConstraint(), and updateAugmentedLagrangianMultiplier().

◆ _nodal_area_var

MooseVariable* MechanicalContactConstraint::_nodal_area_var
protected

Definition at line 118 of file MechanicalContactConstraint.h.

Referenced by nodalArea().

◆ _non_displacement_vars_jacobian

const bool MechanicalContactConstraint::_non_displacement_vars_jacobian
protected

Whether to include coupling terms with non-displacement variables in the Jacobian.

Definition at line 127 of file MechanicalContactConstraint.h.

Referenced by getConnectedDofIndices().

◆ _normalize_penalty

const bool MechanicalContactConstraint::_normalize_penalty
protected

Definition at line 99 of file MechanicalContactConstraint.h.

Referenced by getPenalty(), and getTangentialPenalty().

◆ _old_contact_state

std::set<dof_id_type> MechanicalContactConstraint::_old_contact_state
protected

Definition at line 138 of file MechanicalContactConstraint.h.

Referenced by residualEnd().

◆ _penalty

const Real MechanicalContactConstraint::_penalty
protected

Definition at line 101 of file MechanicalContactConstraint.h.

Referenced by getPenalty(), and MechanicalContactConstraint().

◆ _penalty_tangential

Real MechanicalContactConstraint::_penalty_tangential
protected

Definition at line 102 of file MechanicalContactConstraint.h.

Referenced by getTangentialPenalty(), and MechanicalContactConstraint().

◆ _print_contact_nodes

const bool MechanicalContactConstraint::_print_contact_nodes
protected

Definition at line 140 of file MechanicalContactConstraint.h.

Referenced by residualEnd().

◆ _residual_copy

NumericVector<Number>& MechanicalContactConstraint::_residual_copy
protected

Definition at line 110 of file MechanicalContactConstraint.h.

Referenced by computeContactForce().

◆ _stick_lock_iterations

const unsigned int MechanicalContactConstraint::_stick_lock_iterations
protected

Definition at line 106 of file MechanicalContactConstraint.h.

Referenced by computeContactForce().

◆ _stick_unlock_factor

const Real MechanicalContactConstraint::_stick_unlock_factor
protected

Definition at line 107 of file MechanicalContactConstraint.h.

Referenced by computeContactForce().

◆ _tension_release

const Real MechanicalContactConstraint::_tension_release
protected

Definition at line 104 of file MechanicalContactConstraint.h.

Referenced by computeContactForce().

◆ _update_stateful_data

bool MechanicalContactConstraint::_update_stateful_data
protected

Definition at line 108 of file MechanicalContactConstraint.h.

Referenced by jacobianSetup(), and timestepSetup().

◆ _var_objects

std::vector<MooseVariable *> MechanicalContactConstraint::_var_objects
protected

Definition at line 116 of file MechanicalContactConstraint.h.

Referenced by computeContactForce(), and MechanicalContactConstraint().

◆ _vars

std::vector<unsigned int> MechanicalContactConstraint::_vars
protected

Definition at line 115 of file MechanicalContactConstraint.h.

Referenced by AugmentedLagrangianContactConverged(), computeContactForce(), computeQpJacobian(), computeQpOffDiagJacobian(), getCoupledVarComponent(), MechanicalContactConstraint(), updateAugmentedLagrangianMultiplier(), and updateContactStatefulData().

The documentation for this class was generated from the following files:
MechanicalContactConstraint::_aux_solution
const NumericVector< Number > * _aux_solution
Definition: MechanicalContactConstraint.h:120
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::_capture_tolerance
const Real _capture_tolerance
Definition: MechanicalContactConstraint.h:105
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
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
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::_displaced_problem
MooseSharedPointer< DisplacedProblem > _displaced_problem
Definition: MechanicalContactConstraint.h:91
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::_stick_lock_iterations
const unsigned int _stick_lock_iterations
Definition: MechanicalContactConstraint.h:106
ContactModel::FRICTIONLESS
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::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
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
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
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
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
ContactFormulation::PENALTY
MechanicalContactConstraint::_nodal_area_var
MooseVariable * _nodal_area_var
Definition: MechanicalContactConstraint.h:118
MechanicalContactConstraint::_var_objects
std::vector< MooseVariable * > _var_objects
Definition: MechanicalContactConstraint.h:116
Generated on Sat Jan 25 2020 12:05:07 for www.mooseframework.org by   doxygen 1.8.16