LCOV - code coverage report
Current view: top level - src/constraints - RebarBondSlipConstraint.C (source / functions) Hit Total Coverage
Test: idaholab/blackbear: 9ef4b1 Lines: 227 248 91.5 %
Date: 2025-10-02 23:14:06 Functions: 21 26 80.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/primary/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : // MOOSE includes
      11             : #include "MooseArray.h"
      12             : #include "MooseError.h"
      13             : #include "MooseTypes.h"
      14             : #include "RebarBondSlipConstraint.h"
      15             : #include "Assembly.h"
      16             : #include "SystemBase.h"
      17             : #include "FEProblem.h"
      18             : #include "MathUtils.h"
      19             : #include "MooseTypes.h"
      20             : #include "Coupleable.h"
      21             : #include "Executioner.h"
      22             : 
      23             : #include "metaphysicl/raw_type.h"
      24             : #include "libmesh/utility.h"
      25             : 
      26             : registerMooseObject("BlackBearApp", RebarBondSlipConstraint);
      27             : registerMooseObject("BlackBearApp", ADRebarBondSlipConstraint);
      28             : 
      29             : template <bool is_ad>
      30             : InputParameters
      31         176 : RebarBondSlipConstraintTempl<is_ad>::validParams()
      32             : {
      33         176 :   InputParameters params = EqualValueEmbeddedConstraint::validParams();
      34         176 :   params.addClassDescription("Constraint using a physical model to enforce the bond-slip behavior "
      35             :                              "between concrete and rebar.");
      36         352 :   params.addRequiredCoupledVar(
      37             :       "displacements",
      38             :       "The displacements appropriate for the simulation geometry and coordinate system");
      39         352 :   params.addRequiredParam<Real>("max_bond_stress", "Maximum bond stress");
      40         352 :   params.addRequiredParam<Real>("transitional_slip_value",
      41             :                                 "Transition between loading and frictional slip");
      42         352 :   params.addRequiredParam<Real>("rebar_radius", "Radius of the rebar");
      43         352 :   MooseEnum bondslip_model("concrete_rebar_model elastic_perfect_plastic_model",
      44             :                            "elastic_perfect_plastic_model");
      45         352 :   params.addParam<MooseEnum>("bondslip_model", bondslip_model, "Node-Element bond slip model.");
      46             : 
      47         352 :   params.addCoupledVar("output_axial_slip", "Rebar axial slip output variable.");
      48         352 :   params.addCoupledVar("output_axial_force", "Rebar axial force output variable");
      49         352 :   params.addCoupledVar("output_axial_plastic_slip", "Rebar axial plastic component of bond slip");
      50         176 :   return params;
      51         176 : }
      52             : 
      53             : template <bool is_ad>
      54          88 : RebarBondSlipConstraintTempl<is_ad>::RebarBondSlipConstraintTempl(
      55             :     const InputParameters & parameters)
      56             :   : EqualValueEmbeddedConstraintTempl<is_ad>(parameters),
      57          88 :     _bondslip_model(
      58          88 :         this->template getParam<MooseEnum>("bondslip_model").template getEnum<BondSlipModel>()),
      59          88 :     _component(libMesh::invalid_uint),
      60          88 :     _mesh_dimension(_mesh.dimension()),
      61          88 :     _disp_vars_nums(_mesh_dimension, libMesh::invalid_uint),
      62          88 :     _disp_vars(_mesh_dimension, nullptr),
      63         176 :     _max_bond_stress(this->template getParam<Real>("max_bond_stress")),
      64         176 :     _bar_radius(this->template getParam<Real>("rebar_radius")),
      65         264 :     _transitional_slip(this->template getParam<Real>("transitional_slip_value"))
      66             : {
      67         176 :   if (_mesh_dimension != this->coupledComponents("displacements"))
      68           0 :     mooseError("In RebarBondSlipConstraint, number of displacements must equal the mesh dimension");
      69             : 
      70         264 :   for (unsigned int i = 0; i < _mesh_dimension; ++i)
      71             :   {
      72         176 :     _disp_vars_nums[i] = this->coupled("displacements", i);
      73         176 :     _disp_vars[i] = this->getVar("displacements", i);
      74             :   }
      75             : 
      76         264 :   for (std::size_t i = 0; i < _disp_vars_nums.size(); ++i)
      77         176 :     if (_var.number() == _disp_vars_nums[i])
      78          88 :       _component = i;
      79          88 :   if (_component == libMesh::invalid_uint)
      80           0 :     mooseError("Problem with displacements in " + this->_name);
      81             : 
      82         176 :   if (this->isParamValid("output_axial_slip"))
      83          88 :     _output_axial_slip = &(this->writableVariable("output_axial_slip"));
      84         176 :   if (this->isParamValid("output_axial_force"))
      85          88 :     _output_axial_force = &(this->writableVariable("output_axial_force"));
      86         176 :   if (this->isParamValid("output_axial_plastic_slip"))
      87          88 :     _output_axial_plastic_slip = &(this->writableVariable("output_axial_plastic_slip"));
      88          88 : }
      89             : 
      90             : template <bool is_ad>
      91             : void
      92          88 : RebarBondSlipConstraintTempl<is_ad>::initialSetup()
      93             : {
      94             : 
      95        1584 :   for (auto it = _secondary_to_primary_map.begin(); it != _secondary_to_primary_map.end(); ++it)
      96        1496 :     if (_bondslip.find(it->first) == _bondslip.end())
      97        1496 :       _bondslip.insert(std::pair<dof_id_type, bondSlipData>(it->first, bondSlipData()));
      98          88 : }
      99             : 
     100             : template <bool is_ad>
     101             : void
     102        3840 : RebarBondSlipConstraintTempl<is_ad>::timestepSetup()
     103             : {
     104       69120 :   for (auto iter = _secondary_to_primary_map.begin(); iter != _secondary_to_primary_map.end();
     105             :        ++iter)
     106             :   {
     107       65280 :     dof_id_type node_id = iter->first;
     108             :     mooseAssert(_bondslip.find(node_id) != _bondslip.end(), "Node not found in bond-slip map");
     109       65280 :     if (this->_app.getExecutioner()->lastSolveConverged())
     110             :     {
     111       65280 :       _bondslip[node_id].slip_min_old = _bondslip[node_id].slip_min;
     112       65280 :       _bondslip[node_id].slip_max_old = _bondslip[node_id].slip_max;
     113       65280 :       _bondslip[node_id].bond_stress_min_old = _bondslip[node_id].bond_stress_min;
     114       65280 :       _bondslip[node_id].bond_stress_max_old = _bondslip[node_id].bond_stress_max;
     115             :     }
     116             :   }
     117        3840 : }
     118             : 
     119             : template <bool is_ad>
     120             : void
     121      396576 : RebarBondSlipConstraintTempl<is_ad>::computeTangent()
     122             : {
     123             :   _secondary_node_tangent = 0.0;
     124      396576 :   _secondary_node_length = 0.0;
     125             :   // get normals
     126             :   // get connected elements of the current node
     127      396576 :   const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map = _mesh.nodeToElemMap();
     128      396576 :   auto node_to_elem_pair = node_to_elem_map.find(_current_node->id());
     129             :   mooseAssert(node_to_elem_pair != node_to_elem_map.end(), "Missing entry in node to elem map");
     130             :   const std::vector<dof_id_type> & elems = node_to_elem_pair->second;
     131             :   mooseAssert(_var.currentElem()->dim() == 1, "Elements attached to secondary node must be 1D.");
     132     1143072 :   for (auto & elem : elems)
     133             :   {
     134      746496 :     const Elem * elem_ptr = _mesh.elemPtr(elem);
     135      746496 :     this->_assembly.reinit(elem_ptr);
     136      746496 :     _secondary_node_length += elem_ptr->volume();
     137             :     const std::vector<RealGradient> * tangents =
     138      746496 :         &_subproblem.assembly(Constraint::_tid, _sys.number()).getFE(FEType(), 1)->get_dxyzdxi();
     139             :     const std::vector<Real> * JxW =
     140      746496 :         &_subproblem.assembly(Constraint::_tid, _sys.number()).getFE(FEType(), 1)->get_JxW();
     141     2239488 :     for (std::size_t i = 0; i < tangents->size(); i++)
     142             :       _secondary_node_tangent += (*tangents)[i] * (*JxW)[i];
     143             :   }
     144             : 
     145      563448 :   _secondary_node_tangent = _secondary_node_tangent / _secondary_node_tangent.norm();
     146      563448 :   _secondary_node_length = _secondary_node_length / elems.size();
     147      396576 : }
     148             : 
     149             : template <>
     150             : RealVectorValue
     151      229704 : RebarBondSlipConstraintTempl<false>::computeRelativeDisplacement()
     152             : {
     153             :   RealVectorValue relative_disp;
     154      689112 :   for (unsigned int i = 0; i < _mesh_dimension; ++i)
     155      459408 :     relative_disp(i) = ((_disp_vars[i]->dofValues())[0] - (_disp_vars[i]->slnNeighbor())[0]);
     156             : 
     157      229704 :   return relative_disp;
     158             : }
     159             : 
     160             : template <>
     161             : ADRealVectorValue
     162      166872 : RebarBondSlipConstraintTempl<true>::computeRelativeDisplacement()
     163             : {
     164             :   ADRealVectorValue relative_disp;
     165      500616 :   for (unsigned int i = 0; i < _mesh_dimension; ++i)
     166      667488 :     relative_disp(i) = ((_disp_vars[i]->adDofValues())[0] - (_disp_vars[i]->adSlnNeighbor())[0]);
     167      166872 :   return relative_disp;
     168             : }
     169             : 
     170             : template <bool is_ad>
     171             : void
     172      396576 : RebarBondSlipConstraintTempl<is_ad>::reinitConstraint()
     173             : {
     174      396576 :   computeTangent();
     175             : 
     176      396576 :   GenericRealVectorValue<is_ad> relative_disp = computeRelativeDisplacement();
     177             : 
     178      166872 :   GenericReal<is_ad> slip = relative_disp * _secondary_node_tangent;
     179      166872 :   GenericRealVectorValue<is_ad> slip_axial = slip * _secondary_node_tangent;
     180      166872 :   GenericRealVectorValue<is_ad> slip_normal = relative_disp - slip_axial;
     181             : 
     182      396576 :   const Node * node = _current_node;
     183      396576 :   auto it = _bondslip.find(node->id());
     184             :   mooseAssert(it != _bondslip.end(), "Node not found in bond-slip map");
     185      396576 :   bondSlipData * bond_slip = &(it->second);
     186             : 
     187      396576 :   bond_slip->slip_min = std::min(bond_slip->slip_min_old, MetaPhysicL::raw_value(slip));
     188      396576 :   bond_slip->slip_max = std::max(bond_slip->slip_max_old, MetaPhysicL::raw_value(slip));
     189             : 
     190      396576 :   _bond_stress_deriv = 0.0;
     191      396576 :   _bond_stress = 0.0;
     192             : 
     193      396576 :   Real plastic_slip = 0;
     194      396576 :   switch (_bondslip_model)
     195             :   {
     196      328032 :     case BondSlipModel::CONCRETE_REBAR_MODEL:
     197      328032 :       concreteRebarModel(slip, bond_slip, _bond_stress, _bond_stress_deriv, plastic_slip);
     198      166872 :       break;
     199             : 
     200       68544 :     case BondSlipModel::ELASTIC_PERFECT_PLASTIC_MODEL:
     201       68544 :       elasticPerfectPlasticModel(slip, bond_slip, _bond_stress, _bond_stress_deriv, plastic_slip);
     202           0 :       break;
     203             : 
     204           0 :     default:
     205           0 :       mooseError("Invalid bond slip model");
     206             :       break;
     207             :   }
     208             : 
     209      166872 :   GenericReal<is_ad> bond_force =
     210      396576 :       2.0 * libMesh::pi * _bar_radius * _secondary_node_length * _bond_stress;
     211      166872 :   GenericReal<is_ad> bond_force_deriv =
     212      563448 :       2.0 * libMesh::pi * _bar_radius * _secondary_node_length * _bond_stress_deriv;
     213             : 
     214      166872 :   GenericRealVectorValue<is_ad> constraint_force_axial = bond_force * _secondary_node_tangent;
     215      166872 :   GenericRealVectorValue<is_ad> constraint_force_normal = this->_penalty * slip_normal;
     216             : 
     217      563448 :   _constraint_residual = constraint_force_axial + constraint_force_normal;
     218      396576 :   _constraint_jacobian_axial = bond_force_deriv * _secondary_node_tangent;
     219             : 
     220      396576 :   bond_slip->bond_stress_min =
     221      513822 :       std::min(bond_slip->bond_stress_min_old, MetaPhysicL::raw_value(_bond_stress));
     222      396576 :   bond_slip->bond_stress_max =
     223      396576 :       std::max(bond_slip->bond_stress_max_old, MetaPhysicL::raw_value(_bond_stress));
     224             : 
     225      396576 :   if (_output_axial_slip)
     226      396576 :     _output_axial_slip->setNodalValue(MetaPhysicL::raw_value(slip));
     227      396576 :   if (_output_axial_force)
     228      396576 :     _output_axial_force->setNodalValue(MetaPhysicL::raw_value(_bond_stress));
     229      396576 :   if (_output_axial_plastic_slip)
     230      396576 :     _output_axial_plastic_slip->setNodalValue(plastic_slip);
     231      396576 : }
     232             : 
     233             : template <bool is_ad>
     234             : GenericReal<is_ad>
     235     1347420 : RebarBondSlipConstraintTempl<is_ad>::computeQpResidual(Moose::ConstraintType type)
     236             : {
     237     1347420 :   GenericReal<is_ad> resid = _constraint_residual(_component);
     238     1347420 :   switch (type)
     239             :   {
     240      269484 :     case Moose::Secondary:
     241      269484 :       return resid * _test_secondary[_i][_qp];
     242             : 
     243     1077936 :     case Moose::Primary:
     244     1647504 :       return -resid * _test_primary[_i][_qp];
     245             :   }
     246             : 
     247           0 :   return 0.0;
     248             : }
     249             : 
     250             : template <>
     251             : Real
     252           0 : RebarBondSlipConstraintTempl<true>::computeQpJacobian(Moose::ConstraintJacobianType /*type*/)
     253             : {
     254           0 :   mooseError("In ADRebarBondSlipConstraint, computeQpJacobian should not be called. Check "
     255             :              "computeQpJacobian implementation.");
     256             :   return 0.0;
     257             : }
     258             : 
     259             : template <>
     260             : Real
     261     2688660 : RebarBondSlipConstraintTempl<false>::computeQpJacobian(Moose::ConstraintJacobianType type)
     262             : {
     263             : 
     264     2688660 :   Real jac_axial = _constraint_jacobian_axial(_component);
     265             : 
     266     2688660 :   switch (type)
     267             :   {
     268      225204 :     case Moose::SecondarySecondary:
     269      225204 :       return _phi_secondary[_j][_qp] * jac_axial * _secondary_node_tangent(_component) *
     270      225204 :                  _test_secondary[_i][_qp] +
     271      225204 :              _phi_secondary[_j][_qp] * this->_penalty * _test_secondary[_i][_qp] *
     272      225204 :                  (1.0 - _secondary_node_tangent(_component) * _secondary_node_tangent(_component));
     273             : 
     274      312528 :     case Moose::SecondaryPrimary:
     275      312528 :       return -_phi_primary[_j][_qp] * jac_axial * _secondary_node_tangent(_component) *
     276      312528 :                  _test_secondary[_i][_qp] -
     277      312528 :              _phi_primary[_j][_qp] * this->_penalty * _test_secondary[_i][_qp] *
     278      312528 :                  (1.0 - _secondary_node_tangent(_component) * _secondary_node_tangent(_component));
     279             : 
     280      900816 :     case Moose::PrimarySecondary:
     281     1801632 :       return -_test_primary[_i][_qp] * jac_axial * _secondary_node_tangent(_component) *
     282      900816 :                  _phi_secondary[_j][_qp] -
     283      900816 :              _test_primary[_i][_qp] * this->_penalty * _phi_secondary[_j][_qp] *
     284      900816 :                  (1.0 - _secondary_node_tangent(_component) * _secondary_node_tangent(_component));
     285             : 
     286     1250112 :     case Moose::PrimaryPrimary:
     287     1250112 :       return _test_primary[_i][_qp] * jac_axial * _secondary_node_tangent(_component) *
     288     1250112 :                  _phi_primary[_j][_qp] +
     289     1250112 :              _test_primary[_i][_qp] * this->_penalty * _phi_primary[_j][_qp] *
     290     1250112 :                  (1.0 - _secondary_node_tangent(_component) * _secondary_node_tangent(_component));
     291             : 
     292           0 :     default:
     293           0 :       mooseError("Unsupported type");
     294             :       break;
     295             :   }
     296             :   return 0.0;
     297             : }
     298             : 
     299             : template <>
     300             : Real
     301           0 : RebarBondSlipConstraintTempl<true>::computeQpOffDiagJacobian(Moose::ConstraintJacobianType /*type*/,
     302             :                                                              unsigned int /*jvar*/)
     303             : {
     304           0 :   mooseError("In ADRebarBondSlipConstraint, computeQpOffDiagJacobian should not be called. Check "
     305             :              "computeQpJacobian implementation.");
     306             :   return 0.0;
     307             : }
     308             : 
     309             : template <>
     310             : Real
     311     2688660 : RebarBondSlipConstraintTempl<false>::computeQpOffDiagJacobian(Moose::ConstraintJacobianType type,
     312             :                                                               unsigned int jvar)
     313             : {
     314     2688660 :   Real jac_axial = _constraint_jacobian_axial(_component);
     315             : 
     316             :   unsigned int coupled_component;
     317             :   Real tangent_component_in_coupled_var_dir = 1.0;
     318             :   if (getCoupledVarComponent(jvar, coupled_component))
     319     2688660 :     tangent_component_in_coupled_var_dir = _secondary_node_tangent(coupled_component);
     320             : 
     321     2688660 :   switch (type)
     322             :   {
     323      225204 :     case Moose::SecondarySecondary:
     324      225204 :       return _phi_secondary[_j][_qp] * jac_axial * tangent_component_in_coupled_var_dir *
     325      225204 :                  _test_secondary[_i][_qp] -
     326      225204 :              _phi_secondary[_j][_qp] * this->_penalty * _test_secondary[_i][_qp] *
     327      225204 :                  _secondary_node_tangent(_component) * tangent_component_in_coupled_var_dir;
     328             : 
     329      312528 :     case Moose::SecondaryPrimary:
     330      312528 :       return -_phi_primary[_j][_qp] * jac_axial * tangent_component_in_coupled_var_dir *
     331      312528 :                  _test_secondary[_i][_qp] +
     332      312528 :              _phi_primary[_j][_qp] * this->_penalty * _test_secondary[_i][_qp] *
     333      312528 :                  _secondary_node_tangent(_component) * tangent_component_in_coupled_var_dir;
     334             : 
     335      900816 :     case Moose::PrimarySecondary:
     336      900816 :       return -_test_primary[_i][_qp] * jac_axial * tangent_component_in_coupled_var_dir *
     337      900816 :                  _phi_secondary[_j][_qp] +
     338      900816 :              _test_primary[_i][_qp] * this->_penalty * _phi_secondary[_j][_qp] *
     339      900816 :                  _secondary_node_tangent(_component) * tangent_component_in_coupled_var_dir;
     340             : 
     341     1250112 :     case Moose::PrimaryPrimary:
     342     1250112 :       return _test_primary[_i][_qp] * jac_axial * tangent_component_in_coupled_var_dir *
     343     1250112 :                  _phi_primary[_j][_qp] -
     344     1250112 :              _test_primary[_i][_qp] * this->_penalty * _phi_primary[_j][_qp] *
     345     1250112 :                  _secondary_node_tangent(_component) * tangent_component_in_coupled_var_dir;
     346             : 
     347           0 :     default:
     348           0 :       mooseError("Unsupported type");
     349             :       break;
     350             :   }
     351             :   return 0.0;
     352             : }
     353             : 
     354             : template <bool is_ad>
     355             : bool
     356           0 : RebarBondSlipConstraintTempl<is_ad>::getCoupledVarComponent(unsigned int var_num,
     357             :                                                             unsigned int & component)
     358             : {
     359           0 :   component = std::numeric_limits<unsigned int>::max();
     360             :   bool coupled_var_is_disp_var = false;
     361     4032990 :   for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
     362     4032990 :     if (var_num == _disp_vars_nums[i])
     363             :     {
     364             :       coupled_var_is_disp_var = true;
     365           0 :       component = i;
     366           0 :       break;
     367             :     }
     368             : 
     369           0 :   return coupled_var_is_disp_var;
     370             : }
     371             : 
     372             : template <bool is_ad>
     373             : void
     374      328032 : RebarBondSlipConstraintTempl<is_ad>::concreteRebarModel(const GenericReal<is_ad> slip,
     375             :                                                         const bondSlipData * const bond_slip,
     376             :                                                         GenericReal<is_ad> & bond_stress,
     377             :                                                         GenericReal<is_ad> & bond_stress_deriv,
     378             :                                                         Real & plastic_slip) const
     379             : {
     380      328032 :   const Real slip_min = bond_slip->slip_min;
     381      328032 :   const Real slip_max = bond_slip->slip_max;
     382      328032 :   GenericReal<is_ad> slip_ratio = std::abs(slip) / _transitional_slip;
     383             : 
     384      328032 :   const Real slope = 5.0 * _max_bond_stress / _transitional_slip;
     385      328032 :   const Real plastic_slip_max = slip_max - bond_slip->bond_stress_max_old / slope;
     386      328032 :   const Real plastic_slip_min = slip_min - bond_slip->bond_stress_min_old / slope;
     387             : 
     388      328032 :   plastic_slip = MetaPhysicL::raw_value(_bond_stress) / slope - MetaPhysicL::raw_value(slip);
     389             : 
     390      328032 :   bond_stress_deriv = 0.0;
     391      328032 :   bond_stress = 0.0;
     392      328032 :   if (slip >= slip_max || slip <= slip_min)
     393             :   {
     394      224014 :     if (std::abs(slip) < _transitional_slip)
     395             :     {
     396             :       // elastic load or unload
     397      410474 :       bond_stress = _max_bond_stress * MathUtils::sign(slip) *
     398      315492 :                     (5.0 * slip_ratio - 4.5 * Utility::pow<2>(slip_ratio) +
     399      315492 :                      1.4 * Utility::pow<3>(slip_ratio));
     400      220510 :       bond_stress_deriv =
     401      125528 :           _max_bond_stress *
     402      220510 :           (5.0 / _transitional_slip -
     403      220510 :            MathUtils::sign(slip) * 4.5 * 2.0 * slip_ratio / _transitional_slip +
     404      315492 :            MathUtils::sign(slip) * 1.4 * 3.0 * Utility::pow<2>(slip_ratio) / _transitional_slip);
     405             :     }
     406             :     else
     407             :     {
     408             :       // plastic slip
     409        5364 :       bond_stress = MathUtils::sign(slip) * _max_bond_stress * 1.9;
     410        3504 :       bond_stress_deriv = 0;
     411             :     }
     412             :   }
     413      104018 :   else if (slip > plastic_slip_max && slip < slip_max)
     414             :   {
     415             :     // unload after positive plastic slip
     416       50728 :     bond_stress =
     417       50728 :         (slip - plastic_slip_max) * bond_slip->bond_stress_max / (slip_max - plastic_slip_max);
     418       50728 :     bond_stress_deriv = slope;
     419             :   }
     420       53290 :   else if (slip < plastic_slip_min && slip > slip_min)
     421             :   {
     422             :     // unload after negative plastic slip
     423       24168 :     bond_stress =
     424       24168 :         (slip - plastic_slip_min) * bond_slip->bond_stress_min / (slip_min - plastic_slip_min);
     425       24168 :     bond_stress_deriv = slope;
     426             :   }
     427             :   else
     428             :   {
     429             :     // slip at zero stress
     430       19398 :     bond_stress = 0;
     431       19398 :     bond_stress_deriv = 0;
     432             :   }
     433      328032 : }
     434             : 
     435             : template <bool is_ad>
     436             : void
     437       68544 : RebarBondSlipConstraintTempl<is_ad>::elasticPerfectPlasticModel(
     438             :     const GenericReal<is_ad> slip,
     439             :     const bondSlipData * const bond_slip,
     440             :     GenericReal<is_ad> & bond_stress,
     441             :     GenericReal<is_ad> & bond_stress_deriv,
     442             :     Real & plastic_slip) const
     443             : {
     444       68544 :   const Real slip_min = bond_slip->slip_min;
     445       68544 :   const Real slip_max = bond_slip->slip_max;
     446             : 
     447       68544 :   const Real slope = _max_bond_stress / _transitional_slip;
     448       68544 :   const Real plastic_slip_max = slip_max - bond_slip->bond_stress_max_old / slope;
     449       68544 :   const Real plastic_slip_min = slip_min - bond_slip->bond_stress_min_old / slope;
     450             : 
     451       68544 :   plastic_slip = MetaPhysicL::raw_value(_bond_stress) / slope - MetaPhysicL::raw_value(slip);
     452             : 
     453       68544 :   bond_stress_deriv = 0.0;
     454       68544 :   bond_stress = 0.0;
     455       68544 :   if (slip >= slip_max || slip <= slip_min)
     456             :   {
     457       38996 :     if (std::abs(slip) < _transitional_slip)
     458             :     {
     459             :       // elastic load or unload
     460       36392 :       bond_stress = slip * slope;
     461       36392 :       bond_stress_deriv = slope;
     462             :     }
     463             :     else
     464             :     {
     465             :       // plastic slip
     466        2604 :       bond_stress = MathUtils::sign(slip) * _max_bond_stress;
     467        2604 :       bond_stress_deriv = 0;
     468             :     }
     469             :   }
     470       29548 :   else if (slip > plastic_slip_max && slip < bond_slip->slip_max)
     471             :   {
     472             :     // unload after positive plastic slip
     473       18912 :     bond_stress = (slip - plastic_slip_max) * bond_slip->bond_stress_max /
     474       18912 :                   (bond_slip->slip_max - plastic_slip_max);
     475       18912 :     bond_stress_deriv = slope;
     476             :   }
     477       10636 :   else if (slip < plastic_slip_min && slip > bond_slip->slip_min)
     478             :   {
     479             :     // unload after negative plastic slip
     480        8224 :     bond_stress = (slip - plastic_slip_min) * bond_slip->bond_stress_min /
     481        8224 :                   (bond_slip->slip_min - plastic_slip_min);
     482        8224 :     bond_stress_deriv = slope;
     483             :   }
     484             :   else
     485             :   {
     486             :     // slip at zero stress
     487           0 :     bond_stress = 0;
     488           0 :     bond_stress_deriv = 0;
     489             :   }
     490       68544 : }
     491             : 
     492             : template class RebarBondSlipConstraintTempl<false>;
     493             : template class RebarBondSlipConstraintTempl<true>;

Generated by: LCOV version 1.14