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 "CZMComputeDisplacementJumpTotalLagrangian.h"
      11             : #include "CohesiveZoneModelTools.h"
      12             : #include "FactorizedRankTwoTensor.h"
      13             : 
      14             : registerMooseObject("SolidMechanicsApp", CZMComputeDisplacementJumpTotalLagrangian);
      15             : registerMooseObject("SolidMechanicsApp", ADCZMComputeDisplacementJumpTotalLagrangian);
      16             : 
      17             : template <bool is_ad>
      18             : InputParameters
      19         132 : CZMComputeDisplacementJumpTotalLagrangianTempl<is_ad>::validParams()
      20             : {
      21         132 :   InputParameters params = CZMComputeDisplacementJumpBase<is_ad>::validParams();
      22         132 :   params.addClassDescription(
      23             :       "Compute the displacement jump increment across a czm interface in local "
      24             :       "coordinates for the Total Lagrangian kinematic formulation");
      25             : 
      26         132 :   return params;
      27           0 : }
      28             : 
      29             : template <bool is_ad>
      30          66 : CZMComputeDisplacementJumpTotalLagrangianTempl<
      31             :     is_ad>::CZMComputeDisplacementJumpTotalLagrangianTempl(const InputParameters & parameters)
      32             :   : CZMComputeDisplacementJumpBase<is_ad>(parameters),
      33          66 :     _F(this->template declareGenericPropertyByName<RankTwoTensor, is_ad>(_base_name + "F_czm")),
      34          66 :     _R(this->template declareGenericPropertyByName<RankTwoTensor, is_ad>(_base_name +
      35             :                                                                          "czm_rotation")),
      36          66 :     _czm_reference_rotation(this->template declareGenericPropertyByName<RankTwoTensor, is_ad>(
      37          66 :         _base_name + "czm_reference_rotation"))
      38             : {
      39             :   // initializing the displacement gradient vectors
      40         246 :   for (unsigned int i = 0; i < _ndisp; ++i)
      41             :   {
      42         180 :     _grad_disp.push_back(&this->template coupledGenericGradient<is_ad>("displacements", i));
      43             :     _grad_disp_neighbor.push_back(
      44         540 :         &this->template coupledGenericNeighborGradient<is_ad>("displacements", i));
      45             :   }
      46             : 
      47             :   // All others zero (so this will work naturally for 2D and 1D problems)
      48          84 :   for (unsigned int i = _ndisp; i < 3; i++)
      49             :   {
      50             :     if constexpr (is_ad)
      51             :     {
      52           0 :       _grad_disp.push_back(&_ad_grad_zero);
      53           0 :       _grad_disp_neighbor.push_back(&_ad_grad_zero);
      54             :     }
      55             :     else
      56             :     {
      57          18 :       _grad_disp.push_back(&_grad_zero);
      58          18 :       _grad_disp_neighbor.push_back(&_grad_zero);
      59             :     }
      60             :   }
      61          66 : }
      62             : 
      63             : template <bool is_ad>
      64             : void
      65         640 : CZMComputeDisplacementJumpTotalLagrangianTempl<is_ad>::initQpStatefulProperties()
      66             : {
      67         640 :   CZMComputeDisplacementJumpBase<is_ad>::initQpStatefulProperties();
      68         640 :   _displacement_jump_global[_qp] = 0;
      69         640 :   _F[_qp].setToIdentity();
      70         640 :   _R[_qp].setToIdentity();
      71         640 : }
      72             : 
      73             : template <bool is_ad>
      74             : void
      75       39904 : CZMComputeDisplacementJumpTotalLagrangianTempl<is_ad>::computeLocalDisplacementJump()
      76             : {
      77       47680 :   _interface_displacement_jump[_qp] =
      78       39904 :       _czm_total_rotation[_qp].transpose() * _displacement_jump_global[_qp];
      79       39904 : }
      80             : 
      81             : template <bool is_ad>
      82             : void
      83       39904 : CZMComputeDisplacementJumpTotalLagrangianTempl<is_ad>::computeRotationMatrices()
      84             : {
      85       39904 :   _czm_reference_rotation[_qp] = CohesiveZoneModelTools::computeReferenceRotation<is_ad>(
      86       39904 :       _normals[_qp], this->_mesh.dimension());
      87        7776 :   computeFandR();
      88       39904 :   _czm_total_rotation[_qp] = _R[_qp] * _czm_reference_rotation[_qp];
      89       39904 : }
      90             : 
      91             : template <bool is_ad>
      92             : void
      93       39904 : CZMComputeDisplacementJumpTotalLagrangianTempl<is_ad>::computeFandR()
      94             : {
      95       79808 :   const auto F = (GenericRankTwoTensor<is_ad>::Identity() +
      96             :                   GenericRankTwoTensor<is_ad>::initializeFromRows(
      97       39904 :                       (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]));
      98       39904 :   const auto F_neighbor =
      99             :       (GenericRankTwoTensor<is_ad>::Identity() +
     100       39904 :        GenericRankTwoTensor<is_ad>::initializeFromRows((*_grad_disp_neighbor[0])[_qp],
     101       39904 :                                                        (*_grad_disp_neighbor[1])[_qp],
     102       39904 :                                                        (*_grad_disp_neighbor[2])[_qp]));
     103             : 
     104       47680 :   _F[_qp] = 0.5 * (F + F_neighbor);
     105             :   // According to Cody mooseError are always fatal, so nothing we can do about
     106             :   // them. The norm of the tensor might work, but there is the risk of an
     107             :   // unwanted overflow (I tried and it happens). So checking component by
     108             :   // component remains the only reasonable strategy. If someone finds a  better
     109             :   // way this could be changed in the future
     110      159616 :   for (unsigned int i = 0; i < 3; i++)
     111      478848 :     for (unsigned int j = 0; j < 3; j++)
     112      359136 :       if (!std::isfinite(MetaPhysicL::raw_value(_F[_qp](i, j))))
     113           0 :         throw MooseException("CZMMaterialBaseIncremental _F is not finite, reducing time step");
     114             : 
     115             :   using FR2T = std::conditional_t<is_ad, ADFactorizedRankTwoTensor, FactorizedRankTwoTensor>;
     116       39904 :   const FR2T C = _F[_qp].transpose() * _F[_qp];
     117       39904 :   const auto Uinv = MathUtils::sqrt(C).inverse().get();
     118       39904 :   _R[_qp] = _F[_qp] * Uinv;
     119       39904 : }
     120             : 
     121             : template class CZMComputeDisplacementJumpTotalLagrangianTempl<false>;
     122             : template class CZMComputeDisplacementJumpTotalLagrangianTempl<true>;