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 "ComputeStrainBaseNOSPD.h"
      11             : #include "ElasticityTensorTools.h"
      12             : 
      13             : #include "libmesh/quadrature.h"
      14             : 
      15             : InputParameters
      16         564 : ComputeStrainBaseNOSPD::validParams()
      17             : {
      18         564 :   InputParameters params = MechanicsMaterialBasePD::validParams();
      19         564 :   params.addClassDescription(
      20             :       "Base material strain class for the peridynamic correspondence models");
      21             : 
      22        1128 :   MooseEnum stabilization_option("FORCE BOND_HORIZON_I BOND_HORIZON_II");
      23        1128 :   params.addRequiredParam<MooseEnum>(
      24             :       "stabilization",
      25             :       stabilization_option,
      26             :       "Stabilization techniques used for the peridynamic correspondence model");
      27        1128 :   params.addParam<bool>("plane_strain",
      28        1128 :                         false,
      29             :                         "Plane strain problem or not, this will affect the mechanical stretch "
      30             :                         "calculation for problem with eigenstrains");
      31        1128 :   params.addParam<std::vector<MaterialPropertyName>>(
      32             :       "eigenstrain_names", {}, "List of eigenstrains to be applied in this strain calculation");
      33             : 
      34         564 :   return params;
      35         564 : }
      36             : 
      37         444 : ComputeStrainBaseNOSPD::ComputeStrainBaseNOSPD(const InputParameters & parameters)
      38             :   : DerivativeMaterialInterface<MechanicsMaterialBasePD>(parameters),
      39         444 :     _stabilization(getParam<MooseEnum>("stabilization")),
      40         888 :     _plane_strain(getParam<bool>("plane_strain")),
      41         888 :     _Cijkl(getMaterialProperty<RankFourTensor>("elasticity_tensor")),
      42         888 :     _eigenstrain_names(getParam<std::vector<MaterialPropertyName>>("eigenstrain_names")),
      43         444 :     _eigenstrains(_eigenstrain_names.size()),
      44         444 :     _shape2(declareProperty<RankTwoTensor>("rank_two_shape_tensor")),
      45         444 :     _deformation_gradient(declareProperty<RankTwoTensor>("deformation_gradient")),
      46         444 :     _ddgraddu(declareProperty<RankTwoTensor>("ddeformation_gradient_du")),
      47         444 :     _ddgraddv(declareProperty<RankTwoTensor>("ddeformation_gradient_dv")),
      48         444 :     _ddgraddw(declareProperty<RankTwoTensor>("ddeformation_gradient_dw")),
      49         444 :     _total_strain(declareProperty<RankTwoTensor>("total_strain")),
      50         444 :     _mechanical_strain(declareProperty<RankTwoTensor>("mechanical_strain")),
      51         444 :     _multi(declareProperty<Real>("multi")),
      52         888 :     _sf_coeff(declareProperty<Real>("stabilization_force_coeff"))
      53             : {
      54         615 :   for (unsigned int i = 0; i < _eigenstrains.size(); ++i)
      55             :   {
      56             :     _eigenstrain_names[i] = _eigenstrain_names[i];
      57         171 :     _eigenstrains[i] = &getMaterialProperty<RankTwoTensor>(_eigenstrain_names[i]);
      58             :   }
      59         444 : }
      60             : 
      61             : void
      62        4896 : ComputeStrainBaseNOSPD::initQpStatefulProperties()
      63             : {
      64        4896 :   _mechanical_strain[_qp].zero();
      65        4896 :   _total_strain[_qp].zero();
      66        4896 :   _deformation_gradient[_qp].setToIdentity();
      67             : 
      68        4896 :   if (_qrule->n_points() < 2) // Gauss_Lobatto: order = 2n-3 (n is number of qp)
      69           0 :     mooseError(
      70             :         "NOSPD models require Gauss_Lobatto rule and a minimum of 2 quadrature points, i.e., "
      71             :         "order of FIRST");
      72        4896 : }
      73             : 
      74             : void
      75     2860812 : ComputeStrainBaseNOSPD::computeQpDeformationGradient()
      76             : {
      77     2860812 :   _shape2[_qp].zero();
      78     2860812 :   _deformation_gradient[_qp].zero();
      79     2860812 :   _ddgraddu[_qp].zero();
      80     2860812 :   _ddgraddv[_qp].zero();
      81     2860812 :   _ddgraddw[_qp].zero();
      82     2860812 :   _multi[_qp] = 0.0;
      83     2860812 :   _sf_coeff[_qp] = 0.0;
      84             : 
      85     2860812 :   if (_dim == 2)
      86     2574432 :     _shape2[_qp](2, 2) = _deformation_gradient[_qp](2, 2) = 1.0;
      87             : 
      88     2860812 :   if (_bond_status_var->getElementalValue(_current_elem) > 0.5)
      89             :   {
      90     2834100 :     if (_stabilization == "FORCE")
      91             :     {
      92      101520 :       computeConventionalQpDeformationGradient();
      93             : 
      94      101520 :       _sf_coeff[_qp] = ElasticityTensorTools::getIsotropicYoungsModulus(_Cijkl[_qp]) *
      95      101520 :                        _pdmesh.getNodeAverageSpacing(_current_elem->node_id(_qp)) *
      96      101520 :                        _horizon_radius[_qp] / _origin_vec.norm();
      97      101520 :       _multi[_qp] = _horizon_radius[_qp] / _origin_vec.norm() * _node_vol[0] * _node_vol[1];
      98             :     }
      99     2732580 :     else if (_stabilization == "BOND_HORIZON_I" || _stabilization == "BOND_HORIZON_II")
     100     2732580 :       computeBondHorizonQpDeformationGradient();
     101             :     else
     102           0 :       paramError("stabilization",
     103             :                  "Unknown stabilization scheme for peridynamic correspondence model!");
     104             :   }
     105             :   else
     106             :   {
     107       26712 :     _shape2[_qp].setToIdentity();
     108       26712 :     _deformation_gradient[_qp].setToIdentity();
     109             :   }
     110     2860812 : }
     111             : 
     112             : void
     113      101520 : ComputeStrainBaseNOSPD::computeConventionalQpDeformationGradient()
     114             : {
     115      101520 :   std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(_qp));
     116      101520 :   std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(_qp));
     117             : 
     118             :   // calculate the shape tensor and prepare the deformation gradient tensor
     119             :   Real vol_nb, weight_nb;
     120             :   RealGradient origin_vec_nb, current_vec_nb;
     121             : 
     122     1376280 :   for (unsigned int nb = 0; nb < neighbors.size(); ++nb)
     123     1274760 :     if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[nb])) > 0.5)
     124             :     {
     125     1247040 :       vol_nb = _pdmesh.getNodeVolume(neighbors[nb]);
     126     1247040 :       origin_vec_nb =
     127     1247040 :           _pdmesh.getNodeCoord(neighbors[nb]) - _pdmesh.getNodeCoord(_current_elem->node_id(_qp));
     128             : 
     129     3741120 :       for (unsigned int k = 0; k < _dim; ++k)
     130     2494080 :         current_vec_nb(k) = origin_vec_nb(k) +
     131     4988160 :                             _disp_var[k]->getNodalValue(*_pdmesh.nodePtr(neighbors[nb])) -
     132     2494080 :                             _disp_var[k]->getNodalValue(*_current_elem->node_ptr(_qp));
     133             : 
     134     1247040 :       weight_nb = _horizon_radius[_qp] / origin_vec_nb.norm();
     135     3741120 :       for (unsigned int k = 0; k < _dim; ++k)
     136             :       {
     137     7482240 :         for (unsigned int l = 0; l < _dim; ++l)
     138             :         {
     139     4988160 :           _shape2[_qp](k, l) += weight_nb * origin_vec_nb(k) * origin_vec_nb(l) * vol_nb;
     140     4988160 :           _deformation_gradient[_qp](k, l) +=
     141     4988160 :               weight_nb * current_vec_nb(k) * origin_vec_nb(l) * vol_nb;
     142             :         }
     143             :         // calculate derivatives of deformation_gradient w.r.t displacements of node nb
     144     2494080 :         _ddgraddu[_qp](0, k) += -weight_nb * origin_vec_nb(k) * vol_nb;
     145     2494080 :         _ddgraddv[_qp](1, k) += -weight_nb * origin_vec_nb(k) * vol_nb;
     146     2494080 :         if (_dim == 3)
     147           0 :           _ddgraddw[_qp](2, k) += -weight_nb * origin_vec_nb(k) * vol_nb;
     148             :       }
     149             :     }
     150             : 
     151             :   // finalize the deformation gradient tensor
     152      101520 :   if (_shape2[_qp].det() == 0.)
     153           0 :     mooseError("Singular shape tensor is detected in ComputeStrainBaseNOSPD! Use "
     154             :                "SingularShapeTensorEliminatorUserObjectPD to avoid singular shape tensor!");
     155             : 
     156      101520 :   _deformation_gradient[_qp] *= _shape2[_qp].inverse();
     157      101520 :   _ddgraddu[_qp] *= _shape2[_qp].inverse();
     158      101520 :   _ddgraddv[_qp] *= _shape2[_qp].inverse();
     159      101520 :   _ddgraddw[_qp] *= _shape2[_qp].inverse();
     160      101520 : }
     161             : 
     162             : void
     163     2732580 : ComputeStrainBaseNOSPD::computeBondHorizonQpDeformationGradient()
     164             : {
     165     2732580 :   std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(_qp));
     166     2732580 :   std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(_qp));
     167             : 
     168             :   dof_id_type nb_index =
     169     2732580 :       std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - _qp)) -
     170     2732580 :       neighbors.begin();
     171             :   std::vector<dof_id_type> dg_neighbors =
     172     2732580 :       _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(_qp), nb_index);
     173     2732580 :   Real dg_sub_vol = _pdmesh.getHorizonSubsetVolume(_current_elem->node_id(_qp), nb_index);
     174     2732580 :   Real dg_sub_vol_sum = _pdmesh.getHorizonSubsetVolumeSum(_current_elem->node_id(_qp));
     175             : 
     176             :   // calculate the shape tensor and prepare the deformation gradient tensor
     177             :   Real vol_nb, weight_nb;
     178             :   RealGradient origin_vec_nb, current_vec_nb;
     179             : 
     180    48914754 :   for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
     181    46182174 :     if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
     182             :     {
     183    46035080 :       vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
     184    46035080 :       origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
     185    46035080 :                       _pdmesh.getNodeCoord(_current_elem->node_id(_qp));
     186             : 
     187   142210536 :       for (unsigned int k = 0; k < _dim; ++k)
     188    96175456 :         current_vec_nb(k) =
     189    96175456 :             origin_vec_nb(k) +
     190   192350912 :             _disp_var[k]->getNodalValue(*_pdmesh.nodePtr(neighbors[dg_neighbors[nb]])) -
     191    96175456 :             _disp_var[k]->getNodalValue(*_current_elem->node_ptr(_qp));
     192             : 
     193    46035080 :       weight_nb = _horizon_radius[_qp] / origin_vec_nb.norm();
     194   142210536 :       for (unsigned int k = 0; k < _dim; ++k)
     195             :       {
     196   300842256 :         for (unsigned int l = 0; l < _dim; ++l)
     197             :         {
     198   204666800 :           _shape2[_qp](k, l) += weight_nb * origin_vec_nb(k) * origin_vec_nb(l) * vol_nb;
     199   204666800 :           _deformation_gradient[_qp](k, l) +=
     200   204666800 :               weight_nb * current_vec_nb(k) * origin_vec_nb(l) * vol_nb;
     201             :         }
     202             :         // calculate derivatives of deformation_gradient w.r.t displacements of node i
     203    96175456 :         _ddgraddu[_qp](0, k) += -weight_nb * origin_vec_nb(k) * vol_nb;
     204    96175456 :         _ddgraddv[_qp](1, k) += -weight_nb * origin_vec_nb(k) * vol_nb;
     205    96175456 :         if (_dim == 3)
     206    12315888 :           _ddgraddw[_qp](2, k) += -weight_nb * origin_vec_nb(k) * vol_nb;
     207             :       }
     208             :     }
     209             :   // finalize the deformation gradient and its derivatives
     210     2732580 :   if (_shape2[_qp].det() == 0.)
     211           0 :     computeConventionalQpDeformationGradient(); // this will affect the corresponding jacobian
     212             :                                                 // calculation
     213             :   else
     214             :   {
     215     2732580 :     _deformation_gradient[_qp] *= _shape2[_qp].inverse();
     216     2732580 :     _ddgraddu[_qp] *= _shape2[_qp].inverse();
     217     2732580 :     _ddgraddv[_qp] *= _shape2[_qp].inverse();
     218     2732580 :     _ddgraddw[_qp] *= _shape2[_qp].inverse();
     219             :   }
     220             : 
     221             :   // force state multiplier
     222     2732580 :   if (_stabilization == "BOND_HORIZON_I")
     223     2584200 :     _multi[_qp] = _horizon_radius[_qp] / _origin_vec.norm() * _node_vol[0] * _node_vol[1] *
     224     2584200 :                   dg_sub_vol / _horizon_vol[_qp];
     225      148380 :   else if (_stabilization == "BOND_HORIZON_II")
     226      148380 :     _multi[_qp] = _node_vol[_qp] * dg_sub_vol / dg_sub_vol_sum;
     227     2732580 : }
     228             : 
     229             : void
     230     1430406 : ComputeStrainBaseNOSPD::computeProperties()
     231             : {
     232     1430406 :   setupMeshRelatedData();     // function from base class
     233     1430406 :   computeBondCurrentLength(); // current length of a bond from base class
     234     1430406 :   computeBondStretch();       // stretch of a bond
     235             : 
     236     4291218 :   for (_qp = 0; _qp < _nnodes; ++_qp)
     237     2860812 :     computeQpStrain();
     238     1430406 : }
     239             : 
     240             : void
     241     1430406 : ComputeStrainBaseNOSPD::computeBondStretch()
     242             : {
     243     1430406 :   _total_stretch[0] = _current_len / _origin_vec.norm() - 1.0;
     244     1430406 :   _total_stretch[1] = _total_stretch[0];
     245     1430406 :   _mechanical_stretch[0] = _total_stretch[0];
     246             : 
     247             :   Real factor = 1.0;
     248     1430406 :   if (_plane_strain)
     249      144738 :     factor = 1.0 + ElasticityTensorTools::getIsotropicPoissonsRatio(_Cijkl[0]);
     250             : 
     251     1793096 :   for (auto es : _eigenstrains)
     252      362690 :     _mechanical_stretch[0] -= 0.5 * factor * ((*es)[0](2, 2) + (*es)[1](2, 2));
     253             : 
     254     1430406 :   _mechanical_stretch[1] = _mechanical_stretch[0];
     255     1430406 : }