Kernel class for bond-associated correspondence material model for finite strain. More...

#include <FiniteStrainMechanicsNOSPD.h>

Inheritance diagram for FiniteStrainMechanicsNOSPD:

Public Member Functions
	FiniteStrainMechanicsNOSPD (const InputParameters &parameters)

virtual void	computeOffDiagJacobian (MooseVariableFEBase &jvar) override

virtual void	initialSetup () override

virtual void	prepare () override

Protected Member Functions
virtual void	computeLocalResidual () override

virtual void	computeLocalJacobian () override

virtual void	computeNonlocalJacobian () override

virtual void	computeLocalOffDiagJacobian (unsigned int coupled_component) override
	Function to compute local contribution to the off-diagonal Jacobian at the current nodes. More...

virtual void	computePDNonlocalOffDiagJacobian (unsigned int jvar_num, unsigned int coupled_component) override
	Function to compute nonlocal contribution to the off-diagonal Jacobian at the current nodes. More...

virtual RankTwoTensor	computeDSDU (unsigned int component, unsigned int nd) override
	Function to compute derivative of stress with respect to displacements. More...

RankFourTensor	computeDSDFhat (unsigned int nd)
	Function to compute derivative of stress with respect to derived deformation gradient. More...

Real	computeDJDU (unsigned int component, unsigned int nd)
	Function to compute derivative of determinant of deformation gradient with respect to displacements. More...

RankTwoTensor	computeDinvFTDU (unsigned int component, unsigned int nd)
	Function to compute derivative of deformation gradient inverse with respect to displacements. More...

Protected Attributes
const unsigned int	_component
	The index of displacement component. More...

std::vector< MooseVariable * >	_disp_var
	displacement variables More...

unsigned int	_ndisp
	number of displacement components More...

const std::vector< RealGradient > *	_orientation
	Vector of bond in current configuration. More...

std::vector< dof_id_type >	_ivardofs_ij
	Current variable dof numbers for nodes i and j. More...

RealGradient	_cur_ori_ij
	Vector of bond in current configuration. More...

Real	_cur_len_ij
	Current bond length. More...


const MaterialProperty< RankTwoTensor > &	_dgrad_old
	Material point based material property. More...

const MaterialProperty< RankTwoTensor > &	_E_inc

const MaterialProperty< RankTwoTensor > &	_R_inc


const MaterialProperty< Real > &	_multi
	Material point based material properties. More...

const MaterialProperty< RankTwoTensor > &	_stress

const MaterialProperty< RankTwoTensor > &	_shape2

const MaterialProperty< RankTwoTensor > &	_dgrad

const MaterialProperty< RankTwoTensor > &	_ddgraddu

const MaterialProperty< RankTwoTensor > &	_ddgraddv

const MaterialProperty< RankTwoTensor > &	_ddgraddw

const MaterialProperty< RankFourTensor > &	_Jacobian_mult

const std::vector< MaterialPropertyName >	_eigenstrain_names

std::vector< const MaterialProperty< RankTwoTensor > * >	_deigenstrain_dT


const bool	_temp_coupled
	Temperature variable. More...

MooseVariable *	_temp_var


const bool	_out_of_plane_strain_coupled
	Parameters for out-of-plane strain in weak plane stress formulation. More...

MooseVariable *	_out_of_plane_strain_var

Detailed Description

Kernel class for bond-associated correspondence material model for finite strain.

Definition at line 22 of file FiniteStrainMechanicsNOSPD.h.

Constructor & Destructor Documentation

◆ FiniteStrainMechanicsNOSPD()

FiniteStrainMechanicsNOSPD::FiniteStrainMechanicsNOSPD ( const InputParameters & parameters )

Definition at line 32 of file FiniteStrainMechanicsNOSPD.C.

   : MechanicsBaseNOSPD(parameters),
     _dgrad_old(getMaterialPropertyOld<RankTwoTensor>("deformation_gradient")),
     _E_inc(getMaterialProperty<RankTwoTensor>("strain_increment")),
     _R_inc(getMaterialProperty<RankTwoTensor>("rotation_increment")),
     _component(getParam<unsigned int>("component"))
 {
 }

Member Function Documentation

◆ computeDinvFTDU()

RankTwoTensor FiniteStrainMechanicsNOSPD::computeDinvFTDU	(	unsigned int	component,
		unsigned int	nd
	)

protected

Function to compute derivative of deformation gradient inverse with respect to displacements.

Parameters

component	The index of displacement component
nd	The local index of element node (either 1 or 2 for Edge2 element)

Returns: The calculated derivative

Definition at line 380 of file FiniteStrainMechanicsNOSPD.C.

 {
   // for finite formulation, compute the derivative of transpose of inverse of deformation gradient
   // w.r.t the solution components
   // d(inv(F)_Ji)/du = d(inv(F)_Ji)/dF_kL * dF_kL/du = - inv(F)_Jk * inv(F)_Li * dF_kL/du
   // the bases are gi, GJ, gk, GL, indictates the inverse transpose rather than the inverse
  
   // RankTwoTensor dinvFTdU;
   // dinvFTdU.zero();
   // RankTwoTensor invF = _dgrad[nd].inverse();
   // for (unsigned int i = 0; i < 3; ++i)
   //   for (unsigned int J = 0; J < 3; ++J)
   //     for (unsigned int k = 0; k < 3; ++k)
   //       for (unsigned int L = 0; L < 3; ++L)
   //       {
   //         if (component == 0)
   //           dinvFTdU(i, J) += -invF(J, k) * invF(L, i) * _ddgraddu[nd](k, L);
   //         else if (component == 1)
   //           dinvFTdU(i, J) += -invF(J, k) * invF(L, i) * _ddgraddv[nd](k, L);
   //         else if (component == 2)
   //           dinvFTdU(i, J) += -invF(J, k) * invF(L, i) * _ddgraddw[nd](k, L);
   //       }
  
   RankTwoTensor dinvFTdU;
   dinvFTdU.zero();
   RankTwoTensor invF = _dgrad[nd].inverse();
   if (component == 0)
   {
     dinvFTdU(0, 1) =
         _ddgraddu[nd](0, 2) * _dgrad[nd](2, 1) - _ddgraddu[nd](0, 1) * _dgrad[nd](2, 2);
     dinvFTdU(0, 2) =
         _ddgraddu[nd](0, 1) * _dgrad[nd](1, 2) - _ddgraddu[nd](0, 2) * _dgrad[nd](1, 1);
     dinvFTdU(1, 1) =
         _ddgraddu[nd](0, 0) * _dgrad[nd](2, 2) - _ddgraddu[nd](0, 2) * _dgrad[nd](2, 0);
     dinvFTdU(1, 2) =
         _ddgraddu[nd](0, 2) * _dgrad[nd](1, 0) - _ddgraddu[nd](0, 0) * _dgrad[nd](1, 2);
     dinvFTdU(2, 1) =
         _ddgraddu[nd](0, 1) * _dgrad[nd](2, 0) - _ddgraddu[nd](0, 0) * _dgrad[nd](2, 1);
     dinvFTdU(2, 2) =
         _ddgraddu[nd](0, 0) * _dgrad[nd](1, 1) - _ddgraddu[nd](0, 1) * _dgrad[nd](1, 0);
   }
   else if (component == 1)
   {
     dinvFTdU(0, 0) =
         _ddgraddv[nd](1, 1) * _dgrad[nd](2, 2) - _ddgraddv[nd](1, 2) * _dgrad[nd](2, 1);
     dinvFTdU(0, 2) =
         _ddgraddv[nd](1, 2) * _dgrad[nd](0, 1) - _ddgraddv[nd](0, 2) * _dgrad[nd](1, 1);
     dinvFTdU(1, 0) =
         _ddgraddv[nd](1, 2) * _dgrad[nd](2, 0) - _ddgraddv[nd](1, 0) * _dgrad[nd](2, 2);
     dinvFTdU(1, 2) =
         _ddgraddv[nd](1, 0) * _dgrad[nd](0, 2) - _ddgraddv[nd](1, 2) * _dgrad[nd](0, 0);
     dinvFTdU(2, 0) =
         _ddgraddv[nd](1, 0) * _dgrad[nd](2, 1) - _ddgraddv[nd](1, 1) * _dgrad[nd](2, 0);
     dinvFTdU(2, 2) =
         _ddgraddv[nd](1, 1) * _dgrad[nd](0, 0) - _ddgraddv[nd](1, 0) * _dgrad[nd](0, 1);
   }
   else if (component == 2)
   {
     dinvFTdU(0, 0) =
         _ddgraddw[nd](2, 2) * _dgrad[nd](1, 1) - _ddgraddw[nd](2, 1) * _dgrad[nd](1, 2);
     dinvFTdU(0, 1) =
         _ddgraddw[nd](2, 1) * _dgrad[nd](0, 2) - _ddgraddw[nd](2, 2) * _dgrad[nd](0, 1);
     dinvFTdU(1, 0) =
         _ddgraddw[nd](2, 0) * _dgrad[nd](1, 2) - _ddgraddw[nd](2, 2) * _dgrad[nd](1, 0);
     dinvFTdU(1, 1) =
         _ddgraddw[nd](2, 2) * _dgrad[nd](0, 0) - _ddgraddw[nd](2, 0) * _dgrad[nd](0, 2);
     dinvFTdU(2, 0) =
         _ddgraddw[nd](2, 1) * _dgrad[nd](1, 0) - _ddgraddw[nd](2, 0) * _dgrad[nd](1, 1);
     dinvFTdU(2, 1) =
         _ddgraddw[nd](2, 0) * _dgrad[nd](0, 1) - _ddgraddw[nd](2, 1) * _dgrad[nd](0, 0);
   }
   dinvFTdU /= _dgrad[nd].det();
   for (unsigned int i = 0; i < 3; ++i)
     for (unsigned int J = 0; J < 3; ++J)
       for (unsigned int k = 0; k < 3; ++k)
         for (unsigned int L = 0; L < 3; ++L)
         {
           if (component == 0)
             dinvFTdU(i, J) -= invF(i, J) * invF(L, k) * _ddgraddu[nd](k, L);
           else if (component == 1)
             dinvFTdU(i, J) -= invF(i, J) * invF(L, k) * _ddgraddv[nd](k, L);
           else if (component == 2)
             dinvFTdU(i, J) -= invF(i, J) * invF(L, k) * _ddgraddw[nd](k, L);
         }
  
   return dinvFTdU;
 }

Referenced by computeLocalJacobian(), and computeLocalOffDiagJacobian().

◆ computeDJDU()

Real FiniteStrainMechanicsNOSPD::computeDJDU	(	unsigned int	component,
		unsigned int	nd
	)

protected

Function to compute derivative of determinant of deformation gradient with respect to displacements.

Parameters

component	The index of displacement component
nd	The local index of element node (either 1 or 2 for Edge2 element)

Returns: The calculated derivative

Definition at line 356 of file FiniteStrainMechanicsNOSPD.C.

 {
   // for finite formulation, compute the derivative of determinant of deformation gradient w.r.t the
   // solution components
   // dJ / du = dJ / dF_iJ * dF_iJ / du = J * inv(F)_Ji * dF_iJ / du
  
   Real dJdU = 0.0;
   RankTwoTensor invF = _dgrad[nd].inverse();
   Real detF = _dgrad[nd].det();
   for (unsigned int i = 0; i < 3; ++i)
     for (unsigned int J = 0; J < 3; ++J)
     {
       if (component == 0)
         dJdU += detF * invF(J, i) * _ddgraddu[nd](i, J);
       else if (component == 1)
         dJdU += detF * invF(J, i) * _ddgraddv[nd](i, J);
       else if (component == 2)
         dJdU += detF * invF(J, i) * _ddgraddw[nd](i, J);
     }
  
   return dJdU;
 }

Referenced by computeLocalJacobian(), and computeLocalOffDiagJacobian().

◆ computeDSDFhat()

RankFourTensor FiniteStrainMechanicsNOSPD::computeDSDFhat ( unsigned int nd )

protected

Function to compute derivative of stress with respect to derived deformation gradient.

Parameters

nd	The local index of element node (either 1 or 2 for Edge2 element)

Returns: The calculated derivative

Definition at line 307 of file FiniteStrainMechanicsNOSPD.C.

 {
   // compute the derivative of stress w.r.t the Fhat for finite strain
   RankTwoTensor I(RankTwoTensor::initIdentity);
   RankFourTensor dSdFhat;
   dSdFhat.zero();
  
   // first calculate the derivative of incremental Cauchy stress w.r.t the inverse of Fhat
   // Reference: M. M. Rashid (1993), Incremental Kinematics for finite element applications, IJNME
   RankTwoTensor S_inc = _Jacobian_mult[nd] * _E_inc[nd];
   RankFourTensor Tp1;
   Tp1.zero();
   for (unsigned int i = 0; i < 3; ++i)
     for (unsigned int j = 0; j < 3; ++j)
       for (unsigned int k = 0; k < 3; ++k)
         for (unsigned int l = 0; l < 3; ++l)
           for (unsigned int m = 0; m < 3; ++m)
             for (unsigned int n = 0; n < 3; ++n)
               for (unsigned int r = 0; r < 3; ++r)
                 Tp1(i, j, k, l) +=
                     S_inc(m, n) *
                         (_R_inc[nd](j, n) * (0.5 * I(k, m) * I(i, l) - I(m, l) * _R_inc[nd](i, k) +
                                              0.5 * _R_inc[nd](i, k) * _R_inc[nd](m, l)) +
                          _R_inc[nd](i, m) * (0.5 * I(k, n) * I(j, l) - I(n, l) * _R_inc[nd](j, k) +
                                              0.5 * _R_inc[nd](j, k) * _R_inc[nd](n, l))) -
                     _R_inc[nd](l, m) * _R_inc[nd](i, n) * _R_inc[nd](j, r) *
                         _Jacobian_mult[nd](n, r, m, k);
  
   // second calculate derivative of inverse of Fhat w.r.t Fhat
   // d(inv(Fhat)_kl)/dFhat_mn = - inv(Fhat)_km * inv(Fhat)_nl
   // the bases are gk, gl, gm, gn, indictates the inverse rather than the inverse transpose
  
   RankFourTensor Tp2;
   Tp2.zero();
   RankTwoTensor invFhat = (_dgrad[nd] * _dgrad_old[nd].inverse()).inverse();
   for (unsigned int k = 0; k < 3; ++k)
     for (unsigned int l = 0; l < 3; ++l)
       for (unsigned int m = 0; m < 3; ++m)
         for (unsigned int n = 0; n < 3; ++n)
           Tp2(k, l, m, n) += -invFhat(k, m) * invFhat(n, l);
  
   // assemble two calculated quantities to form the derivative of Cauchy stress w.r.t
   // Fhat
   dSdFhat = Tp1 * Tp2;
  
   return dSdFhat;
 }

Referenced by computeDSDU(), computeNonlocalJacobian(), and computePDNonlocalOffDiagJacobian().

◆ computeDSDU()

RankTwoTensor FiniteStrainMechanicsNOSPD::computeDSDU	(	unsigned int	component,
		unsigned int	nd
	)

overrideprotectedvirtual

Function to compute derivative of stress with respect to displacements.

Parameters

component	The index of displacement component
nd	The local index of element node (either 1 or 2 for Edge2 element)

Returns: The calculated derivative

Reimplemented from MechanicsBaseNOSPD.

Definition at line 282 of file FiniteStrainMechanicsNOSPD.C.

 {
   // compute the derivative of stress w.r.t the solution components for finite strain
   RankTwoTensor dSdU;
  
   // fetch the derivative of stress w.r.t the Fhat
   RankFourTensor DSDFhat = computeDSDFhat(nd);
  
   // third calculate derivative of Fhat w.r.t solution components
   RankTwoTensor Tp3;
   if (component == 0)
     Tp3 = _dgrad_old[nd].inverse() * _ddgraddu[nd];
   else if (component == 1)
     Tp3 = _dgrad_old[nd].inverse() * _ddgraddv[nd];
   else if (component == 2)
     Tp3 = _dgrad_old[nd].inverse() * _ddgraddw[nd];
  
   // assemble the fetched and calculated quantities to form the derivative of Cauchy stress w.r.t
   // solution components
   dSdU = DSDFhat * Tp3;
  
   return dSdU;
 }

Referenced by computeLocalJacobian(), and computeLocalOffDiagJacobian().

◆ computeLocalJacobian()

void FiniteStrainMechanicsNOSPD::computeLocalJacobian ( )

overrideprotectedvirtual

Definition at line 64 of file FiniteStrainMechanicsNOSPD.C.

 {
   // excludes dTi/dUj and dTj/dUi contribution which was considered as nonlocal contribution
   std::vector<RankTwoTensor> dPxdUx(_nnodes);
   for (unsigned int nd = 0; nd < _nnodes; ++nd)
     dPxdUx[nd] = computeDJDU(_component, nd) * _stress[nd] * _dgrad[nd].inverse().transpose() +
                  _dgrad[nd].det() * computeDSDU(_component, nd) * _dgrad[nd].inverse().transpose() +
                  _dgrad[nd].det() * _stress[nd] * computeDinvFTDU(_component, nd);
  
   for (_i = 0; _i < _test.size(); ++_i)
     for (_j = 0; _j < _phi.size(); ++_j)
       _local_ke(_i, _j) += (_i == 0 ? -1 : 1) * _multi[_j] *
                            (dPxdUx[_j] * _shape2[_j].inverse()).row(_component) * _origin_vec_ij *
                            _bond_status_ij;
 }

◆ computeLocalOffDiagJacobian()

void FiniteStrainMechanicsNOSPD::computeLocalOffDiagJacobian ( unsigned int )

overrideprotectedvirtual

Function to compute local contribution to the off-diagonal Jacobian at the current nodes.

Parameters

coupled_component The coupled variable number

Reimplemented from MechanicsBasePD.

Definition at line 169 of file FiniteStrainMechanicsNOSPD.C.

 {
   _local_ke.zero();
   if (coupled_component == 3)
   {
     std::vector<RankTwoTensor> dSdT(_nnodes);
     for (unsigned int nd = 0; nd < _nnodes; ++nd)
       for (unsigned int es = 0; es < _deigenstrain_dT.size(); ++es)
         dSdT[nd] = -_dgrad[nd].det() * _Jacobian_mult[nd] * (*_deigenstrain_dT[es])[nd] *
                    _dgrad[nd].inverse().transpose();
  
     for (_i = 0; _i < _test.size(); ++_i)
       for (_j = 0; _j < _phi.size(); ++_j)
         _local_ke(_i, _j) += (_i == 0 ? -1 : 1) * _multi[_j] *
                              (dSdT[_j] * _shape2[_j].inverse()).row(_component) * _origin_vec_ij *
                              _bond_status_ij;
   }
   else
   {
     std::vector<RankTwoTensor> dPxdUy(_nnodes);
     for (unsigned int nd = 0; nd < _nnodes; ++nd)
       dPxdUy[nd] =
           computeDJDU(coupled_component, nd) * _stress[nd] * _dgrad[nd].inverse().transpose() +
           _dgrad[nd].det() * computeDSDU(coupled_component, nd) * _dgrad[nd].inverse().transpose() +
           _dgrad[nd].det() * _stress[nd] * computeDinvFTDU(coupled_component, nd);
  
     for (_i = 0; _i < _test.size(); ++_i)
       for (_j = 0; _j < _phi.size(); ++_j)
         _local_ke(_i, _j) += (_i == 0 ? -1 : 1) * _multi[_j] *
                              (dPxdUy[_j] * _shape2[_j].inverse()).row(_component) * _origin_vec_ij *
                              _bond_status_ij;
   }
 }

◆ computeLocalResidual()

void FiniteStrainMechanicsNOSPD::computeLocalResidual ( )

overrideprotectedvirtual

Definition at line 42 of file FiniteStrainMechanicsNOSPD.C.

 {
   // For finite strain formulation, the _stress tensor gotten from material class is the
   // Cauchy stress (Sigma). the first Piola-Kirchhoff stress (P) is then obtained as
   // P = J * Sigma * inv(F)^T.
   // Nodal force states are based on the first Piola-Kirchhoff stress tensors (P).
   // i.e., T = (J * Sigma * inv(F)^T) * inv(Shape) * xi * multi.
   // Cauchy stress is calculated as Sigma_n+1 = Sigma_n + R * (C * dt * D) * R^T
  
   std::vector<Real> nodal_force_comp(_nnodes);
   for (unsigned int nd = 0; nd < _nnodes; ++nd)
     nodal_force_comp[nd] = _multi[nd] *
                            ((_dgrad[nd].det() * _stress[nd] * _dgrad[nd].inverse().transpose()) *
                             _shape2[nd].inverse())
                                .row(_component) *
                            (nd == 0 ? 1 : -1) * _origin_vec_ij;
  
   _local_re(0) = -(nodal_force_comp[0] - nodal_force_comp[1]) * _bond_status_ij;
   _local_re(1) = -_local_re(0);
 }

◆ computeNonlocalJacobian()

void FiniteStrainMechanicsNOSPD::computeNonlocalJacobian ( )

overrideprotectedvirtual

Definition at line 81 of file FiniteStrainMechanicsNOSPD.C.

 {
   // includes dTi/dUj and dTj/dUi contributions
   for (unsigned int cur_nd = 0; cur_nd < _nnodes; ++cur_nd)
   {
     RankFourTensor dSdFhat = computeDSDFhat(cur_nd);
     RankTwoTensor invF = _dgrad[cur_nd].inverse();
     Real detF = _dgrad[cur_nd].det();
     // calculation of jacobian contribution to current_node's neighbors
     std::vector<dof_id_type> dof(_nnodes);
     dof[0] = _current_elem->node_ptr(cur_nd)->dof_number(_sys.number(), _var.number(), 0);
     std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(cur_nd));
     unsigned int nb =
         std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - cur_nd)) -
         neighbors.begin();
     std::vector<dof_id_type> dg_neighbors =
         _pdmesh.getDefGradNeighbors(_current_elem->node_id(cur_nd), nb);
     std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(cur_nd));
     for (unsigned int k = 0; k < dg_neighbors.size(); ++k)
     {
       Node * node_k = _pdmesh.nodePtr(neighbors[dg_neighbors[k]]);
       dof[1] = node_k->dof_number(_sys.number(), _var.number(), 0);
       Real vol_k = _pdmesh.getPDNodeVolume(neighbors[dg_neighbors[k]]);
  
       // obtain bond ik's origin vector
       RealGradient origin_vec_ijk = *node_k - *_pdmesh.nodePtr(_current_elem->node_id(cur_nd));
  
       RankTwoTensor dFdUk;
       dFdUk.zero();
       for (unsigned int j = 0; j < _dim; ++j)
         dFdUk(_component, j) =
             _horiz_rad[cur_nd] / origin_vec_ijk.norm() * origin_vec_ijk(j) * vol_k;
  
       dFdUk *= _shape2[cur_nd].inverse();
  
       RankTwoTensor dPxdUkx;
       // calculate dJ/du
       Real dJdU = 0.0;
       for (unsigned int i = 0; i < 3; ++i)
         for (unsigned int J = 0; J < 3; ++J)
           dJdU += detF * invF(J, i) * dFdUk(i, J);
  
       // calculate dS/du
       RankTwoTensor dSdU = dSdFhat * dFdUk * _dgrad_old[cur_nd].inverse();
  
       // calculate dinv(F)Tdu
       RankTwoTensor dinvFTdU;
       dinvFTdU.zero();
       for (unsigned int i = 0; i < 3; ++i)
         for (unsigned int J = 0; J < 3; ++J)
           for (unsigned int k = 0; k < 3; ++k)
             for (unsigned int L = 0; L < 3; ++L)
               dinvFTdU(i, J) += -invF(J, k) * invF(L, i) * dFdUk(k, L);
  
       // calculate the derivative of first Piola-Kirchhoff stress w.r.t displacements
       dPxdUkx = dJdU * _stress[cur_nd] * invF.transpose() + detF * dSdU * invF.transpose() +
                 detF * _stress[cur_nd] * dinvFTdU;
  
       // bond status for bond k
       Real bond_status_ijk =
           _bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[k]]));
  
       _local_ke.resize(_test.size(), _phi.size());
       _local_ke.zero();
       for (_i = 0; _i < _test.size(); ++_i)
         for (_j = 0; _j < _phi.size(); ++_j)
           _local_ke(_i, _j) = (_i == 0 ? -1 : 1) * (_j == 0 ? 0 : 1) * _multi[cur_nd] *
                               (dPxdUkx * _shape2[cur_nd].inverse()).row(_component) *
                               _origin_vec_ij * _bond_status_ij * bond_status_ijk;
  
       _assembly.cacheJacobianBlock(_local_ke, _ivardofs_ij, dof, _var.scalingFactor());
  
       if (_has_diag_save_in)
       {
         unsigned int rows = _test.size();
         DenseVector<Real> diag(rows);
         for (unsigned int i = 0; i < rows; ++i)
           diag(i) = _local_ke(i, i);
  
         Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
         for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
           _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
       }
     }
   }
 }

◆ computeOffDiagJacobian()

void MechanicsBasePD::computeOffDiagJacobian ( MooseVariableFEBase & jvar )

overridevirtualinherited

Definition at line 71 of file MechanicsBasePD.C.

 {
   prepare();
  
   if (jvar.number() == _var.number())
     computeJacobian();
   else
   {
     unsigned int coupled_component = 0;
     bool active = false;
  
     for (unsigned int i = 0; i < _dim; ++i)
       if (jvar.number() == _disp_var[i]->number())
       {
         coupled_component = i;
         active = true;
       }
  
     if (_temp_coupled && jvar.number() == _temp_var->number())
     {
       coupled_component = 3;
       active = true;
     }
  
     if (_out_of_plane_strain_coupled && jvar.number() == _out_of_plane_strain_var->number())
     {
       coupled_component = 4;
       active = true;
     }
  
     if (active)
     {
       DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar.number());
       _local_ke.resize(ke.m(), ke.n());
       _local_ke.zero();
  
       computeLocalOffDiagJacobian(coupled_component);
  
       ke += _local_ke;
  
       if (_use_full_jacobian)
         computePDNonlocalOffDiagJacobian(jvar.number(), coupled_component);
     }
   }
 }

◆ computePDNonlocalOffDiagJacobian()

void FiniteStrainMechanicsNOSPD::computePDNonlocalOffDiagJacobian	(	unsigned int	,
		unsigned int
	)

overrideprotectedvirtual

Function to compute nonlocal contribution to the off-diagonal Jacobian at the current nodes.

Parameters

jvar_num	The number of the first coupled variable
coupled_component	The component number of the second coupled variable

Reimplemented from MechanicsBasePD.

Definition at line 204 of file FiniteStrainMechanicsNOSPD.C.

 {
   if (coupled_component != 3 && coupled_component != 4)
   {
     for (unsigned int cur_nd = 0; cur_nd < _nnodes; ++cur_nd)
     {
       RankFourTensor dSdFhat = computeDSDFhat(cur_nd);
       RankTwoTensor invF = _dgrad[cur_nd].inverse();
       Real detF = _dgrad[cur_nd].det();
       // calculation of jacobian contribution to current_node's neighbors
       std::vector<dof_id_type> jvardofs_ijk(_nnodes);
       jvardofs_ijk[0] = _current_elem->node_ptr(cur_nd)->dof_number(_sys.number(), jvar_num, 0);
       std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(cur_nd));
       unsigned int nb =
           std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - cur_nd)) -
           neighbors.begin();
       std::vector<dof_id_type> dg_neighbors =
           _pdmesh.getDefGradNeighbors(_current_elem->node_id(cur_nd), nb);
       std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(cur_nd));
       for (unsigned int k = 0; k < dg_neighbors.size(); ++k)
       {
         Node * node_k = _pdmesh.nodePtr(neighbors[dg_neighbors[k]]);
         jvardofs_ijk[1] = node_k->dof_number(_sys.number(), jvar_num, 0);
         Real vol_k = _pdmesh.getPDNodeVolume(neighbors[dg_neighbors[k]]);
  
         // obtain bond k's origin vector
         RealGradient origin_vec_ijk = *node_k - *_pdmesh.nodePtr(_current_elem->node_id(cur_nd));
  
         RankTwoTensor dFdUk;
         dFdUk.zero();
         for (unsigned int j = 0; j < _dim; ++j)
           dFdUk(coupled_component, j) =
               _horiz_rad[cur_nd] / origin_vec_ijk.norm() * origin_vec_ijk(j) * vol_k;
  
         dFdUk *= _shape2[cur_nd].inverse();
  
         RankTwoTensor dPxdUky;
         // calculate dJ/du
         Real dJdU = 0.0;
         for (unsigned int i = 0; i < 3; ++i)
           for (unsigned int J = 0; J < 3; ++J)
             dJdU += detF * invF(J, i) * dFdUk(i, J);
  
         // calculate dS/du
         RankTwoTensor dSdU = dSdFhat * dFdUk * _dgrad_old[cur_nd].inverse();
  
         // calculate dinv(F)Tdu
         RankTwoTensor dinvFTdU;
         dinvFTdU.zero();
         for (unsigned int i = 0; i < 3; ++i)
           for (unsigned int J = 0; J < 3; ++J)
             for (unsigned int k = 0; k < 3; ++k)
               for (unsigned int L = 0; L < 3; ++L)
                 dinvFTdU(i, J) += -invF(J, k) * invF(L, i) * dFdUk(k, L);
  
         // calculate the derivative of first Piola-Kirchhoff stress w.r.t displacements
         dPxdUky = dJdU * _stress[cur_nd] * invF.transpose() + detF * dSdU * invF.transpose() +
                   detF * _stress[cur_nd] * dinvFTdU;
  
         // bond status for bond k
         Real bond_status_ijk =
             _bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[k]]));
  
         _local_ke.zero();
         for (_i = 0; _i < _test.size(); ++_i)
           for (_j = 0; _j < _phi.size(); ++_j)
             _local_ke(_i, _j) = (_i == 0 ? -1 : 1) * (_j == 0 ? 0 : 1) * _multi[cur_nd] *
                                 (dPxdUky * _shape2[cur_nd].inverse()).row(_component) *
                                 _origin_vec_ij * _bond_status_ij * bond_status_ijk;
  
         _assembly.cacheJacobianBlock(_local_ke, _ivardofs_ij, jvardofs_ijk, _var.scalingFactor());
       }
     }
   }
 }

◆ initialSetup()

void MechanicsBasePD::initialSetup ( )

overridevirtualinherited

Definition at line 47 of file MechanicsBasePD.C.

 {
   _orientation = &_assembly.getFE(FEType(), 1)->get_dxyzdxi();
 }

◆ prepare()

void MechanicsBasePD::prepare ( )

overridevirtualinherited

Definition at line 53 of file MechanicsBasePD.C.

 {
   PeridynamicsKernelBase::prepare();
  
   _ivardofs_ij.resize(_nnodes);
  
   for (unsigned int i = 0; i < _nnodes; ++i)
     _ivardofs_ij[i] = _current_elem->node_ptr(i)->dof_number(_sys.number(), _var.number(), 0);
  
   for (unsigned int i = 0; i < _dim; ++i)
     _cur_ori_ij(i) = _origin_vec_ij(i) + _disp_var[i]->getNodalValue(*_current_elem->node_ptr(1)) -
                      _disp_var[i]->getNodalValue(*_current_elem->node_ptr(0));
  
   _cur_len_ij = _cur_ori_ij.norm();
   _cur_ori_ij /= _cur_len_ij;
 }