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 "WeakPlaneStressNOSPD.h"
      11             : #include "MooseVariable.h"
      12             : #include "PeridynamicsMesh.h"
      13             : 
      14             : registerMooseObject("PeridynamicsApp", WeakPlaneStressNOSPD);
      15             : 
      16             : InputParameters
      17          29 : WeakPlaneStressNOSPD::validParams()
      18             : {
      19          29 :   InputParameters params = MechanicsBaseNOSPD::validParams();
      20          29 :   params.addClassDescription(
      21             :       "Class for calculating the residual and the Jacobian for the peridynamic plane "
      22             :       "stress model using weak formulation based on peridynamic correspondence models");
      23             : 
      24          29 :   return params;
      25           0 : }
      26             : 
      27          16 : WeakPlaneStressNOSPD::WeakPlaneStressNOSPD(const InputParameters & parameters)
      28          16 :   : MechanicsBaseNOSPD(parameters)
      29             : {
      30          16 : }
      31             : 
      32             : void
      33      162680 : WeakPlaneStressNOSPD::computeLocalResidual()
      34             : {
      35      162680 :   _local_re(0) = _stress[0](2, 2) * _weights[0] * _node_vol[0] * _bond_status;
      36      162680 :   _local_re(1) = _stress[1](2, 2) * _weights[1] * _node_vol[1] * _bond_status;
      37      162680 : }
      38             : 
      39             : void
      40        6848 : WeakPlaneStressNOSPD::computeLocalJacobian()
      41             : {
      42       20544 :   for (unsigned int i = 0; i < _nnodes; ++i)
      43       41088 :     for (unsigned int j = 0; j < _nnodes; ++j)
      44       41088 :       _local_ke(i, j) += (i == j ? 1 : 0) * _Jacobian_mult[i](2, 2, 2, 2) * _weights[i] *
      45       27392 :                          _node_vol[i] * _bond_status;
      46        6848 : }
      47             : 
      48             : void
      49       14284 : WeakPlaneStressNOSPD::computeLocalOffDiagJacobian(unsigned int jvar_num,
      50             :                                                   unsigned int coupled_component)
      51             : {
      52             :   _local_ke.zero();
      53       14284 :   if (_temp_coupled && jvar_num == _temp_var->number()) // for coupled temperature variable
      54             :   {
      55         588 :     std::vector<RankTwoTensor> dSdT(_nnodes);
      56        1764 :     for (unsigned int nd = 0; nd < _nnodes; nd++)
      57        2352 :       for (unsigned int es = 0; es < _deigenstrain_dT.size(); es++)
      58        1176 :         dSdT[nd] = -_Jacobian_mult[nd] * (*_deigenstrain_dT[es])[nd];
      59             : 
      60        1764 :     for (unsigned int i = 0; i < _nnodes; ++i)
      61        3528 :       for (unsigned int j = 0; j < _nnodes; ++j)
      62        2352 :         _local_ke(i, j) +=
      63        3528 :             (i == j ? 1 : 0) * dSdT[i](2, 2) * _weights[i] * _node_vol[i] * _bond_status;
      64         588 :   }
      65             :   else // for coupled displacement variables
      66             :   {
      67             :     // dSidUi and dSjdUj are considered as local off-diagonal Jacobian
      68             :     // contributions to their neighbors are considered as nonlocal off-diagonal
      69       41088 :     for (unsigned int i = 0; i < _nnodes; ++i)
      70       82176 :       for (unsigned int j = 0; j < _nnodes; ++j)
      71       82176 :         _local_ke(i, j) += (i == j ? 1 : 0) * computeDSDU(coupled_component, j)(2, 2) *
      72       54784 :                            _weights[j] * _node_vol[j] * _bond_status;
      73             :   }
      74       14284 : }
      75             : 
      76             : void
      77        1764 : WeakPlaneStressNOSPD::computePDNonlocalOffDiagJacobian(unsigned int jvar_num,
      78             :                                                        unsigned int coupled_component)
      79             : {
      80        1764 :   if (_temp_coupled && jvar_num == _temp_var->number())
      81             :   {
      82             :     // no nonlocal contribution from temperature
      83             :   }
      84             :   else
      85             :   {
      86        3528 :     for (unsigned int nd = 0; nd < _nnodes; nd++)
      87             :     {
      88             :       // calculation of jacobian contribution to current_node's neighbors
      89             :       // NOT including the contribution to nodes i and j, which is considered as local
      90             :       // off-diagonal
      91        2352 :       std::vector<dof_id_type> jvardofs(_nnodes);
      92        2352 :       jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
      93        2352 :       std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
      94        2352 :       std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
      95             : 
      96             :       dof_id_type nb_index =
      97        2352 :           std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
      98        2352 :           neighbors.begin();
      99             :       std::vector<dof_id_type> dg_neighbors =
     100        2352 :           _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
     101             : 
     102             :       Real vol_nb, weight_nb;
     103             :       RealGradient origin_vec_nb;
     104        2352 :       RankTwoTensor dFdUk, dPdUk;
     105             : 
     106       31776 :       for (unsigned int nb = 0; nb < dg_neighbors.size(); nb++)
     107       29424 :         if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
     108             :         {
     109       28328 :           jvardofs[1] =
     110       28328 :               _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])->dof_number(_sys.number(), jvar_num, 0);
     111       28328 :           vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
     112             : 
     113       28328 :           origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
     114       56656 :                           _pdmesh.getNodeCoord(_current_elem->node_id(nd));
     115             : 
     116       28328 :           weight_nb = _horizon_radius[nd] / origin_vec_nb.norm();
     117             : 
     118             :           dFdUk.zero();
     119       84984 :           for (unsigned int i = 0; i < _dim; ++i)
     120       56656 :             dFdUk(coupled_component, i) = weight_nb * origin_vec_nb(i) * vol_nb;
     121             : 
     122       28328 :           dFdUk *= _shape2[nd].inverse();
     123             : 
     124       28328 :           dPdUk = _Jacobian_mult[nd] * 0.5 *
     125       56656 :                   (dFdUk.transpose() * _dgrad[nd] + _dgrad[nd].transpose() * dFdUk);
     126             : 
     127             :           _local_ke.zero();
     128       28328 :           _local_ke(nd, 1) = dPdUk(2, 2) * _weights[nd] * _node_vol[nd] * _bond_status;
     129             : 
     130       28328 :           addJacobian(_assembly, _local_ke, _ivardofs, jvardofs, _var.scalingFactor());
     131             :         }
     132        2352 :     }
     133             :   }
     134        1764 : }