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
 
FEProblem & _fe_problem
 
const unsigned int _component
 
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...
 
std::shared_ptr< 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 28 of file MechanicalContactConstraint.h.

Constructor & Destructor Documentation

◆ MechanicalContactConstraint()

MechanicalContactConstraint::MechanicalContactConstraint ( const InputParameters &  parameters)

Definition at line 119 of file MechanicalContactConstraint.C.

120  : NodeFaceConstraint(parameters),
121  _displaced_problem(parameters.get<FEProblemBase *>("_fe_problem_base")->getDisplacedProblem()),
122  _fe_problem(*parameters.get<FEProblem *>("_fe_problem")),
123  _component(getParam<unsigned int>("component")),
124  _model(ContactMaster::contactModel(getParam<std::string>("model"))),
125  _formulation(ContactMaster::contactFormulation(getParam<std::string>("formulation"))),
126  _normalize_penalty(getParam<bool>("normalize_penalty")),
127  _penalty(getParam<Real>("penalty")),
128  _friction_coefficient(getParam<Real>("friction_coefficient")),
129  _tension_release(getParam<Real>("tension_release")),
130  _capture_tolerance(getParam<Real>("capture_tolerance")),
131  _stick_lock_iterations(getParam<unsigned int>("stick_lock_iterations")),
132  _stick_unlock_factor(getParam<Real>("stick_unlock_factor")),
133  _update_stateful_data(true),
134  _residual_copy(_sys.residualGhosted()),
135  _mesh_dimension(_mesh.dimension()),
136  _vars(3, libMesh::invalid_uint),
137  _var_objects(3, nullptr),
138  _nodal_area_var(getVar("nodal_area", 0)),
139  _aux_system(_nodal_area_var->sys()),
140  _aux_solution(_aux_system.currentSolution()),
141  _master_slave_jacobian(getParam<bool>("master_slave_jacobian")),
142  _connected_slave_nodes_jacobian(getParam<bool>("connected_slave_nodes_jacobian")),
143  _non_displacement_vars_jacobian(getParam<bool>("non_displacement_variables_jacobian")),
144  _contact_linesearch(
145  std::dynamic_pointer_cast<ContactLineSearchBase>(_fe_problem.getLineSearch())),
146  _print_contact_nodes(getParam<bool>("print_contact_nodes"))
147 {
148  _overwrite_slave_residual = false;
149 
150  if (isParamValid("displacements"))
151  {
152  // modern parameter scheme for displacements
153  for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
154  {
155  _vars[i] = coupled("displacements", i);
156  _var_objects[i] = getVar("displacements", i);
157  }
158  }
159  else
160  {
161  // Legacy parameter scheme for displacements
162  if (isParamValid("disp_x"))
163  {
164  _vars[0] = coupled("disp_x");
165  _var_objects[0] = getVar("disp_x", 0);
166  }
167  if (isParamValid("disp_y"))
168  {
169  _vars[1] = coupled("disp_y");
170  _var_objects[1] = getVar("disp_y", 0);
171  }
172  if (isParamValid("disp_z"))
173  {
174  _vars[2] = coupled("disp_z");
175  _var_objects[2] = getVar("disp_z", 0);
176  }
177 
178  mooseDeprecated("use the `displacements` parameter rather than the `disp_*` parameters (those "
179  "will go away with the deprecation of the Solid Mechanics module).");
180  }
181 
182  if (parameters.isParamValid("tangential_tolerance"))
183  _penetration_locator.setTangentialTolerance(getParam<Real>("tangential_tolerance"));
184 
185  if (parameters.isParamValid("normal_smoothing_distance"))
186  _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
187 
188  if (parameters.isParamValid("normal_smoothing_method"))
189  _penetration_locator.setNormalSmoothingMethod(
190  parameters.get<std::string>("normal_smoothing_method"));
191 
192  if (_formulation == CF_TANGENTIAL_PENALTY && _model != CM_COULOMB)
193  mooseError("The 'tangential_penalty' formulation can only be used with the 'coulomb' model");
194 
195  if (_model == CM_GLUED)
196  _penetration_locator.setUpdate(false);
197 
198  if (_friction_coefficient < 0)
199  mooseError("The friction coefficient must be nonnegative");
200 
201  // set _penalty_tangential to the value of _penalty for now
202  _penalty_tangential = _penalty;
203 
204  if (_formulation == CF_AUGMENTED_LAGRANGE)
205  {
206  if (_model == CM_GLUED)
207  mooseError("The Augmented Lagrangian contact formulation does not support GLUED case.");
208 
209  FEProblemBase * fe_problem = getParam<FEProblemBase *>("_fe_problem_base");
210  if (dynamic_cast<AugmentedLagrangianContactProblem *>(fe_problem) == NULL)
211  mooseError("The Augmented Lagrangian contact formulation must use "
212  "AugmentedLagrangianContactProblem.");
213 
214  if (!parameters.isParamValid("al_penetration_tolerance"))
215  mooseError("For Augmented Lagrangian contact, al_penetration_tolerance must be provided.");
216  else
217  _al_penetration_tolerance = parameters.get<Real>("al_penetration_tolerance");
218 
219  if (_model != CM_FRICTIONLESS)
220  {
221  if (!parameters.isParamValid("al_incremental_slip_tolerance") ||
222  !parameters.isParamValid("al_frictional_force_tolerance"))
223  {
224  mooseError("For the Augmented Lagrangian frictional contact formualton, "
225  "al_incremental_slip_tolerance and "
226  "al_frictional_force_tolerance must be provided.");
227  }
228  else
229  {
230  _al_incremental_slip_tolerance = parameters.get<Real>("al_incremental_slip_tolerance");
231  _al_frictional_force_tolerance = parameters.get<Real>("al_frictional_force_tolerance");
232  }
233  }
234  }
235 }
MechanicalContactConstraint::_residual_copy
NumericVector< Number > & _residual_copy
Definition: MechanicalContactConstraint.h:110
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::_contact_linesearch
std::shared_ptr< ContactLineSearchBase > _contact_linesearch
Definition: MechanicalContactConstraint.h:136
MechanicalContactConstraint::_var_objects
std::vector< MooseVariable * > _var_objects
Definition: MechanicalContactConstraint.h:116
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
MechanicalContactConstraint::_displaced_problem
MooseSharedPointer< DisplacedProblem > _displaced_problem
Definition: MechanicalContactConstraint.h:90
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
ContactMaster::contactFormulation
static ContactFormulation contactFormulation(std::string name)
Definition: ContactMaster.C:435
MechanicalContactConstraint::_vars
std::vector< unsigned int > _vars
Definition: MechanicalContactConstraint.h:115
MechanicalContactConstraint::_al_incremental_slip_tolerance
Real _al_incremental_slip_tolerance
The tolerance of the incremental slip for augmented Lagrangian method.
Definition: MechanicalContactConstraint.h:132
ContactMaster::contactModel
static ContactModel contactModel(std::string name)
Definition: ContactMaster.C:417
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::_aux_solution
const NumericVector< Number > * _aux_solution
Definition: MechanicalContactConstraint.h:120
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

Member Function Documentation

◆ addCouplingEntriesToJacobian()

virtual bool MechanicalContactConstraint::addCouplingEntriesToJacobian ( )
inlineoverridevirtual

Definition at line 84 of file MechanicalContactConstraint.h.

84 { return _master_slave_jacobian; }
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

◆ AugmentedLagrangianContactConverged()

bool MechanicalContactConstraint::AugmentedLagrangianContactConverged ( )
virtual

Definition at line 341 of file MechanicalContactConstraint.C.

342 {
343  Real contactResidual = 0.0;
344  unsigned int converged = 0;
345 
346  for (auto & pinfo_pair : _penetration_locator._penetration_info)
347  {
348  const dof_id_type slave_node_num = pinfo_pair.first;
349  PenetrationInfo * pinfo = pinfo_pair.second;
350 
351  // Skip this pinfo if there are no DOFs on this node.
352  if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
353  continue;
354 
355  const Real distance = pinfo->_normal * (pinfo->_closest_point - _mesh.nodeRef(slave_node_num));
356 
357  if (pinfo->isCaptured())
358  {
359  if (contactResidual < std::abs(distance))
360  contactResidual = std::abs(distance);
361 
362  // penetration < tol
363  if (contactResidual > _al_penetration_tolerance)
364  {
365  converged = 1;
366  break;
367  }
368 
369  if (_model == CM_COULOMB)
370  {
371  RealVectorValue contact_force_normal((pinfo->_contact_force * pinfo->_normal) *
372  pinfo->_normal);
373  RealVectorValue contact_force_tangential(pinfo->_contact_force - contact_force_normal);
374 
375  RealVectorValue tangential_inc_slip =
376  pinfo->_incremental_slip - (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
377 
378  const Real tan_mag(contact_force_tangential.norm());
379  const Real tangential_inc_slip_mag = tangential_inc_slip.norm();
380 
381  RealVectorValue distance_vec =
382  (pinfo->_normal * (_mesh.nodeRef(slave_node_num) - pinfo->_closest_point)) *
383  pinfo->_normal;
384 
385  Real penalty = getPenalty(*pinfo);
386  RealVectorValue pen_force_normal =
387  penalty * distance_vec + pinfo->_lagrange_multiplier * pinfo->_normal;
388 
389  // Frictional capacity
390  Real capacity(_friction_coefficient * (pen_force_normal * pinfo->_normal < 0
391  ? -pen_force_normal * pinfo->_normal
392  : 0.0));
393 
394  // incremental slip <= tol for all pinfo_pair such that tan_mag < capacity
395  if (MooseUtils::absoluteFuzzyLessThan(tan_mag, capacity) &&
396  pinfo->_mech_status == PenetrationInfo::MS_STICKING)
397  {
398  if (MooseUtils::absoluteFuzzyGreaterThan(tangential_inc_slip_mag,
399  _al_incremental_slip_tolerance))
400  {
401  converged = 2;
402  break;
403  }
404  }
405 
406  // for all pinfo_pair, tag_mag should be less than (1 + tol) * (capacity + tol)
407  if (tan_mag >
408  (1 + _al_frictional_force_tolerance) * (capacity + _al_frictional_force_tolerance))
409  {
410  converged = 3;
411  break;
412  }
413  }
414  }
415  }
416 
417  _communicator.max(converged);
418 
419  if (converged == 1)
420  _console
421  << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied \n";
422  else if (converged == 2)
423  _console
424  << "The Augmented Lagrangian contact tangential sliding enforcement is NOT satisfied \n";
425  else if (converged == 3)
426  _console << "The Augmented Lagrangian contact frictional force enforcement is NOT satisfied \n";
427  else
428  return true;
429 
430  return false;
431 }
MechanicalContactConstraint::_al_penetration_tolerance
Real _al_penetration_tolerance
The tolerance of the penetration for augmented Lagrangian method.
Definition: MechanicalContactConstraint.h:130
MechanicalContactConstraint::getPenalty
Real getPenalty(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1664
MechanicalContactConstraint::_vars
std::vector< unsigned int > _vars
Definition: MechanicalContactConstraint.h:115
MechanicalContactConstraint::_al_incremental_slip_tolerance
Real _al_incremental_slip_tolerance
The tolerance of the incremental slip for augmented Lagrangian method.
Definition: MechanicalContactConstraint.h:132
MechanicalContactConstraint::_al_frictional_force_tolerance
Real _al_frictional_force_tolerance
The tolerance of the frictional force for augmented Lagrangian method.
Definition: MechanicalContactConstraint.h:134

◆ computeContactForce()

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

Definition at line 507 of file MechanicalContactConstraint.C.

Referenced by shouldApply().

508 {
509  const Node * node = pinfo->_node;
510 
511  // Build up residual vector
512  RealVectorValue res_vec;
513  for (unsigned int i = 0; i < _mesh_dimension; ++i)
514  {
515  dof_id_type dof_number = node->dof_number(0, _vars[i], 0);
516  res_vec(i) = _residual_copy(dof_number) / _var_objects[i]->scalingFactor();
517  }
518 
519  RealVectorValue distance_vec(_mesh.nodeRef(node->id()) - pinfo->_closest_point);
520 
521  const Real gap_size = -1.0 * pinfo->_normal * distance_vec;
522 
523  // This is for preventing an increment of pinfo->_locked_this_step for nodes that are
524  // captured and released in this function
525  bool newly_captured = false;
526 
527  // Capture nodes that are newly in contact
528  if (update_contact_set && !pinfo->isCaptured() &&
529  MooseUtils::absoluteFuzzyGreaterEqual(gap_size, 0.0, _capture_tolerance))
530  {
531  newly_captured = true;
532  pinfo->capture();
533 
534  // Increment the lock count every time the node comes back into contact from not being in
535  // contact.
536  if (_formulation == CF_KINEMATIC || _formulation == CF_TANGENTIAL_PENALTY)
537  ++pinfo->_locked_this_step;
538  }
539 
540  if (!pinfo->isCaptured())
541  return;
542 
543  const Real penalty = getPenalty(*pinfo);
544  const Real penalty_slip = getTangentialPenalty(*pinfo);
545 
546  RealVectorValue pen_force(penalty * distance_vec);
547 
548  switch (_model)
549  {
550  case CM_FRICTIONLESS:
551  switch (_formulation)
552  {
553  case CF_KINEMATIC:
554  pinfo->_contact_force = -pinfo->_normal * (pinfo->_normal * res_vec);
555  break;
556 
557  case CF_PENALTY:
558  pinfo->_contact_force = pinfo->_normal * (pinfo->_normal * pen_force);
559  break;
560 
561  case CF_AUGMENTED_LAGRANGE:
562  pinfo->_contact_force =
563  (pinfo->_normal *
564  (pinfo->_normal * (pen_force + pinfo->_lagrange_multiplier * pinfo->_normal)));
565  break;
566 
567  default:
568  mooseError("Invalid contact formulation");
569  break;
570  }
571  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
572  break;
573  case CM_COULOMB:
574  switch (_formulation)
575  {
576  case CF_KINEMATIC:
577  {
578  // Frictional capacity
579  const Real capacity(_friction_coefficient *
580  (res_vec * pinfo->_normal > 0 ? res_vec * pinfo->_normal : 0));
581 
582  // Normal and tangential components of predictor force
583  pinfo->_contact_force = -res_vec;
584  RealVectorValue contact_force_normal((pinfo->_contact_force * pinfo->_normal) *
585  pinfo->_normal);
586  RealVectorValue contact_force_tangential(pinfo->_contact_force - contact_force_normal);
587 
588  RealVectorValue tangential_inc_slip =
589  pinfo->_incremental_slip -
590  (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
591 
592  // Magnitude of tangential predictor force
593  const Real tan_mag(contact_force_tangential.size());
594  const Real tangential_inc_slip_mag = tangential_inc_slip.size();
595  const Real slip_tol = capacity / penalty;
596  pinfo->_slip_tol = slip_tol;
597 
598  if ((tangential_inc_slip_mag > slip_tol || tan_mag > capacity) &&
599  (pinfo->_stick_locked_this_step < _stick_lock_iterations ||
600  tan_mag > capacity * _stick_unlock_factor))
601  {
602  if (pinfo->_stick_locked_this_step >= _stick_lock_iterations)
603  pinfo->_stick_locked_this_step = 0;
604 
605  // Check for scenario where node has slipped too far in a previous iteration
606  // for special treatment (only done if the incremental slip is non-trivial)
607  bool slipped_too_far = false;
608  RealVectorValue slip_inc_direction;
609  if (tangential_inc_slip_mag > slip_tol)
610  {
611  slip_inc_direction = tangential_inc_slip / tangential_inc_slip_mag;
612  Real slip_dot_tang_force = slip_inc_direction * contact_force_tangential;
613  if (slip_dot_tang_force < capacity)
614  slipped_too_far = true;
615  }
616 
617  if (slipped_too_far) // slip back along slip increment
618  pinfo->_contact_force = contact_force_normal + capacity * slip_inc_direction;
619  else
620  {
621  if (tan_mag > 0) // slip along tangential force direction
622  pinfo->_contact_force =
623  contact_force_normal + capacity * contact_force_tangential / tan_mag;
624  else // treat as frictionless
625  pinfo->_contact_force = contact_force_normal;
626  }
627  if (capacity == 0)
628  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
629  else
630  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
631  }
632  else
633  {
634  if (pinfo->_mech_status != PenetrationInfo::MS_STICKING &&
635  pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
636  ++pinfo->_stick_locked_this_step;
637  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
638  }
639  break;
640  }
641 
642  case CF_PENALTY:
643  {
644  distance_vec = pinfo->_incremental_slip +
645  (pinfo->_normal * (_mesh.nodeRef(node->id()) - pinfo->_closest_point)) *
646  pinfo->_normal;
647  pen_force = penalty * distance_vec;
648 
649  // Frictional capacity
650  // const Real capacity( _friction_coefficient * (pen_force * pinfo->_normal < 0 ?
651  // -pen_force * pinfo->_normal : 0) );
652  const Real capacity(_friction_coefficient *
653  (res_vec * pinfo->_normal > 0 ? res_vec * pinfo->_normal : 0));
654 
655  // Elastic predictor
656  pinfo->_contact_force =
657  pen_force + (pinfo->_contact_force_old -
658  pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old));
659  RealVectorValue contact_force_normal((pinfo->_contact_force * pinfo->_normal) *
660  pinfo->_normal);
661  RealVectorValue contact_force_tangential(pinfo->_contact_force - contact_force_normal);
662 
663  // Tangential magnitude of elastic predictor
664  const Real tan_mag(contact_force_tangential.norm());
665 
666  if (tan_mag > capacity)
667  {
668  pinfo->_contact_force =
669  contact_force_normal + capacity * contact_force_tangential / tan_mag;
670  if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
671  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
672  else
673  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
674  }
675  else
676  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
677  break;
678  }
679 
680  case CF_AUGMENTED_LAGRANGE:
681  {
682  distance_vec = (pinfo->_normal * (_mesh.nodeRef(node->id()) - pinfo->_closest_point)) *
683  pinfo->_normal;
684 
685  RealVectorValue contact_force_normal =
686  penalty * distance_vec + pinfo->_lagrange_multiplier * pinfo->_normal;
687 
688  RealVectorValue tangential_inc_slip =
689  pinfo->_incremental_slip -
690  (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
691 
692  RealVectorValue contact_force_tangential =
693  pinfo->_lagrange_multiplier_slip +
694  (pinfo->_contact_force_old -
695  pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old));
696 
697  RealVectorValue inc_pen_force_tangential = penalty_slip * tangential_inc_slip;
698 
699  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
700  pinfo->_contact_force =
701  contact_force_normal + contact_force_tangential + inc_pen_force_tangential;
702  else
703  pinfo->_contact_force = contact_force_normal + contact_force_tangential;
704 
705  break;
706  }
707 
708  case CF_TANGENTIAL_PENALTY:
709  {
710  // Frictional capacity (kinematic formulation)
711  const Real capacity = _friction_coefficient * std::max(res_vec * pinfo->_normal, 0.0);
712 
713  // Normal component of contact force (kinematic formulation)
714  RealVectorValue contact_force_normal((-res_vec * pinfo->_normal) * pinfo->_normal);
715 
716  // Predictor tangential component of contact force (penalty formulation)
717  RealVectorValue inc_pen_force_tangential = penalty * pinfo->_incremental_slip;
718  RealVectorValue contact_force_tangential =
719  inc_pen_force_tangential +
720  (pinfo->_contact_force_old -
721  pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old));
722 
723  // Magnitude of tangential predictor
724  const Real tan_mag(contact_force_tangential.norm());
725 
726  if (tan_mag > capacity)
727  {
728  pinfo->_contact_force =
729  contact_force_normal + capacity * contact_force_tangential / tan_mag;
730  if (MooseUtils::absoluteFuzzyEqual(capacity, 0))
731  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
732  else
733  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
734  }
735  else
736  {
737  pinfo->_contact_force = contact_force_normal + contact_force_tangential;
738  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
739  }
740  break;
741  }
742 
743  default:
744  mooseError("Invalid contact formulation");
745  break;
746  }
747  break;
748 
749  case CM_GLUED:
750  switch (_formulation)
751  {
752  case CF_KINEMATIC:
753  pinfo->_contact_force = -res_vec;
754  break;
755 
756  case CF_PENALTY:
757  pinfo->_contact_force = pen_force;
758  break;
759 
760  case CF_AUGMENTED_LAGRANGE:
761  pinfo->_contact_force =
762  pen_force + pinfo->_lagrange_multiplier * distance_vec / distance_vec.norm();
763  break;
764 
765  default:
766  mooseError("Invalid contact formulation");
767  break;
768  }
769  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
770  break;
771 
772  default:
773  mooseError("Invalid or unavailable contact model");
774  break;
775  }
776 
777  // Release
778  if (update_contact_set && _model != CM_GLUED && pinfo->isCaptured() && !newly_captured &&
779  _tension_release >= 0.0 && (_contact_linesearch ? true : pinfo->_locked_this_step < 2))
780  {
781  const Real contact_pressure = -(pinfo->_normal * pinfo->_contact_force) / nodalArea(*pinfo);
782  if (-contact_pressure >= _tension_release)
783  {
784  pinfo->release();
785  pinfo->_contact_force.zero();
786  }
787  }
788 }
MechanicalContactConstraint::_residual_copy
NumericVector< Number > & _residual_copy
Definition: MechanicalContactConstraint.h:110
MechanicalContactConstraint::nodalArea
Real nodalArea(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1646
MechanicalContactConstraint::_formulation
const ContactFormulation _formulation
Definition: MechanicalContactConstraint.h:98
MechanicalContactConstraint::getPenalty
Real getPenalty(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1664
MechanicalContactConstraint::_contact_linesearch
std::shared_ptr< ContactLineSearchBase > _contact_linesearch
Definition: MechanicalContactConstraint.h:136
MechanicalContactConstraint::_var_objects
std::vector< MooseVariable * > _var_objects
Definition: MechanicalContactConstraint.h:116
MechanicalContactConstraint::getTangentialPenalty
Real getTangentialPenalty(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1674
MechanicalContactConstraint::_vars
std::vector< unsigned int > _vars
Definition: MechanicalContactConstraint.h:115

◆ 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 }
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::computeQpJacobian
virtual Real computeQpJacobian(Moose::ConstraintJacobianType type) override
Definition: MechanicalContactConstraint.C:843
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

◆ 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 }
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::computeQpOffDiagJacobian
virtual Real computeQpOffDiagJacobian(Moose::ConstraintJacobianType type, unsigned int jvar) override
Compute off-diagonal Jacobian entries.
Definition: MechanicalContactConstraint.C:1291
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

◆ computeQpJacobian()

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

Definition at line 843 of file MechanicalContactConstraint.C.

Referenced by computeJacobian().

844 {
845  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
846 
847  const Real penalty = getPenalty(*pinfo);
848  const Real penalty_slip = getTangentialPenalty(*pinfo);
849 
850  switch (type)
851  {
852  case Moose::SlaveSlave:
853  switch (_model)
854  {
855  case CM_FRICTIONLESS:
856  switch (_formulation)
857  {
858  case CF_KINEMATIC:
859  {
860  RealVectorValue jac_vec;
861  for (unsigned int i = 0; i < _mesh_dimension; ++i)
862  {
863  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
864  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
865  }
866  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
867  (_phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp]) *
868  pinfo->_normal(_component) * pinfo->_normal(_component);
869  }
870 
871  case CF_PENALTY:
872  case CF_AUGMENTED_LAGRANGE:
873  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
874  pinfo->_normal(_component) * pinfo->_normal(_component);
875 
876  default:
877  mooseError("Invalid contact formulation");
878  }
879 
880  case CM_COULOMB:
881  switch (_formulation)
882  {
883  case CF_KINEMATIC:
884  {
885  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
886  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
887  {
888  RealVectorValue jac_vec;
889  for (unsigned int i = 0; i < _mesh_dimension; ++i)
890  {
891  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
892  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
893  }
894  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
895  (_phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp]) *
896  pinfo->_normal(_component) * pinfo->_normal(_component);
897  }
898  else
899  {
900  const Real curr_jac = (*_jacobian)(
901  _current_node->dof_number(0, _vars[_component], 0), _connected_dof_indices[_j]);
902  return -curr_jac + _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp];
903  }
904  }
905 
906  case CF_PENALTY:
907  {
908  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
909  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
910  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
911  pinfo->_normal(_component) * pinfo->_normal(_component);
912  else
913  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp];
914  }
915  case CF_AUGMENTED_LAGRANGE:
916  {
917  Real normal_comp = _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
918  pinfo->_normal(_component) * pinfo->_normal(_component);
919 
920  Real tang_comp = 0.0;
921  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
922  tang_comp = _phi_slave[_j][_qp] * penalty_slip * _test_slave[_i][_qp] *
923  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
924  return normal_comp + tang_comp;
925  }
926 
927  case CF_TANGENTIAL_PENALTY:
928  {
929  RealVectorValue jac_vec;
930  for (unsigned int i = 0; i < _mesh_dimension; ++i)
931  {
932  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
933  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
934  }
935  Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
936  (_phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp]) *
937  pinfo->_normal(_component) * pinfo->_normal(_component);
938 
939  Real tang_comp = 0.0;
940  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
941  tang_comp = _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
942  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
943 
944  return normal_comp + tang_comp;
945  }
946 
947  default:
948  mooseError("Invalid contact formulation");
949  }
950 
951  case CM_GLUED:
952  switch (_formulation)
953  {
954  case CF_KINEMATIC:
955  {
956  const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
957  _connected_dof_indices[_j]);
958  return -curr_jac + _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp];
959  }
960 
961  case CF_PENALTY:
962  case CF_AUGMENTED_LAGRANGE:
963  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp];
964 
965  default:
966  mooseError("Invalid contact formulation");
967  }
968  default:
969  mooseError("Invalid or unavailable contact model");
970  }
971 
972  case Moose::SlaveMaster:
973  switch (_model)
974  {
975  case CM_FRICTIONLESS:
976  switch (_formulation)
977  {
978  case CF_KINEMATIC:
979  {
980  Node * curr_master_node = _current_master->get_node(_j);
981 
982  RealVectorValue jac_vec;
983  for (unsigned int i = 0; i < _mesh_dimension; ++i)
984  {
985  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
986  jac_vec(i) =
987  (*_jacobian)(dof_number, curr_master_node->dof_number(0, _vars[_component], 0));
988  }
989  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
990  (_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp]) *
991  pinfo->_normal(_component) * pinfo->_normal(_component);
992  }
993 
994  case CF_PENALTY:
995  case CF_AUGMENTED_LAGRANGE:
996  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
997  pinfo->_normal(_component) * pinfo->_normal(_component);
998 
999  default:
1000  mooseError("Invalid contact formulation");
1001  }
1002 
1003  case CM_COULOMB:
1004  switch (_formulation)
1005  {
1006  case CF_KINEMATIC:
1007  {
1008  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1009  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1010  {
1011  Node * curr_master_node = _current_master->get_node(_j);
1012 
1013  RealVectorValue jac_vec;
1014  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1015  {
1016  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1017  jac_vec(i) = (*_jacobian)(dof_number,
1018  curr_master_node->dof_number(0, _vars[_component], 0));
1019  }
1020  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1021  (_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp]) *
1022  pinfo->_normal(_component) * pinfo->_normal(_component);
1023  }
1024  else
1025  {
1026  Node * curr_master_node = _current_master->get_node(_j);
1027  const Real curr_jac =
1028  (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1029  curr_master_node->dof_number(0, _vars[_component], 0));
1030  return -curr_jac - _phi_master[_j][_qp] * penalty * _test_slave[_i][_qp];
1031  }
1032  }
1033 
1034  case CF_PENALTY:
1035  {
1036  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1037  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1038  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
1039  pinfo->_normal(_component) * pinfo->_normal(_component);
1040  else
1041  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp];
1042  }
1043  case CF_AUGMENTED_LAGRANGE:
1044  {
1045  Real normal_comp = -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
1046  pinfo->_normal(_component) * pinfo->_normal(_component);
1047 
1048  Real tang_comp = 0.0;
1049  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1050  tang_comp = -_phi_master[_j][_qp] * penalty_slip * _test_slave[_i][_qp] *
1051  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1052  return normal_comp + tang_comp;
1053  }
1054 
1055  case CF_TANGENTIAL_PENALTY:
1056  {
1057  Node * curr_master_node = _current_master->get_node(_j);
1058 
1059  RealVectorValue jac_vec;
1060  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1061  {
1062  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1063  jac_vec(i) =
1064  (*_jacobian)(dof_number, curr_master_node->dof_number(0, _vars[_component], 0));
1065  }
1066  Real normal_comp = -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1067  (_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp]) *
1068  pinfo->_normal(_component) * pinfo->_normal(_component);
1069 
1070  Real tang_comp = 0.0;
1071  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1072  tang_comp = -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
1073  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1074 
1075  return normal_comp + tang_comp;
1076  }
1077 
1078  default:
1079  mooseError("Invalid contact formulation");
1080  }
1081  case CM_GLUED:
1082  switch (_formulation)
1083  {
1084  case CF_KINEMATIC:
1085  {
1086  Node * curr_master_node = _current_master->get_node(_j);
1087  const Real curr_jac =
1088  (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1089  curr_master_node->dof_number(0, _vars[_component], 0));
1090  return -curr_jac - _phi_master[_j][_qp] * penalty * _test_slave[_i][_qp];
1091  }
1092 
1093  case CF_PENALTY:
1094  case CF_AUGMENTED_LAGRANGE:
1095  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp];
1096 
1097  default:
1098  mooseError("Invalid contact formulation");
1099  }
1100 
1101  default:
1102  mooseError("Invalid or unavailable contact model");
1103  }
1104 
1105  case Moose::MasterSlave:
1106  switch (_model)
1107  {
1108  case CM_FRICTIONLESS:
1109  switch (_formulation)
1110  {
1111  case CF_KINEMATIC:
1112  {
1113  RealVectorValue jac_vec;
1114  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1115  {
1116  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1117  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1118  }
1119  return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1120  _test_master[_i][_qp];
1121  }
1122 
1123  case CF_PENALTY:
1124  case CF_AUGMENTED_LAGRANGE:
1125  return -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp] *
1126  pinfo->_normal(_component) * pinfo->_normal(_component);
1127 
1128  default:
1129  mooseError("Invalid contact formulation");
1130  }
1131  case CM_COULOMB:
1132  switch (_formulation)
1133  {
1134  case CF_KINEMATIC:
1135  {
1136  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1137  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1138  {
1139  RealVectorValue jac_vec;
1140  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1141  {
1142  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1143  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1144  }
1145  return pinfo->_normal(_component) * (pinfo->_normal * jac_vec) *
1146  _test_master[_i][_qp];
1147  }
1148  else
1149  {
1150  const Real slave_jac = (*_jacobian)(
1151  _current_node->dof_number(0, _vars[_component], 0), _connected_dof_indices[_j]);
1152  return slave_jac * _test_master[_i][_qp];
1153  }
1154  }
1155 
1156  case CF_PENALTY:
1157  {
1158  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1159  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1160  return -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp] *
1161  pinfo->_normal(_component) * pinfo->_normal(_component);
1162  else
1163  return -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp];
1164  }
1165  case CF_AUGMENTED_LAGRANGE:
1166  {
1167  Real normal_comp = -_phi_slave[_j][_qp] * penalty * _test_master[_i][_qp] *
1168  pinfo->_normal(_component) * pinfo->_normal(_component);
1169 
1170  Real tang_comp = 0.0;
1171  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1172  tang_comp = -_phi_slave[_j][_qp] * penalty_slip * _test_master[_i][_qp] *
1173  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1174  return normal_comp + tang_comp;
1175  }
1176 
1177  case CF_TANGENTIAL_PENALTY:
1178  {
1179  RealVectorValue jac_vec;
1180  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1181  {
1182  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1183  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1184  }
1185  Real normal_comp =
1186  pinfo->_normal(_component) * (pinfo->_normal * jac_vec) * _test_master[_i][_qp];
1187 
1188  Real tang_comp = 0.0;
1189  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1190  tang_comp = -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp] *
1191  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1192 
1193  return normal_comp + tang_comp;
1194  }
1195 
1196  default:
1197  mooseError("Invalid contact formulation");
1198  }
1199 
1200  case CM_GLUED:
1201  switch (_formulation)
1202  {
1203  case CF_KINEMATIC:
1204  {
1205  const Real slave_jac = (*_jacobian)(
1206  _current_node->dof_number(0, _vars[_component], 0), _connected_dof_indices[_j]);
1207  return slave_jac * _test_master[_i][_qp];
1208  }
1209 
1210  case CF_PENALTY:
1211  case CF_AUGMENTED_LAGRANGE:
1212  return -_test_master[_i][_qp] * penalty * _phi_slave[_j][_qp];
1213 
1214  default:
1215  mooseError("Invalid contact formulation");
1216  }
1217 
1218  default:
1219  mooseError("Invalid or unavailable contact model");
1220  }
1221 
1222  case Moose::MasterMaster:
1223  switch (_model)
1224  {
1225  case CM_FRICTIONLESS:
1226  switch (_formulation)
1227  {
1228  case CF_KINEMATIC:
1229  return 0.0;
1230 
1231  case CF_PENALTY:
1232  case CF_AUGMENTED_LAGRANGE:
1233  return _test_master[_i][_qp] * penalty * _phi_master[_j][_qp] *
1234  pinfo->_normal(_component) * pinfo->_normal(_component);
1235 
1236  default:
1237  mooseError("Invalid contact formulation");
1238  }
1239 
1240  case CM_COULOMB:
1241  case CM_GLUED:
1242  switch (_formulation)
1243  {
1244  case CF_KINEMATIC:
1245  return 0.0;
1246 
1247  case CF_PENALTY:
1248  {
1249  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1250  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
1251  return _test_master[_i][_qp] * penalty * _phi_master[_j][_qp] *
1252  pinfo->_normal(_component) * pinfo->_normal(_component);
1253  else
1254  return _test_master[_i][_qp] * penalty * _phi_master[_j][_qp];
1255  }
1256 
1257  case CF_TANGENTIAL_PENALTY:
1258  {
1259  Real tang_comp = 0.0;
1260  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1261  tang_comp = _test_master[_i][_qp] * penalty * _phi_master[_j][_qp] *
1262  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1263  return tang_comp; // normal component is zero
1264  }
1265 
1266  case CF_AUGMENTED_LAGRANGE:
1267  {
1268  Real normal_comp = _phi_master[_j][_qp] * penalty * _test_master[_i][_qp] *
1269  pinfo->_normal(_component) * pinfo->_normal(_component);
1270 
1271  Real tang_comp = 0.0;
1272  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1273  tang_comp = _phi_master[_j][_qp] * penalty_slip * _test_master[_i][_qp] *
1274  (1.0 - pinfo->_normal(_component) * pinfo->_normal(_component));
1275  return normal_comp + tang_comp;
1276  }
1277 
1278  default:
1279  mooseError("Invalid contact formulation");
1280  }
1281 
1282  default:
1283  mooseError("Invalid or unavailable contact model");
1284  }
1285  }
1286 
1287  return 0.0;
1288 }
MechanicalContactConstraint::_formulation
const ContactFormulation _formulation
Definition: MechanicalContactConstraint.h:98
MechanicalContactConstraint::getPenalty
Real getPenalty(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1664
MechanicalContactConstraint::getTangentialPenalty
Real getTangentialPenalty(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1674
MechanicalContactConstraint::_vars
std::vector< unsigned int > _vars
Definition: MechanicalContactConstraint.h:115

◆ 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 1291 of file MechanicalContactConstraint.C.

Referenced by computeOffDiagJacobian().

1293 {
1294  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
1295 
1296  const Real penalty = getPenalty(*pinfo);
1297  const Real penalty_slip = getTangentialPenalty(*pinfo);
1298 
1299  unsigned int coupled_component;
1300  Real normal_component_in_coupled_var_dir = 1.0;
1301  if (getCoupledVarComponent(jvar, coupled_component))
1302  normal_component_in_coupled_var_dir = pinfo->_normal(coupled_component);
1303 
1304  switch (type)
1305  {
1306  case Moose::SlaveSlave:
1307  switch (_model)
1308  {
1309  case CM_FRICTIONLESS:
1310  switch (_formulation)
1311  {
1312  case CF_KINEMATIC:
1313  {
1314  RealVectorValue jac_vec;
1315  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1316  {
1317  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1318  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1319  }
1320  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1321  (_phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp]) *
1322  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1323  }
1324 
1325  case CF_PENALTY:
1326  case CF_AUGMENTED_LAGRANGE:
1327  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
1328  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1329 
1330  default:
1331  mooseError("Invalid contact formulation");
1332  }
1333 
1334  case CM_COULOMB:
1335  {
1336  if ((_formulation == CF_KINEMATIC || _formulation == CF_TANGENTIAL_PENALTY) &&
1337  (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1338  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1339  {
1340  RealVectorValue jac_vec;
1341  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1342  {
1343  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1344  jac_vec(i) = (*_jacobian)(dof_number, _connected_dof_indices[_j]);
1345  }
1346  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) +
1347  (_phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp]) *
1348  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1349  }
1350  else if ((_formulation == CF_PENALTY) &&
1351  (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1352  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1353  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
1354  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1355 
1356  else if (_formulation == CF_AUGMENTED_LAGRANGE)
1357  {
1358  Real normal_comp = _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp] *
1359  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1360 
1361  Real tang_comp = 0.0;
1362  if (pinfo->_mech_status == PenetrationInfo::MS_STICKING)
1363  tang_comp = _phi_slave[_j][_qp] * penalty_slip * _test_slave[_i][_qp] *
1364  (-pinfo->_normal(_component) * normal_component_in_coupled_var_dir);
1365  return normal_comp + tang_comp;
1366  }
1367  else
1368  {
1369  const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1370  _connected_dof_indices[_j]);
1371  return -curr_jac;
1372  }
1373  }
1374 
1375  case CM_GLUED:
1376  {
1377  const Real curr_jac = (*_jacobian)(_current_node->dof_number(0, _vars[_component], 0),
1378  _connected_dof_indices[_j]);
1379  return -curr_jac;
1380  }
1381  default:
1382  mooseError("Invalid or unavailable contact model");
1383  }
1384 
1385  case Moose::SlaveMaster:
1386  switch (_model)
1387  {
1388  case CM_FRICTIONLESS:
1389  switch (_formulation)
1390  {
1391  case CF_KINEMATIC:
1392  {
1393  Node * curr_master_node = _current_master->get_node(_j);
1394 
1395  RealVectorValue jac_vec;
1396  for (unsigned int i = 0; i < _mesh_dimension; ++i)
1397  {
1398  dof_id_type dof_number = _current_node->dof_number(0, _vars[i], 0);
1399  jac_vec(i) =
1400  (*_jacobian)(dof_number, curr_master_node->dof_number(0, _vars[_component], 0));
1401  }
1402  return -pinfo->_normal(_component) * (pinfo->_normal * jac_vec) -
1403  (_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp]) *
1404  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1405  }
1406 
1407  case CF_PENALTY:
1408  case CF_AUGMENTED_LAGRANGE:
1409  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp] *
1410  pinfo->_normal(_component) * normal_component_in_coupled_var_dir;
1411 
1412  default:
1413  mooseError("Invalid contact formulation");
1414  }
1415 
1416  case CM_COULOMB:
1417  if ((_formulation == CF_KINEMATIC || _formulation == CF_TANGENTIAL_PENALTY) &&
1418  (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
1419  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION))
1420  {
1421  Node * curr_master_node = _current_master->get_node(_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 == CF_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 == CF_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 CM_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 CM_FRICTIONLESS:
1465  switch (_formulation)
1466  {
1467  case CF_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 CF_PENALTY:
1480  case CF_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 CM_COULOMB:
1488  switch (_formulation)
1489  {
1490  case CF_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 CF_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 CF_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 CF_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 CM_GLUED:
1555  switch (_formulation)
1556  {
1557  case CF_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 CF_PENALTY:
1565  case CF_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 CM_FRICTIONLESS:
1580  switch (_formulation)
1581  {
1582  case CF_KINEMATIC:
1583  return 0.0;
1584 
1585  case CF_PENALTY:
1586  case CF_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 CM_COULOMB:
1595  {
1596  if (_formulation == CF_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 == CF_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 CM_GLUED:
1618  if (_formulation == CF_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 == CF_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 }
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::_formulation
const ContactFormulation _formulation
Definition: MechanicalContactConstraint.h:98
MechanicalContactConstraint::getPenalty
Real getPenalty(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1664
MechanicalContactConstraint::getTangentialPenalty
Real getTangentialPenalty(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1674
MechanicalContactConstraint::_vars
std::vector< unsigned int > _vars
Definition: MechanicalContactConstraint.h:115

◆ computeQpResidual()

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

Definition at line 797 of file MechanicalContactConstraint.C.

798 {
799  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
800  Real resid = pinfo->_contact_force(_component);
801  switch (type)
802  {
803  case Moose::Slave:
804  if (_formulation == CF_KINEMATIC)
805  {
806  RealVectorValue distance_vec(*_current_node - pinfo->_closest_point);
807  const Real penalty = getPenalty(*pinfo);
808  RealVectorValue pen_force(penalty * distance_vec);
809 
810  if (_model == CM_FRICTIONLESS)
811  resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
812 
813  else if (_model == CM_COULOMB)
814  {
815  distance_vec = distance_vec - pinfo->_incremental_slip;
816  pen_force = penalty * distance_vec;
817  if (pinfo->_mech_status == PenetrationInfo::MS_SLIPPING ||
818  pinfo->_mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
819  resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
820  else
821  resid += pen_force(_component);
822  }
823  else if (_model == CM_GLUED)
824  resid += pen_force(_component);
825  }
826  else if (_formulation == CF_TANGENTIAL_PENALTY && _model == CM_COULOMB)
827  {
828  RealVectorValue distance_vec(*_current_node - pinfo->_closest_point);
829  const Real penalty = getPenalty(*pinfo);
830  RealVectorValue pen_force(penalty * distance_vec);
831  resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
832  }
833  return _test_slave[_i][_qp] * resid;
834 
835  case Moose::Master:
836  return _test_master[_i][_qp] * -resid;
837  }
838 
839  return 0.0;
840 }
MechanicalContactConstraint::_formulation
const ContactFormulation _formulation
Definition: MechanicalContactConstraint.h:98
MechanicalContactConstraint::getPenalty
Real getPenalty(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1664

◆ computeQpSlaveValue()

Real MechanicalContactConstraint::computeQpSlaveValue ( )
overridevirtual

Definition at line 791 of file MechanicalContactConstraint.C.

792 {
793  return _u_slave[_qp];
794 }

◆ 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.

Referenced by computeJacobian(), and computeOffDiagJacobian().

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 }
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
MaterialTensorCalculatorTools::component
Real component(const SymmTensor &symm_tensor, unsigned int index)
Definition: MaterialTensorCalculatorTools.C:16
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
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::_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

◆ 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.

Referenced by computeQpOffDiagJacobian(), and getConnectedDofIndices().

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 }
MaterialTensorCalculatorTools::component
Real component(const SymmTensor &symm_tensor, unsigned int index)
Definition: MaterialTensorCalculatorTools.C:16
MechanicalContactConstraint::_vars
std::vector< unsigned int > _vars
Definition: MechanicalContactConstraint.h:115

◆ getPenalty()

Real MechanicalContactConstraint::getPenalty ( PenetrationInfo &  pinfo)
protected

Definition at line 1664 of file MechanicalContactConstraint.C.

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

1665 {
1666  Real penalty = _penalty;
1667  if (_normalize_penalty)
1668  penalty *= nodalArea(pinfo);
1669 
1670  return penalty;
1671 }
MechanicalContactConstraint::nodalArea
Real nodalArea(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1646

◆ getTangentialPenalty()

Real MechanicalContactConstraint::getTangentialPenalty ( PenetrationInfo &  pinfo)
protected

Definition at line 1674 of file MechanicalContactConstraint.C.

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

1675 {
1676  Real penalty = _penalty_tangential;
1677  if (_normalize_penalty)
1678  penalty *= nodalArea(pinfo);
1679 
1680  return penalty;
1681 }
MechanicalContactConstraint::nodalArea
Real nodalArea(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1646

◆ jacobianSetup()

void MechanicalContactConstraint::jacobianSetup ( )
overridevirtual

Definition at line 254 of file MechanicalContactConstraint.C.

255 {
256  if (_component == 0)
257  {
258  if (_update_stateful_data)
259  updateContactStatefulData();
260  _update_stateful_data = true;
261  }
262 }
MechanicalContactConstraint::updateContactStatefulData
virtual void updateContactStatefulData(bool beginning_of_step=false)
Definition: MechanicalContactConstraint.C:434

◆ nodalArea()

Real MechanicalContactConstraint::nodalArea ( PenetrationInfo &  pinfo)
protected

Definition at line 1646 of file MechanicalContactConstraint.C.

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

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 }

◆ 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 }
MechanicalContactConstraint::_current_contact_state
std::set< dof_id_type > _current_contact_state
Definition: MechanicalContactConstraint.h:137
MechanicalContactConstraint::_old_contact_state
std::set< dof_id_type > _old_contact_state
Definition: MechanicalContactConstraint.h:138
MechanicalContactConstraint::_contact_linesearch
std::shared_ptr< ContactLineSearchBase > _contact_linesearch
Definition: MechanicalContactConstraint.h:136

◆ shouldApply()

bool MechanicalContactConstraint::shouldApply ( )
override

Definition at line 473 of file MechanicalContactConstraint.C.

474 {
475  bool in_contact = false;
476 
477  std::map<dof_id_type, PenetrationInfo *>::iterator found =
478  _penetration_locator._penetration_info.find(_current_node->id());
479  if (found != _penetration_locator._penetration_info.end())
480  {
481  PenetrationInfo * pinfo = found->second;
482  if (pinfo != NULL)
483  {
484  bool is_nonlinear = _fe_problem.computingNonlinearResid();
485 
486  // This computes the contact force once per constraint, rather than once per quad point
487  // and for both master and slave cases.
488  if (_component == 0)
489  computeContactForce(pinfo, is_nonlinear);
490 
491  if (pinfo->isCaptured())
492  {
493  in_contact = true;
494  if (is_nonlinear)
495  {
496  Threads::spin_mutex::scoped_lock lock(_contact_set_mutex);
497  _current_contact_state.insert(pinfo->_node->id());
498  }
499  }
500  }
501  }
502 
503  return in_contact;
504 }
MechanicalContactConstraint::_current_contact_state
std::set< dof_id_type > _current_contact_state
Definition: MechanicalContactConstraint.h:137
MechanicalContactConstraint::_contact_set_mutex
static Threads::spin_mutex _contact_set_mutex
Definition: MechanicalContactConstraint.h:141
MechanicalContactConstraint::computeContactForce
void computeContactForce(PenetrationInfo *pinfo, bool update_contact_set)
Definition: MechanicalContactConstraint.C:507

◆ timestepSetup()

void MechanicalContactConstraint::timestepSetup ( )
overridevirtual

Definition at line 238 of file MechanicalContactConstraint.C.

239 {
240  if (_component == 0)
241  {
242  updateContactStatefulData(true);
243  if (_formulation == CF_AUGMENTED_LAGRANGE)
244  updateAugmentedLagrangianMultiplier(true);
245 
246  _update_stateful_data = false;
247 
248  if (_contact_linesearch)
249  _contact_linesearch->reset();
250  }
251 }
MechanicalContactConstraint::_formulation
const ContactFormulation _formulation
Definition: MechanicalContactConstraint.h:98
MechanicalContactConstraint::updateAugmentedLagrangianMultiplier
virtual void updateAugmentedLagrangianMultiplier(bool beginning_of_step=false)
Definition: MechanicalContactConstraint.C:265
MechanicalContactConstraint::_contact_linesearch
std::shared_ptr< ContactLineSearchBase > _contact_linesearch
Definition: MechanicalContactConstraint.h:136
MechanicalContactConstraint::updateContactStatefulData
virtual void updateContactStatefulData(bool beginning_of_step=false)
Definition: MechanicalContactConstraint.C:434

◆ updateAugmentedLagrangianMultiplier()

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

Definition at line 265 of file MechanicalContactConstraint.C.

Referenced by timestepSetup().

266 {
267  for (auto & pinfo_pair : _penetration_locator._penetration_info)
268  {
269  const dof_id_type slave_node_num = pinfo_pair.first;
270  PenetrationInfo * pinfo = pinfo_pair.second;
271 
272  if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
273  continue;
274 
275  const Real distance = pinfo->_normal * (pinfo->_closest_point - _mesh.nodeRef(slave_node_num));
276 
277  if (beginning_of_step && _model == CM_COULOMB)
278  {
279  pinfo->_lagrange_multiplier_slip.zero();
280  if (pinfo->isCaptured())
281  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
282  }
283 
284  if (pinfo->isCaptured())
285  {
286  if (_model == CM_FRICTIONLESS)
287  pinfo->_lagrange_multiplier -= getPenalty(*pinfo) * distance;
288 
289  if (_model == CM_COULOMB)
290  {
291  if (!beginning_of_step)
292  {
293  Real penalty = getPenalty(*pinfo);
294  RealVectorValue pen_force_normal =
295  penalty * (-distance) * pinfo->_normal + pinfo->_lagrange_multiplier * pinfo->_normal;
296 
297  // update normal lagrangian multiplier
298  pinfo->_lagrange_multiplier += penalty * (-distance);
299 
300  // Frictional capacity
301  const Real capacity(_friction_coefficient * (pen_force_normal * pinfo->_normal < 0
302  ? -pen_force_normal * pinfo->_normal
303  : 0));
304 
305  RealVectorValue tangential_inc_slip =
306  pinfo->_incremental_slip -
307  (pinfo->_incremental_slip * pinfo->_normal) * pinfo->_normal;
308 
309  Real penalty_slip = getTangentialPenalty(*pinfo);
310 
311  RealVectorValue inc_pen_force_tangential =
312  pinfo->_lagrange_multiplier_slip + penalty_slip * tangential_inc_slip;
313 
314  RealVectorValue tau_old = pinfo->_contact_force_old -
315  pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old);
316 
317  RealVectorValue contact_force_tangential = inc_pen_force_tangential + tau_old;
318  const Real tan_mag(contact_force_tangential.norm());
319 
320  if (tan_mag > capacity * (_al_frictional_force_tolerance + 1.0))
321  {
322  pinfo->_lagrange_multiplier_slip =
323  -tau_old + capacity * contact_force_tangential / tan_mag;
324  if (MooseUtils::absoluteFuzzyEqual(capacity, 0.0))
325  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
326  else
327  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING_FRICTION;
328  }
329  else
330  {
331  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
332  pinfo->_lagrange_multiplier_slip += penalty_slip * tangential_inc_slip;
333  }
334  }
335  }
336  }
337  }
338 }
MechanicalContactConstraint::getPenalty
Real getPenalty(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1664
MechanicalContactConstraint::getTangentialPenalty
Real getTangentialPenalty(PenetrationInfo &pinfo)
Definition: MechanicalContactConstraint.C:1674
MechanicalContactConstraint::_vars
std::vector< unsigned int > _vars
Definition: MechanicalContactConstraint.h:115
MechanicalContactConstraint::_al_frictional_force_tolerance
Real _al_frictional_force_tolerance
The tolerance of the frictional force for augmented Lagrangian method.
Definition: MechanicalContactConstraint.h:134

◆ updateContactStatefulData()

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

Definition at line 434 of file MechanicalContactConstraint.C.

Referenced by jacobianSetup(), and timestepSetup().

435 {
436  for (auto & pinfo_pair : _penetration_locator._penetration_info)
437  {
438  PenetrationInfo * pinfo = pinfo_pair.second;
439 
440  // Skip this pinfo if there are no DOFs on this node.
441  if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
442  continue;
443 
444  if (beginning_of_step)
445  {
446  if (_app.getExecutioner()->lastSolveConverged())
447  {
448  pinfo->_contact_force_old = pinfo->_contact_force;
449  pinfo->_accumulated_slip_old = pinfo->_accumulated_slip;
450  pinfo->_frictional_energy_old = pinfo->_frictional_energy;
451  pinfo->_mech_status_old = pinfo->_mech_status;
452  }
453  else if (pinfo->_mech_status_old == PenetrationInfo::MS_NO_CONTACT &&
454  pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
455  {
456  // The penetration info object could be based on a bad state so delete it
457  delete pinfo_pair.second;
458  pinfo_pair.second = NULL;
459  continue;
460  }
461 
462  pinfo->_locked_this_step = 0;
463  pinfo->_stick_locked_this_step = 0;
464  pinfo->_starting_elem = pinfo->_elem;
465  pinfo->_starting_side_num = pinfo->_side_num;
466  pinfo->_starting_closest_point_ref = pinfo->_closest_point_ref;
467  }
468  pinfo->_incremental_slip_prev_iter = pinfo->_incremental_slip;
469  }
470 }
MechanicalContactConstraint::_vars
std::vector< unsigned int > _vars
Definition: MechanicalContactConstraint.h:115

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

std::shared_ptr<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 90 of file MechanicalContactConstraint.h.

◆ _fe_problem

FEProblem& MechanicalContactConstraint::_fe_problem
protected

Definition at line 91 of file MechanicalContactConstraint.h.

Referenced by shouldApply().

◆ _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

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:
Generated on Mon Dec 10 2018 13:32:40 for www.mooseframework.org by   doxygen 1.8.14