Merge pull request idaholab#28473 from Kavan-Khaledi/shell_contravari…

…ant_base_vector

Modifing base vectors of shell element and added test example

modules/solid_mechanics/doc/content/modules/solid_mechanics/ShellElements.md

@@ -51,7 +51,19 @@ The shell element assumes a plane stress condition, i.e., the normal stress in t
 \tilde{C}_{ijkl} = (g^i \cdot \hat{e}_m)(g^j \cdot \hat{e}_n)(g^k \cdot \hat{e}_o)(g^l \cdot \hat{e}_p) \hat{C}_{mnop}
 \end{equation}
 
-where $g^i$ are the contravariant base vectors and $\hat{e}_i$ are the local cartesian coordinate system computed using the covariant base vectors, i.e.,
+where $g^i$ are the contravariant base vectors. The contravariant base vectors $g^1$, $g^2$ and $g^3$ are determined from the covariant base vectors $g_1$, $g_2$ and $g_3$ using the following equations:
+\begin{equation}
+g^1=\dfrac{g_2 \times g_3}{g_1 \cdot (g_2 \times g_3)}
+\end{equation}
+
+\begin{equation}
+g^2=\dfrac{g_3 \times g_1}{g_1 \cdot (g_2 \times g_3)}
+\end{equation}
+
+\begin{equation}
+g^3=\dfrac{g_1 \times g_2}{g_1 \cdot (g_2 \times g_3)}
+\end{equation}
+These contravariant base vectors satisfy the orthogonality condition: $g^i \cdot g_j =\delta_{ij}$, where $\delta_{ij}$ is Kronecker delta. The $\hat{e}_i$ vectors are the local cartesian coordinate system computed using the covariant base vectors, i.e.,
 
 \begin{equation}
 \hat{e}_3 = \frac{g_3}{|g_3|}; \hat{e}_1 = \frac{g_2 \times \hat{e}^3}{|g_2 \times \hat{e}^3|}; \hat{e}_2 = \hat{e}_3 \times \hat{e}_1

modules/solid_mechanics/src/materials/ADComputeIncrementalShellStrain.C

@@ -353,6 +353,7 @@ ADComputeIncrementalShellStrain::computeGMatrix()
         gmn(0, 2) += (*_dxyz_dxi[j])[i](component) * (*_dxyz_dzeta[j])[i](component);
         gmn(1, 2) += (*_dxyz_deta[j])[i](component) * (*_dxyz_dzeta[j])[i](component);
 
+        // why are we dropping derivatives here?!
         J(0, component) = MetaPhysicL::raw_value((*_dxyz_dxi[j])[i](component));
         J(1, component) = MetaPhysicL::raw_value((*_dxyz_deta[j])[i](component));
         J(2, component) = MetaPhysicL::raw_value((*_dxyz_dzeta[j])[i](component));
@@ -384,6 +385,7 @@ ADComputeIncrementalShellStrain::computeGMatrix()
       J(2, 1) /= normz;
       J(2, 2) /= normz;
 
+      // _transformation_matrix is an AD property, but we're not setting the derivatives!
       (*_transformation_matrix)[i] = J;
 
       // calculate ge
@@ -393,28 +395,30 @@ ADComputeIncrementalShellStrain::computeGMatrix()
       ADRealVectorValue e2 = e3.cross(e1);
       e2 /= e2.norm();
 
-      ADRankTwoTensor local_rotation_mat;
-      local_rotation_mat(0, 0) = e1(0);
-      local_rotation_mat(0, 1) = e1(1);
-      local_rotation_mat(0, 2) = e1(2);
-      local_rotation_mat(1, 0) = e2(0);
-      local_rotation_mat(1, 1) = e2(1);
-      local_rotation_mat(1, 2) = e2(2);
-      local_rotation_mat(2, 0) = e3(0);
-      local_rotation_mat(2, 1) = e3(1);
-      local_rotation_mat(2, 2) = e3(2);
-
-      ADRankTwoTensor gmninv = local_rotation_mat.transpose() * gmninv_temp * local_rotation_mat;
-
-      (*_ge[j])[i](0, 0) = (gmninv * (*_dxyz_dxi[j])[i]) * e1;
-      (*_ge[j])[i](0, 1) = (gmninv * (*_dxyz_dxi[j])[i]) * e2;
-      (*_ge[j])[i](0, 2) = (gmninv * (*_dxyz_dxi[j])[i]) * e3;
-      (*_ge[j])[i](1, 0) = (gmninv * (*_dxyz_deta[j])[i]) * e1;
-      (*_ge[j])[i](1, 1) = (gmninv * (*_dxyz_deta[j])[i]) * e2;
-      (*_ge[j])[i](1, 2) = (gmninv * (*_dxyz_deta[j])[i]) * e3;
-      (*_ge[j])[i](2, 0) = (gmninv * (*_dxyz_dzeta[j])[i]) * e1;
-      (*_ge[j])[i](2, 1) = (gmninv * (*_dxyz_dzeta[j])[i]) * e2;
-      (*_ge[j])[i](2, 2) = (gmninv * (*_dxyz_dzeta[j])[i]) * e3;
+      // Calculate the contravariant base vectors g^0, g^1, g^2
+      // The base vectors for the strain tensor in the convected coordinate
+      // are g_0, g_1, g_2 (g_i=dx/dr_i)
+      // The following contravariant base vectors have the property:
+      // g^i*g_j= 1 if {i=j} otherwise g^i*g_j= 0
+
+      const auto denom = (*_dxyz_dxi[j])[i] * (*_dxyz_deta[j])[i].cross((*_dxyz_dzeta[j])[i]);
+      const auto gi0 = (*_dxyz_deta[j])[i].cross((*_dxyz_dzeta[j])[i]) / denom;
+      const auto gi1 = (*_dxyz_dzeta[j])[i].cross((*_dxyz_dxi[j])[i]) / denom;
+      const auto gi2 = (*_dxyz_dxi[j])[i].cross((*_dxyz_deta[j])[i]) / denom;
+
+      // Calculate the rotation matrix for the elasticity tensor that maps
+      //  the strain tensor (with g_1, g2_, g_3 base vectors) to
+      //  the stress tensor (with base vectors e1, e2, e3)
+
+      (*_ge[j])[i](0, 0) = gi0 * e1;
+      (*_ge[j])[i](0, 1) = gi0 * e2;
+      (*_ge[j])[i](0, 2) = gi0 * e3;
+      (*_ge[j])[i](1, 0) = gi1 * e1;
+      (*_ge[j])[i](1, 1) = gi1 * e2;
+      (*_ge[j])[i](1, 2) = gi1 * e3;
+      (*_ge[j])[i](2, 0) = gi2 * e1;
+      (*_ge[j])[i](2, 1) = gi2 * e2;
+      (*_ge[j])[i](2, 2) = gi2 * e3;
     }
   }
 }