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 "ForceStabilizedSmallStrainMechanicsNOSPD.h"
      11             : #include "PeridynamicsMesh.h"
      12             : 
      13             : registerMooseObject("PeridynamicsApp", ForceStabilizedSmallStrainMechanicsNOSPD);
      14             : 
      15             : InputParameters
      16          56 : ForceStabilizedSmallStrainMechanicsNOSPD::validParams()
      17             : {
      18          56 :   InputParameters params = MechanicsBaseNOSPD::validParams();
      19          56 :   params.addClassDescription(
      20             :       "Class for calculating the residual and Jacobian for the force-stabilized peridynamic "
      21             :       "correspondence model under small strain assumptions");
      22             : 
      23         112 :   params.addRequiredParam<unsigned int>(
      24             :       "component",
      25             :       "An integer corresponding to the variable this kernel acts on (0 for x, 1 for y, 2 for z)");
      26             : 
      27          56 :   return params;
      28           0 : }
      29             : 
      30          40 : ForceStabilizedSmallStrainMechanicsNOSPD::ForceStabilizedSmallStrainMechanicsNOSPD(
      31          40 :     const InputParameters & parameters)
      32             :   : MechanicsBaseNOSPD(parameters),
      33          40 :     _sf_coeff(getMaterialProperty<Real>("stabilization_force_coeff")),
      34         120 :     _component(getParam<unsigned int>("component"))
      35             : {
      36          40 : }
      37             : 
      38             : void
      39      103488 : ForceStabilizedSmallStrainMechanicsNOSPD::computeLocalResidual()
      40             : {
      41      103488 :   Real sforce = 0.5 * (_sf_coeff[0] + _sf_coeff[1]) *
      42      103488 :                 (_disp_var[_component]->getNodalValue(*_current_elem->node_ptr(1)) -
      43      103488 :                  _disp_var[_component]->getNodalValue(*_current_elem->node_ptr(0))) *
      44      103488 :                 _bond_status;
      45             : 
      46             :   // For small strain assumptions, stress measures, i.e., Cauchy stress (Sigma), the first
      47             :   // Piola-Kirchhoff stress (P), and the second Piola-Kirchhoff stress (S) are approximately the
      48             :   // same. Thus, the nodal force state tensors are calculated using the Cauchy stresses,
      49             :   // i.e., T = Sigma * inv(Shape) * xi * multi.
      50             :   // Cauchy stress is calculated as Sigma = C * E.
      51             : 
      52      310464 :   for (unsigned int nd = 0; nd < _nnodes; ++nd)
      53      620928 :     for (unsigned int i = 0; i < _nnodes; ++i)
      54      620928 :       _local_re(i) += (i == 0 ? -1 : 1) * _multi[nd] *
      55      413952 :                       (_stress[nd] * _shape2[nd].inverse()).row(_component) * _origin_vec *
      56      413952 :                       _bond_status;
      57             : 
      58             :   // add fictitious force for model stabilization
      59      103488 :   _local_re(0) += -sforce;
      60      103488 :   _local_re(1) += sforce;
      61      103488 : }
      62             : 
      63             : void
      64        2352 : ForceStabilizedSmallStrainMechanicsNOSPD::computeLocalJacobian()
      65             : {
      66             :   // excludes dTi/dUj and dTj/dUi contribution which was considered as nonlocal contribution
      67        7056 :   for (unsigned int i = 0; i < _nnodes; ++i)
      68       14112 :     for (unsigned int j = 0; j < _nnodes; ++j)
      69       14112 :       _local_ke(i, j) += (i == 0 ? -1 : 1) * _multi[j] *
      70       18816 :                          (computeDSDU(_component, j) * _shape2[j].inverse()).row(_component) *
      71        9408 :                          _origin_vec * _bond_status;
      72             : 
      73             :   // compute local stabilization force jacobian
      74        2352 :   Real diag = 0.5 * (_sf_coeff[0] + _sf_coeff[1]);
      75        7056 :   for (unsigned int i = 0; i < _nnodes; ++i)
      76       14112 :     for (unsigned int j = 0; j < _nnodes; ++j)
      77       14112 :       _local_ke(i, j) += (i == j ? 1 : -1) * diag * _bond_status;
      78        2352 : }
      79             : 
      80             : void
      81        2352 : ForceStabilizedSmallStrainMechanicsNOSPD::computeNonlocalJacobian()
      82             : {
      83             :   // includes dTi/dUj and dTj/dUi contributions
      84        7056 :   for (unsigned int nd = 0; nd < _nnodes; ++nd)
      85             :   {
      86             :     // calculation of jacobian contribution to current_node's neighbors
      87        4704 :     std::vector<dof_id_type> ivardofs(_nnodes);
      88        4704 :     ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
      89        4704 :     std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
      90        4704 :     std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
      91             : 
      92             :     Real vol_nb;
      93             :     RealGradient origin_vec_nb;
      94        4704 :     RankTwoTensor dFdUk, dSxdUkx;
      95             : 
      96       63552 :     for (unsigned int nb = 0; nb < neighbors.size(); ++nb)
      97       58848 :       if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[nb])) > 0.5)
      98             :       {
      99       56656 :         Node * neighbor_nb = _pdmesh.nodePtr(neighbors[nb]);
     100       56656 :         ivardofs[1] = neighbor_nb->dof_number(_sys.number(), _var.number(), 0);
     101       56656 :         vol_nb = _pdmesh.getNodeVolume(neighbors[nb]);
     102             : 
     103             :         // obtain bond ndnb's origin vector
     104       56656 :         origin_vec_nb = _pdmesh.getNodeCoord(neighbor_nb->id()) -
     105      113312 :                         _pdmesh.getNodeCoord(_current_elem->node_id(nd));
     106             : 
     107             :         dFdUk.zero();
     108      169968 :         for (unsigned int i = 0; i < _dim; ++i)
     109      113312 :           dFdUk(_component, i) =
     110      113312 :               vol_nb * _horizon_radius[nd] / origin_vec_nb.norm() * origin_vec_nb(i);
     111             : 
     112       56656 :         dFdUk *= _shape2[nd].inverse();
     113      113312 :         dSxdUkx = _Jacobian_mult[nd] * 0.5 * (dFdUk.transpose() + dFdUk);
     114             : 
     115             :         _local_ke.zero();
     116      169968 :         for (unsigned int i = 0; i < _nnodes; ++i)
     117      339936 :           for (unsigned int j = 0; j < _nnodes; ++j)
     118      339936 :             _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 0 : 1) * _multi[nd] *
     119      226624 :                               (dSxdUkx * _shape2[nd].inverse()).row(_component) * _origin_vec *
     120      226624 :                               _bond_status;
     121             : 
     122       56656 :         addJacobian(_assembly, _local_ke, _ivardofs, ivardofs, _var.scalingFactor());
     123             : 
     124       56656 :         if (_has_diag_save_in)
     125             :         {
     126           0 :           unsigned int rows = _nnodes;
     127           0 :           DenseVector<Real> diag(rows);
     128           0 :           for (unsigned int i = 0; i < rows; ++i)
     129           0 :             diag(i) = _local_ke(i, i);
     130             : 
     131             :           Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     132           0 :           for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
     133           0 :             _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
     134             :         }
     135             :       }
     136        4704 :   }
     137        2352 : }
     138             : 
     139             : void
     140        3528 : ForceStabilizedSmallStrainMechanicsNOSPD::computeLocalOffDiagJacobian(
     141             :     unsigned int jvar_num, unsigned int coupled_component)
     142             : {
     143             :   _local_ke.zero();
     144        3528 :   if (_temp_coupled && jvar_num == _temp_var->number())
     145             :   {
     146        1176 :     std::vector<RankTwoTensor> dSdT(_nnodes);
     147        3528 :     for (unsigned int nd = 0; nd < _nnodes; ++nd)
     148        4704 :       for (unsigned int es = 0; es < _deigenstrain_dT.size(); ++es)
     149        2352 :         dSdT[nd] = -(_Jacobian_mult[nd] * (*_deigenstrain_dT[es])[nd]);
     150             : 
     151        3528 :     for (unsigned int i = 0; i < _nnodes; ++i)
     152        7056 :       for (unsigned int j = 0; j < _nnodes; ++j)
     153        7056 :         _local_ke(i, j) += (i == 0 ? -1 : 1) * _multi[j] *
     154        4704 :                            (dSdT[j] * _shape2[j].inverse()).row(_component) * _origin_vec *
     155        4704 :                            _bond_status;
     156        1176 :   }
     157             :   else
     158             :   {
     159        7056 :     for (unsigned int i = 0; i < _nnodes; ++i)
     160       14112 :       for (unsigned int j = 0; j < _nnodes; ++j)
     161        9408 :         _local_ke(i, j) +=
     162       14112 :             (i == 0 ? -1 : 1) * _multi[j] *
     163       18816 :             (computeDSDU(coupled_component, j) * _shape2[j].inverse()).row(_component) *
     164        9408 :             _origin_vec * _bond_status;
     165             :   }
     166        3528 : }
     167             : 
     168             : void
     169        3528 : ForceStabilizedSmallStrainMechanicsNOSPD::computePDNonlocalOffDiagJacobian(
     170             :     unsigned int jvar_num, unsigned int coupled_component)
     171             : {
     172        3528 :   if (_temp_coupled && jvar_num == _temp_var->number())
     173             :   {
     174             :     // no nonlocal contribution from temperature
     175             :   }
     176             :   else
     177             :   {
     178        7056 :     for (unsigned int nd = 0; nd < _nnodes; ++nd)
     179             :     {
     180             :       // calculation of jacobian contribution to current_node's neighbors
     181        4704 :       std::vector<dof_id_type> jvardofs(_nnodes);
     182        4704 :       jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
     183        4704 :       std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
     184        4704 :       std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
     185             : 
     186             :       Real vol_nb;
     187             :       RealGradient origin_vec_nb;
     188        4704 :       RankTwoTensor dFdUk, dSxdUky;
     189             : 
     190       63552 :       for (unsigned int nb = 0; nb < neighbors.size(); ++nb)
     191       58848 :         if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[nb])) > 0.5)
     192             : 
     193             :         {
     194       56656 :           Node * neighbor_nb = _pdmesh.nodePtr(neighbors[nb]);
     195       56656 :           jvardofs[1] = neighbor_nb->dof_number(_sys.number(), jvar_num, 0);
     196       56656 :           vol_nb = _pdmesh.getNodeVolume(neighbors[nb]);
     197             : 
     198             :           // obtain bond ndnb's origin vector
     199       56656 :           origin_vec_nb = _pdmesh.getNodeCoord(neighbor_nb->id()) -
     200      113312 :                           _pdmesh.getNodeCoord(_current_elem->node_id(nd));
     201             : 
     202             :           dFdUk.zero();
     203      169968 :           for (unsigned int i = 0; i < _dim; ++i)
     204      113312 :             dFdUk(coupled_component, i) =
     205      113312 :                 _horizon_radius[nd] / origin_vec_nb.norm() * origin_vec_nb(i) * vol_nb;
     206             : 
     207       56656 :           dFdUk *= _shape2[nd].inverse();
     208      113312 :           dSxdUky = _Jacobian_mult[nd] * 0.5 * (dFdUk.transpose() + dFdUk);
     209             : 
     210             :           _local_ke.zero();
     211      169968 :           for (unsigned int i = 0; i < _nnodes; ++i)
     212      339936 :             for (unsigned int j = 0; j < _nnodes; ++j)
     213      339936 :               _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 0 : 1) * _multi[nd] *
     214      226624 :                                 (dSxdUky * _shape2[nd].inverse()).row(_component) * _origin_vec *
     215      226624 :                                 _bond_status;
     216             : 
     217       56656 :           addJacobian(_assembly, _local_ke, _ivardofs, jvardofs, _var.scalingFactor());
     218             :         }
     219        4704 :     }
     220             :   }
     221        3528 : }