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             : #pragma once
      11             : 
      12             : #include "RankFourTensor.h"
      13             : #include "RankTwoTensor.h"
      14             : #include "libmesh/vector_value.h"
      15             : #include "MooseTypes.h"
      16             : #include "RotationMatrix.h"
      17             : 
      18             : namespace CohesiveZoneModelTools
      19             : {
      20             : template <bool is_ad = false>
      21             : GenericRealVectorValue<is_ad>
      22             : computedVnormdV(const GenericRealVectorValue<is_ad> & V)
      23             : {
      24             :   return V / V.norm();
      25             : }
      26             : 
      27             : /// compute the derivative of the rotation matrix, R=FU^-1, using Chen and Wheeler 1993
      28             : template <bool is_ad = false>
      29             : GenericRankFourTensor<is_ad>
      30       32128 : computedRdF(const GenericRankTwoTensor<is_ad> & R, const GenericRankTwoTensor<is_ad> & U)
      31             : {
      32       32128 :   const auto Uhat = U.trace() * RankTwoTensor::Identity() - U;
      33             :   unsigned int k, l, m, n, p, q;
      34       32128 :   const auto Uhat_det = Uhat.det();
      35             : 
      36       32128 :   GenericRankFourTensor<is_ad> dR_dF;
      37      128512 :   for (k = 0; k < 3; k++)
      38      385536 :     for (l = 0; l < 3; l++)
      39     1156608 :       for (m = 0; m < 3; m++)
      40     3469824 :         for (n = 0; n < 3; n++)
      41             :         {
      42     2602368 :           dR_dF(k, l, m, n) = 0.;
      43    10409472 :           for (p = 0; p < 3; p++)
      44    31228416 :             for (q = 0; q < 3; q++)
      45    23421312 :               dR_dF(k, l, m, n) +=
      46    23421312 :                   R(k, p) * (Uhat(p, q) * R(m, q) * Uhat(n, l) - Uhat(p, n) * R(m, q) * Uhat(q, l));
      47             : 
      48     2602368 :           dR_dF(k, l, m, n) /= Uhat_det;
      49             :         }
      50             : 
      51       32128 :   return dR_dF;
      52             : }
      53             : 
      54             : /// compute the derivative of F^-1 w.r.t. to F
      55             : RankFourTensor computedFinversedF(const RankTwoTensor & F_inv);
      56             : 
      57             : /// compute the area ratio betweeen the deformed and undeformed configuration, and its derivatives w.r.t. the deformation gradient, F
      58             : ///@{
      59             : template <bool is_ad = false>
      60             : GenericReal<is_ad>
      61       32128 : computeAreaRatio(const RankTwoTensor & FinvT, const Real & J, const RealVectorValue & N)
      62             : {
      63       32128 :   return J * (FinvT * N).norm();
      64             : }
      65             : 
      66             : template <bool is_ad = false>
      67             : GenericRankTwoTensor<is_ad>
      68       32128 : computeDAreaRatioDF(const GenericRankTwoTensor<is_ad> & FinvT,
      69             :                     const GenericRealVectorValue<is_ad> & N,
      70             :                     const GenericReal<is_ad> & J,
      71             :                     const GenericRankFourTensor<is_ad> & DFinv_DF)
      72             : {
      73       32128 :   const auto Fitr_N = FinvT * N;
      74       32128 :   const auto Fitr_N_norm = Fitr_N.norm();
      75       32128 :   GenericRankTwoTensor<is_ad> R2temp;
      76      128512 :   for (unsigned int l = 0; l < 3; l++)
      77      385536 :     for (unsigned int m = 0; m < 3; m++)
      78             :     {
      79      289152 :       R2temp(l, m) = 0;
      80     1156608 :       for (unsigned int i = 0; i < 3; i++)
      81     3469824 :         for (unsigned int j = 0; j < 3; j++)
      82     2602368 :           R2temp(l, m) += Fitr_N(i) * DFinv_DF(j, i, l, m) * N(j);
      83             : 
      84      289152 :       R2temp(l, m) *= J / Fitr_N_norm;
      85             :     }
      86             : 
      87       32128 :   return J * FinvT * Fitr_N_norm + R2temp;
      88             : }
      89             : ///@}
      90             : 
      91             : /// compute the normal componets of a vector
      92             : template <bool is_ad = false>
      93             : GenericRealVectorValue<is_ad>
      94      176896 : computeNormalComponents(const GenericRealVectorValue<is_ad> & normal,
      95             :                         const GenericRealVectorValue<is_ad> & vector)
      96             : {
      97      176896 :   return (normal * vector) * normal;
      98             : }
      99             : 
     100             : /// compute the tangent componets of a vector
     101             : template <bool is_ad = false>
     102             : GenericRealVectorValue<is_ad>
     103       88448 : computeTangentComponents(const GenericRealVectorValue<is_ad> & normal,
     104             :                          const GenericRealVectorValue<is_ad> & vector)
     105             : {
     106       88448 :   return vector - computeNormalComponents<is_ad>(normal, vector);
     107             : }
     108             : 
     109             : /// compute the czm reference rotation based on the normal in the undeformed configuration and the mesh dimension
     110             : template <bool is_ad = false>
     111             : GenericRankTwoTensor<is_ad>
     112       81104 : computeReferenceRotation(const GenericRealVectorValue<is_ad> & normal,
     113             :                          const unsigned int mesh_dimension)
     114             : {
     115       81104 :   GenericRankTwoTensor<is_ad> rot;
     116       81104 :   switch (mesh_dimension)
     117             :   {
     118       60480 :     case 3:
     119       76032 :       rot = RotationMatrix::rotVec1ToVec2<is_ad>(GenericRealVectorValue<false>(1, 0, 0), normal);
     120       76032 :       break;
     121        5040 :     case 2:
     122        5040 :       rot = RotationMatrix::rotVec2DToX<is_ad>(normal).transpose();
     123        5040 :       break;
     124          32 :     case 1:
     125          32 :       rot = RankTwoTensor::Identity();
     126          32 :       break;
     127           0 :     default:
     128           0 :       mooseError("computeReferenceRotation: mesh_dimension value should be 1, 2 or, 3. You "
     129             :                  "provided " +
     130             :                  std::to_string(mesh_dimension));
     131             :   }
     132       81104 :   return rot;
     133             : }
     134             : } // namespace CohesiveZoneModelTools