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

Generated by: LCOV version 1.14