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 "GeneralizedPlaneStrainUserObjectOSPD.h"
      11             : #include "RankFourTensor.h"
      12             : #include "Function.h"
      13             : 
      14             : registerMooseObject("PeridynamicsApp", GeneralizedPlaneStrainUserObjectOSPD);
      15             : 
      16             : InputParameters
      17          80 : GeneralizedPlaneStrainUserObjectOSPD::validParams()
      18             : {
      19          80 :   InputParameters params = GeneralizedPlaneStrainUserObjectBasePD::validParams();
      20          80 :   params.addClassDescription("Class for calculating the scalar residual and diagonal Jacobian "
      21             :                              "entry of generalized plane strain in OSPD formulation");
      22             : 
      23         160 :   params.addCoupledVar("out_of_plane_stress_variable",
      24             :                        "Auxiliary variable name for out-of-plane stress in GPS simulation");
      25             : 
      26          80 :   return params;
      27           0 : }
      28             : 
      29          44 : GeneralizedPlaneStrainUserObjectOSPD::GeneralizedPlaneStrainUserObjectOSPD(
      30          44 :     const InputParameters & parameters)
      31             :   : GeneralizedPlaneStrainUserObjectBasePD(parameters),
      32          44 :     _out_of_plane_stress_var(getVar("out_of_plane_stress_variable", 0))
      33             : {
      34          44 : }
      35             : 
      36             : void
      37      185450 : GeneralizedPlaneStrainUserObjectOSPD::execute()
      38             : {
      39             :   // dof_id_type for node i and j
      40      185450 :   dof_id_type node_i = _current_elem->node_id(0);
      41             :   dof_id_type node_j = _current_elem->node_id(1);
      42             : 
      43             :   // coordinates for node i and j
      44      185450 :   Point coord_i = _pdmesh.getNodeCoord(node_i);
      45      185450 :   Point coord_j = _pdmesh.getNodeCoord(node_j);
      46             : 
      47             :   // nodal area for node i and j
      48      185450 :   Real nv_i = _pdmesh.getNodeVolume(node_i);
      49      185450 :   Real nv_j = _pdmesh.getNodeVolume(node_j);
      50             : 
      51             :   // number of neighbors for node i and j, used to avoid repeated accounting nodal stress in
      52             :   // element-wise loop
      53             : 
      54             :   // calculate number of active neighbors for node i and j
      55      185450 :   std::vector<unsigned int> active_neighbors(_nnodes, 0);
      56      556350 :   for (unsigned int nd = 0; nd < _nnodes; nd++)
      57             :   {
      58      370900 :     std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
      59     4740304 :     for (unsigned int nb = 0; nb < bonds.size(); ++nb)
      60     4369404 :       if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[nb])) > 0.5)
      61     4320084 :         active_neighbors[nd]++;
      62             : 
      63      370900 :     if (active_neighbors[nd] == 0) // avoid dividing by zero
      64        1800 :       active_neighbors[nd] = 1;
      65      370900 :   }
      66             : 
      67      185450 :   Real bond_status = _bond_status_var->getElementalValue(_current_elem);
      68             : 
      69             :   // residual
      70      185450 :   _residual += (_out_of_plane_stress_var->getNodalValue(*_current_elem->node_ptr(0)) -
      71      185450 :                 _pressure.value(_t, coord_i) * _factor) *
      72      185450 :                nv_i / active_neighbors[0] * bond_status;
      73      185450 :   _residual += (_out_of_plane_stress_var->getNodalValue(*_current_elem->node_ptr(1)) -
      74      185450 :                 _pressure.value(_t, coord_j) * _factor) *
      75      185450 :                nv_j / active_neighbors[1] * bond_status;
      76             : 
      77             :   // diagonal jacobian
      78      185450 :   _jacobian += (_Cijkl[0](2, 2, 2, 2) * nv_i / active_neighbors[0] +
      79      185450 :                 _Cijkl[0](2, 2, 2, 2) * nv_j / active_neighbors[1]) *
      80             :                bond_status;
      81      185450 : }