LCOV - code coverage report
Current view: top level - src/userobjects - BilinearMixedModeCohesiveZoneModel.C (source / functions) Hit Total Coverage
Test: idaholab/moose contact: #31653 (1b668c) with base bb0a08 Lines: 246 257 95.7 %
Date: 2025-11-03 17:02:57 Functions: 28 29 96.6 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "BilinearMixedModeCohesiveZoneModel.h"
      11             : #include "MooseVariableFE.h"
      12             : #include "SystemBase.h"
      13             : #include "MortarUtils.h"
      14             : #include "MooseUtils.h"
      15             : #include "MathUtils.h"
      16             : 
      17             : #include "MortarContactUtils.h"
      18             : #include "FactorizedRankTwoTensor.h"
      19             : 
      20             : #include "ADReal.h"
      21             : #include "CohesiveZoneModelTools.h"
      22             : #include <Eigen/Core>
      23             : 
      24             : registerMooseObject("ContactApp", BilinearMixedModeCohesiveZoneModel);
      25             : 
      26             : InputParameters
      27         152 : BilinearMixedModeCohesiveZoneModel::validParams()
      28             : {
      29         152 :   InputParameters params = PenaltySimpleCohesiveZoneModel::validParams();
      30             : 
      31         152 :   params.addClassDescription("Computes the bilinear mixed mode cohesive zone model.");
      32         304 :   params.addRequiredCoupledVar("displacements",
      33             :                                "The string of displacements suitable for the problem statement");
      34             : 
      35             :   // Input parameters for bilinear mixed mode traction.
      36         304 :   params.addParam<MaterialPropertyName>("GI_c",
      37             :                                         "Critical energy release rate in normal direction.");
      38         304 :   params.addParam<MaterialPropertyName>("GII_c",
      39             :                                         "Critical energy release rate in shear direction.");
      40         304 :   params.addParam<MaterialPropertyName>("normal_strength", "Tensile strength in normal direction.");
      41         304 :   params.addParam<MaterialPropertyName>("shear_strength", "Tensile strength in shear direction.");
      42         304 :   params.addParam<Real>("power_law_parameter", "The power law parameter.");
      43         304 :   MooseEnum criterion("POWER_LAW BK", "BK");
      44         304 :   params.addParam<Real>("viscosity", 0.0, "Viscosity for damage model.");
      45         304 :   params.addParam<MooseEnum>(
      46             :       "mixed_mode_criterion", criterion, "Option for mixed mode propagation criterion.");
      47         304 :   params.addParam<bool>(
      48             :       "lag_displacement_jump",
      49         304 :       false,
      50             :       "Whether to use old displacement jumps to compute the effective displacement jump.");
      51         304 :   params.addParam<Real>(
      52         304 :       "regularization_alpha", 1e-10, "Regularization parameter for the Macaulay bracket.");
      53         304 :   params.addRangeCheckedParam<Real>(
      54             :       "penalty_stiffness", "penalty_stiffness > 0.0", "Penalty stiffness for CZM.");
      55         304 :   params.addParamNamesToGroup(
      56             :       "GI_c GII_c normal_strength shear_strength power_law_parameter viscosity "
      57             :       "mixed_mode_criterion lag_displacement_jump regularization_alpha "
      58             :       "penalty_stiffness",
      59             :       "Bilinear mixed mode traction");
      60             :   // End of input parameters for bilinear mixed mode traction.
      61             : 
      62             :   // Suppress augmented Lagrange parameters. AL implementation for CZM remains to be done.
      63         152 :   params.suppressParameter<Real>("max_penalty_multiplier");
      64         152 :   params.suppressParameter<Real>("penalty_multiplier");
      65         152 :   params.suppressParameter<Real>("penetration_tolerance");
      66         152 :   return params;
      67         152 : }
      68             : 
      69          76 : BilinearMixedModeCohesiveZoneModel::BilinearMixedModeCohesiveZoneModel(
      70          76 :     const InputParameters & parameters)
      71             :   : WeightedGapUserObject(parameters),
      72             :     PenaltyWeightedGapUserObject(parameters),
      73             :     WeightedVelocitiesUserObject(parameters),
      74             :     PenaltySimpleCohesiveZoneModel(parameters),
      75          76 :     _ndisp(coupledComponents("displacements")),
      76         152 :     _normal_strength(getMaterialProperty<Real>("normal_strength")),
      77         152 :     _shear_strength(getMaterialProperty<Real>("shear_strength")),
      78         152 :     _GI_c(getMaterialProperty<Real>("GI_c")),
      79         152 :     _GII_c(getMaterialProperty<Real>("GII_c")),
      80         152 :     _penalty_stiffness_czm(getParam<Real>("penalty_stiffness")),
      81         152 :     _mix_mode_criterion(getParam<MooseEnum>("mixed_mode_criterion").getEnum<MixedModeCriterion>()),
      82         152 :     _power_law_parameter(getParam<Real>("power_law_parameter")),
      83         152 :     _viscosity(getParam<Real>("viscosity")),
      84         304 :     _regularization_alpha(getParam<Real>("regularization_alpha"))
      85             : {
      86          76 :   _czm_interpolated_traction.resize(_ndisp);
      87             : 
      88         228 :   for (unsigned int i = 0; i < _ndisp; ++i)
      89             :   {
      90         152 :     _grad_disp.push_back(&adCoupledGradient("displacements", i));
      91         304 :     _grad_disp_neighbor.push_back(&adCoupledGradient("displacements", i));
      92             :   }
      93             : 
      94             :   // Set non-intervening components to zero
      95         152 :   for (unsigned int i = _ndisp; i < 3; i++)
      96             :   {
      97          76 :     _grad_disp.push_back(&_ad_grad_zero);
      98          76 :     _grad_disp_neighbor.push_back(&_ad_grad_zero);
      99             :   }
     100             : 
     101             :   // Checks for bilinear traction input.
     102         532 :   if (!isParamValid("GI_c") || !isParamValid("GII_c") || !isParamValid("normal_strength") ||
     103         532 :       !isParamValid("shear_strength") || !isParamValid("power_law_parameter") ||
     104         228 :       !isParamValid("penalty_stiffness"))
     105           0 :     paramError("GI_c",
     106             :                "The CZM bilinear mixed mode traction parameters GI_c, GII_c, normal_strength, "
     107             :                "shear_strength, and power_law_parameter are required. Revise your input and add "
     108             :                "those parameters if you want to use the bilinear mixed mode traction model. ");
     109          76 : }
     110             : 
     111             : void
     112       34158 : BilinearMixedModeCohesiveZoneModel::computeQpProperties()
     113             : {
     114       34158 :   WeightedVelocitiesUserObject::computeQpProperties();
     115             : 
     116             :   // Called after computeQpProperties() within the same algorithmic step.
     117             : 
     118             :   // Compute F and R.
     119       34158 :   const auto F = (ADRankTwoTensor::Identity() +
     120       68316 :                   ADRankTwoTensor::initializeFromRows(
     121       68316 :                       (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]));
     122       34158 :   const auto F_neighbor = (ADRankTwoTensor::Identity() +
     123       68316 :                            ADRankTwoTensor::initializeFromRows((*_grad_disp_neighbor[0])[_qp],
     124       34158 :                                                                (*_grad_disp_neighbor[1])[_qp],
     125       68316 :                                                                (*_grad_disp_neighbor[2])[_qp]));
     126             : 
     127             :   // TODO in follow-on PRs: Trim interior node variable derivatives
     128       34158 :   _F_interpolation = F * (_JxW_msm[_qp] * _coord[_qp]);
     129       34158 :   _F_neighbor_interpolation = F_neighbor * (_JxW_msm[_qp] * _coord[_qp]);
     130       34158 :   _normal_strength_interpolation = _normal_strength[_qp] * _JxW_msm[_qp] * _coord[_qp];
     131       34158 :   _shear_strength_interpolation = _shear_strength[_qp] * _JxW_msm[_qp] * _coord[_qp];
     132       34158 :   _GI_c_interpolation = _GI_c[_qp] * _JxW_msm[_qp] * _coord[_qp];
     133       34158 :   _GII_c_interpolation = _GII_c[_qp] * _JxW_msm[_qp] * _coord[_qp];
     134       34158 : }
     135             : 
     136             : void
     137       68316 : BilinearMixedModeCohesiveZoneModel::computeQpIProperties()
     138             : {
     139       68316 :   WeightedVelocitiesUserObject::computeQpIProperties();
     140             : 
     141             :   // Get the _dof_to_weighted_gap map
     142       68316 :   const auto * const dof = static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i));
     143             : 
     144             :   // TODO: Probably better to interpolate the deformation gradients.
     145      136632 :   _dof_to_F[dof] += (*_test)[_i][_qp] * _F_interpolation;
     146      136632 :   _dof_to_F_neighbor[dof] += (*_test)[_i][_qp] * _F_neighbor_interpolation;
     147      136632 :   _dof_to_normal_strength[dof] += (*_test)[_i][_qp] * _normal_strength_interpolation;
     148      136632 :   _dof_to_shear_strength[dof] += (*_test)[_i][_qp] * _shear_strength_interpolation;
     149             : 
     150      136632 :   _dof_to_GI_c[dof] += (*_test)[_i][_qp] * _GI_c_interpolation;
     151      136632 :   _dof_to_GII_c[dof] += (*_test)[_i][_qp] * _GII_c_interpolation;
     152       68316 : }
     153             : 
     154             : template <class T>
     155             : T
     156      409896 : BilinearMixedModeCohesiveZoneModel::normalizeQuantity(
     157             :     const std::unordered_map<const DofObject *, T> & map, const Node * const node)
     158             : {
     159      409896 :   return libmesh_map_find(map, node) / _dof_to_weighted_gap[node].second;
     160             : }
     161             : 
     162             : void
     163         634 : BilinearMixedModeCohesiveZoneModel::timestepSetup()
     164             : {
     165             :   // instead we call it explicitly here
     166         634 :   PenaltySimpleCohesiveZoneModel::timestepSetup();
     167             : 
     168             :   // save off tangential traction from the last timestep
     169        2235 :   for (auto & map_pr : _dof_to_damage)
     170             :   {
     171             :     auto & [damage, old_damage] = map_pr.second;
     172        1601 :     old_damage = {MetaPhysicL::raw_value(damage)};
     173        1601 :     damage = {0.0};
     174             :   }
     175             : 
     176         634 :   for (auto & [dof_object, delta_tangential_lm] : _dof_to_frictional_lagrange_multipliers)
     177             :     delta_tangential_lm.setZero();
     178         634 : }
     179             : 
     180             : void
     181        7680 : BilinearMixedModeCohesiveZoneModel::initialize()
     182             : {
     183             :   // instead we call it explicitly here
     184        7680 :   PenaltySimpleCohesiveZoneModel::initialize();
     185             : 
     186             :   // Avoid accumulating interpolation over the time step
     187             :   _dof_to_F.clear();
     188             :   _dof_to_F_neighbor.clear();
     189             :   _dof_to_mode_mixity_ratio.clear();
     190             :   _dof_to_normal_strength.clear();
     191             :   _dof_to_shear_strength.clear();
     192             :   _dof_to_GI_c.clear();
     193             :   _dof_to_GII_c.clear();
     194             :   _dof_to_delta_initial.clear();
     195             :   _dof_to_delta_final.clear();
     196             :   _dof_to_delta_max.clear();
     197             :   _dof_to_czm_traction.clear();
     198             :   _dof_to_interface_displacement_jump.clear();
     199        7680 :   _dof_to_czm_normal_traction.clear();
     200             :   _dof_to_interface_F.clear();
     201             :   _dof_to_interface_R.clear();
     202             : 
     203       24063 :   for (auto & map_pr : _dof_to_rotation_matrix)
     204       16383 :     map_pr.second.setToIdentity();
     205        7680 : }
     206             : 
     207             : void
     208       17079 : BilinearMixedModeCohesiveZoneModel::prepareJumpKinematicQuantities()
     209             : {
     210             :   // Compute rotation matrix from secondary surface
     211             :   // Rotation matrices and local interface displacement jump.
     212       51237 :   for (const auto i : make_range(_test->size()))
     213             :   {
     214       34158 :     const Node * const node = _lower_secondary_elem->node_ptr(i);
     215             : 
     216             :     // First call does not have maps available
     217       34158 :     const bool return_boolean = _dof_to_weighted_gap.find(node) == _dof_to_weighted_gap.end();
     218       34158 :     if (return_boolean)
     219           0 :       return;
     220             : 
     221       68316 :     _dof_to_rotation_matrix[node] = CohesiveZoneModelTools::computeReferenceRotation<true>(
     222       34158 :         _normals[i], _subproblem.mesh().dimension());
     223             : 
     224             :     // Every time we used quantities from a map we "denormalize" it from the mortar integral.
     225             :     // See normalizing member functions.
     226       34158 :     if (_dof_to_weighted_displacements.find(node) != _dof_to_weighted_displacements.end())
     227       68316 :       _dof_to_interface_displacement_jump[node] =
     228       34158 :           (libmesh_map_find(_dof_to_weighted_displacements, node));
     229             :   }
     230             : }
     231             : 
     232             : void
     233       34158 : BilinearMixedModeCohesiveZoneModel::computeFandR(const Node * const node)
     234             : {
     235             :   using std::isfinite;
     236             : 
     237             :   // First call does not have maps available
     238       34158 :   const bool return_boolean = _dof_to_F.find(node) == _dof_to_F.end();
     239       34158 :   if (return_boolean)
     240           0 :     return;
     241             : 
     242       34158 :   const auto normalized_F = normalizeQuantity(_dof_to_F, node);
     243       34158 :   const auto normalized_F_neighbor = normalizeQuantity(_dof_to_F_neighbor, node);
     244             : 
     245             :   // This 'averaging' assumption below can probably be improved upon.
     246      102474 :   _dof_to_interface_F[node] = 0.5 * (normalized_F + normalized_F_neighbor);
     247             : 
     248      136632 :   for (const auto i : make_range(3))
     249      409896 :     for (const auto j : make_range(3))
     250      307422 :       if (!isfinite(MetaPhysicL::raw_value(normalized_F(i, j))))
     251             :         throw MooseException(
     252             :             "The deformation gradient on the secondary surface is not finite in "
     253           0 :             "PenaltySimpleCohesiveZoneModel. MOOSE needs to cut the time step size.");
     254             : 
     255       34158 :   const auto dof_to_interface_F_node = libmesh_map_find(_dof_to_interface_F, node);
     256             : 
     257       34158 :   const ADFactorizedRankTwoTensor C = dof_to_interface_F_node.transpose() * dof_to_interface_F_node;
     258       34158 :   const auto Uinv = MathUtils::sqrt(C).inverse().get();
     259       68316 :   _dof_to_interface_R[node] = dof_to_interface_F_node * Uinv;
     260             : 
     261             :   // Transform interface jump according to two rotation matrices
     262       34158 :   const auto global_interface_displacement = _dof_to_interface_displacement_jump[node];
     263       68316 :   _dof_to_interface_displacement_jump[node] =
     264       34158 :       (_dof_to_interface_R[node] * _dof_to_rotation_matrix[node]).transpose() *
     265       34158 :       global_interface_displacement;
     266             : }
     267             : 
     268             : void
     269       34158 : BilinearMixedModeCohesiveZoneModel::computeBilinearMixedModeTraction(const Node * const node)
     270             : {
     271             :   using std::max, std::min;
     272             : 
     273             :   // First call does not have maps available
     274       34158 :   const bool return_boolean = _dof_to_weighted_gap.find(node) == _dof_to_weighted_gap.end();
     275       34158 :   if (return_boolean)
     276           0 :     return;
     277             : 
     278       34158 :   computeModeMixity(node);
     279       34158 :   computeCriticalDisplacementJump(node);
     280       34158 :   computeFinalDisplacementJump(node);
     281       34158 :   computeEffectiveDisplacementJump(node);
     282       34158 :   computeDamage(node);
     283             : 
     284             :   // Split displacement jump into active and inactive parts
     285             :   const auto interface_displacement_jump =
     286       34158 :       normalizeQuantity(_dof_to_interface_displacement_jump, node);
     287             : 
     288       34158 :   const ADRealVectorValue delta_active(max(interface_displacement_jump(0), 0.0),
     289             :                                        interface_displacement_jump(1),
     290             :                                        interface_displacement_jump(2));
     291       34158 :   const ADRealVectorValue delta_inactive(min(interface_displacement_jump(0), 0.0), 0.0, 0.0);
     292             : 
     293             :   // This traction vector is local at this point.
     294       68316 :   _dof_to_czm_traction[node] =
     295      102474 :       -(1.0 - _dof_to_damage[node].first) * _penalty_stiffness_czm * delta_active -
     296       68316 :       _penalty_stiffness_czm * delta_inactive;
     297             : }
     298             : 
     299             : void
     300       34158 : BilinearMixedModeCohesiveZoneModel::computeGlobalTraction(const Node * const node)
     301             : {
     302             :   // First call does not have maps available
     303       34158 :   const bool return_boolean = _dof_to_czm_traction.find(node) == _dof_to_czm_traction.end();
     304       34158 :   if (return_boolean)
     305           0 :     return;
     306             : 
     307       34158 :   const auto local_traction_vector = libmesh_map_find(_dof_to_czm_traction, node);
     308       34158 :   const auto rotation_matrix = libmesh_map_find(_dof_to_rotation_matrix, node);
     309             : 
     310       68316 :   _dof_to_czm_traction[node] = rotation_matrix * local_traction_vector;
     311             : }
     312             : 
     313             : void
     314       34158 : BilinearMixedModeCohesiveZoneModel::computeModeMixity(const Node * const node)
     315             : {
     316             :   using std::sqrt;
     317             : 
     318             :   const auto interface_displacement_jump =
     319       34158 :       normalizeQuantity(_dof_to_interface_displacement_jump, node);
     320             : 
     321       34158 :   if (interface_displacement_jump(0) > _epsilon_tolerance)
     322             :   {
     323             :     const auto delta_s =
     324       30252 :         sqrt(interface_displacement_jump(1) * interface_displacement_jump(1) +
     325       30252 :              interface_displacement_jump(2) * interface_displacement_jump(2) + _epsilon_tolerance);
     326             : 
     327       30252 :     _dof_to_mode_mixity_ratio[node] = delta_s / interface_displacement_jump(0);
     328             :   }
     329             :   else
     330        3906 :     _dof_to_mode_mixity_ratio[node] = 0;
     331       34158 : }
     332             : 
     333             : void
     334       34158 : BilinearMixedModeCohesiveZoneModel::computeCriticalDisplacementJump(const Node * const node)
     335             : {
     336             :   using std::sqrt;
     337             : 
     338             :   const auto interface_displacement_jump =
     339       34158 :       normalizeQuantity(_dof_to_interface_displacement_jump, node);
     340             : 
     341       34158 :   const auto mixity_ratio = libmesh_map_find(_dof_to_mode_mixity_ratio, node);
     342             : 
     343             :   const auto delta_normal_knot =
     344       34158 :       normalizeQuantity(_dof_to_normal_strength, node) / _penalty_stiffness_czm;
     345             :   const auto delta_shear_knot =
     346       34158 :       normalizeQuantity(_dof_to_shear_strength, node) / _penalty_stiffness_czm;
     347             : 
     348       34158 :   _dof_to_delta_initial[node] = delta_shear_knot;
     349             : 
     350       34158 :   if (interface_displacement_jump(0) > _epsilon_tolerance)
     351             :   {
     352       30252 :     const auto delta_mixed = sqrt(delta_shear_knot * delta_shear_knot +
     353       30252 :                                   Utility::pow<2>(mixity_ratio * delta_normal_knot));
     354             : 
     355       60504 :     _dof_to_delta_initial[node] = delta_normal_knot * delta_shear_knot *
     356       90756 :                                   sqrt(1.0 + mixity_ratio * mixity_ratio) / delta_mixed;
     357             :   }
     358       34158 : }
     359             : 
     360             : void
     361       34158 : BilinearMixedModeCohesiveZoneModel::computeFinalDisplacementJump(const Node * const node)
     362             : {
     363             :   using std::sqrt, std::pow;
     364             : 
     365             :   const auto interface_displacement_jump =
     366       34158 :       normalizeQuantity(_dof_to_interface_displacement_jump, node);
     367             : 
     368       34158 :   const auto mixity_ratio = libmesh_map_find(_dof_to_mode_mixity_ratio, node);
     369             : 
     370       34158 :   const auto normalized_GI_c = normalizeQuantity(_dof_to_GI_c, node);
     371       34158 :   const auto normalized_GII_c = normalizeQuantity(_dof_to_GII_c, node);
     372             : 
     373       68316 :   _dof_to_delta_final[node] =
     374       68316 :       sqrt(2.0) * 2.0 * normalized_GII_c / normalizeQuantity(_dof_to_shear_strength, node);
     375             : 
     376       34158 :   if (interface_displacement_jump(0) > _epsilon_tolerance)
     377             :   {
     378       30252 :     if (_mix_mode_criterion == MixedModeCriterion::BK)
     379             :     {
     380       21264 :       _dof_to_delta_final[node] =
     381       10632 :           2.0 / _penalty_stiffness_czm / libmesh_map_find(_dof_to_delta_initial, node) *
     382       10632 :           (normalized_GI_c +
     383       10632 :            (normalized_GII_c - normalized_GI_c) *
     384       31896 :                pow(mixity_ratio * mixity_ratio / (1 + mixity_ratio * mixity_ratio),
     385       21264 :                    _power_law_parameter));
     386             :     }
     387       19620 :     else if (_mix_mode_criterion == MixedModeCriterion::POWER_LAW)
     388             :     {
     389             :       const auto Gc_mixed =
     390       39240 :           pow(1.0 / normalized_GI_c, _power_law_parameter) +
     391       39240 :           pow(mixity_ratio * mixity_ratio / normalized_GII_c, _power_law_parameter);
     392       58860 :       _dof_to_delta_final[node] = (2.0 + 2.0 * mixity_ratio * mixity_ratio) /
     393       19620 :                                   _penalty_stiffness_czm /
     394       19620 :                                   libmesh_map_find(_dof_to_delta_initial, node) *
     395       39240 :                                   pow(Gc_mixed, -1.0 / _power_law_parameter);
     396             :     }
     397             :   }
     398       34158 : }
     399             : 
     400             : void
     401       34158 : BilinearMixedModeCohesiveZoneModel::computeEffectiveDisplacementJump(const Node * const node)
     402             : {
     403             :   using std::sqrt;
     404             : 
     405             :   const auto interface_displacement_jump =
     406       34158 :       normalizeQuantity(_dof_to_interface_displacement_jump, node);
     407             : 
     408             :   const auto delta_normal_pos =
     409       34158 :       MathUtils::regularizedHeavyside(interface_displacement_jump(0), _regularization_alpha) *
     410             :       interface_displacement_jump(0);
     411             : 
     412       34158 :   _dof_to_delta_max[node] = sqrt(Utility::pow<2>(interface_displacement_jump(1)) +
     413       34158 :                                  Utility::pow<2>(interface_displacement_jump(2)) +
     414       68316 :                                  Utility::pow<2>(delta_normal_pos) + _epsilon_tolerance);
     415       34158 : }
     416             : 
     417             : void
     418       34158 : BilinearMixedModeCohesiveZoneModel::computeDamage(const Node * const node)
     419             : {
     420       34158 :   const auto delta_max = libmesh_map_find(_dof_to_delta_max, node);
     421       34158 :   const auto delta_initial = libmesh_map_find(_dof_to_delta_initial, node);
     422       34158 :   const auto delta_final = libmesh_map_find(_dof_to_delta_final, node);
     423             : 
     424       34158 :   if (delta_max < delta_initial)
     425       34158 :     _dof_to_damage[node].first = 0;
     426           0 :   else if (delta_max > delta_final)
     427           0 :     _dof_to_damage[node].first = 1.0;
     428             :   else
     429           0 :     _dof_to_damage[node].first =
     430           0 :         delta_final * (delta_max - delta_initial) / delta_max / (delta_final - delta_initial);
     431             : 
     432       34158 :   if (_dof_to_damage[node].first < _dof_to_damage[node].second)
     433             :     // Irreversibility
     434           0 :     _dof_to_damage[node].first = _dof_to_damage[node].second;
     435             : 
     436             :   // Viscous regularization
     437       34158 :   _dof_to_damage[node].first =
     438       34158 :       (_dof_to_damage[node].first + _viscosity * _dof_to_damage[node].second / _dt) /
     439       68316 :       (_viscosity / _dt + 1.0);
     440       34158 : }
     441             : 
     442             : void
     443       17079 : BilinearMixedModeCohesiveZoneModel::reinit()
     444             : {
     445             :   // Normal contact pressure with penalty
     446       17079 :   PenaltySimpleCohesiveZoneModel::reinit();
     447             : 
     448       51237 :   for (const auto i : index_range(_czm_interpolated_traction))
     449             :   {
     450       34158 :     _czm_interpolated_traction[i].resize(_qrule_msm->n_points());
     451             : 
     452      102474 :     for (const auto qp : make_range(_qrule_msm->n_points()))
     453       68316 :       _czm_interpolated_traction[i][qp] = 0.0;
     454             :   }
     455             : 
     456             :   // iterate over nodes
     457       51237 :   for (const auto i : make_range(_test->size()))
     458             :   {
     459             :     // current node
     460       34158 :     const Node * const node = _lower_secondary_elem->node_ptr(i);
     461             : 
     462             :     // End of CZM bilinear computations
     463       34158 :     if (_dof_to_czm_traction.find(node) == _dof_to_czm_traction.end())
     464             :       return;
     465             : 
     466       34158 :     const auto & test_i = (*_test)[i];
     467      102474 :     for (const auto qp : make_range(_qrule_msm->n_points()))
     468      204948 :       for (const auto idx : index_range(_czm_interpolated_traction))
     469      273264 :         _czm_interpolated_traction[idx][qp] += test_i[qp] * _dof_to_czm_traction[node](idx);
     470             :   }
     471             : }
     472             : 
     473             : void
     474        7680 : BilinearMixedModeCohesiveZoneModel::finalize()
     475             : {
     476        7680 :   PenaltySimpleCohesiveZoneModel::finalize();
     477             : 
     478        7680 :   const bool send_data_back = !constrainedByOwner();
     479             : 
     480       15360 :   Moose::Mortar::Contact::communicateR2T(
     481        7680 :       _dof_to_F, _subproblem.mesh(), _nodal, _communicator, send_data_back);
     482             : 
     483       15360 :   Moose::Mortar::Contact::communicateR2T(
     484        7680 :       _dof_to_F_neighbor, _subproblem.mesh(), _nodal, _communicator, send_data_back);
     485             : 
     486       15360 :   Moose::Mortar::Contact::communicateRealObject(
     487        7680 :       _dof_to_normal_strength, _subproblem.mesh(), _nodal, _communicator, send_data_back);
     488             : 
     489       15360 :   Moose::Mortar::Contact::communicateRealObject(
     490        7680 :       _dof_to_shear_strength, _subproblem.mesh(), _nodal, _communicator, send_data_back);
     491             : 
     492       15360 :   Moose::Mortar::Contact::communicateRealObject(
     493        7680 :       _dof_to_GI_c, _subproblem.mesh(), _nodal, _communicator, send_data_back);
     494             : 
     495       15360 :   Moose::Mortar::Contact::communicateRealObject(
     496        7680 :       _dof_to_GII_c, _subproblem.mesh(), _nodal, _communicator, send_data_back);
     497             : 
     498       15360 :   Moose::Mortar::Contact::communicateRealObject(
     499        7680 :       _dof_to_weighted_displacements, _subproblem.mesh(), _nodal, _communicator, send_data_back);
     500        7680 : }
     501             : 
     502             : const ADVariableValue &
     503      273264 : BilinearMixedModeCohesiveZoneModel::czmGlobalTraction(const unsigned int i) const
     504             : {
     505             :   mooseAssert(i <= 3,
     506             :               "Internal error in czmGlobalTraction. Index exceeds the traction vector size.");
     507             : 
     508      273264 :   return _czm_interpolated_traction[i];
     509             : }
     510             : 
     511             : Real
     512        1086 : BilinearMixedModeCohesiveZoneModel::getModeMixityRatio(const Node * const node) const
     513             : {
     514        1086 :   const auto it = _dof_to_mode_mixity_ratio.find(_subproblem.mesh().nodePtr(node->id()));
     515             : 
     516        1086 :   if (it != _dof_to_mode_mixity_ratio.end())
     517        1086 :     return MetaPhysicL::raw_value(it->second);
     518             :   else
     519             :     return 0.0;
     520             : }
     521             : 
     522             : Real
     523        1086 : BilinearMixedModeCohesiveZoneModel::getCohesiveDamage(const Node * const node) const
     524             : {
     525        1086 :   const auto it = _dof_to_damage.find(_subproblem.mesh().nodePtr(node->id()));
     526             : 
     527        1086 :   if (it != _dof_to_damage.end())
     528        1086 :     return MetaPhysicL::raw_value(it->second.first);
     529             :   else
     530             :     return 0.0;
     531             : }
     532             : 
     533             : Real
     534        1086 : BilinearMixedModeCohesiveZoneModel::getLocalDisplacementNormal(const Node * const node) const
     535             : {
     536        1086 :   const auto it = _dof_to_interface_displacement_jump.find(_subproblem.mesh().nodePtr(node->id()));
     537        1086 :   const auto it2 = _dof_to_weighted_gap.find(_subproblem.mesh().nodePtr(node->id()));
     538             : 
     539        1086 :   if (it != _dof_to_interface_displacement_jump.end() && it2 != _dof_to_weighted_gap.end())
     540        2172 :     return MetaPhysicL::raw_value(it->second(0) / it2->second.second);
     541             :   else
     542             :     return 0.0;
     543             : }
     544             : 
     545             : Real
     546        1086 : BilinearMixedModeCohesiveZoneModel::getLocalDisplacementTangential(const Node * const node) const
     547             : {
     548        1086 :   const auto it = _dof_to_interface_displacement_jump.find(_subproblem.mesh().nodePtr(node->id()));
     549        1086 :   const auto it2 = _dof_to_weighted_gap.find(_subproblem.mesh().nodePtr(node->id()));
     550             : 
     551        1086 :   if (it != _dof_to_interface_displacement_jump.end() && it2 != _dof_to_weighted_gap.end())
     552        2172 :     return MetaPhysicL::raw_value(it->second(1) / it2->second.second);
     553             :   else
     554             :     return 0.0;
     555             : }

Generated by: LCOV version 1.14