$\newcommand{\Dev}{\operatorname{Dev}}$ $\newcommand{\dev}{\operatorname{dev}}$ $\newcommand{\grad}{\operatorname{grad}}$ $\newcommand{\Grad}{\operatorname{Grad}}$ $\newcommand{\divg}{\operatorname{div}}$ $\newcommand{\Ei}{\operatorname{Ei}}$ $\newcommand{\tr}{\operatorname{tr}}$ $\newcommand{\Divg}{\operatorname{Div}}$ $\newcommand{\cay}{\operatorname{cay}}$ $\newcommand{\rot}{\operatorname{rot}}$

We create a shell formulation by reducing a 3D element interpolation which is linear along $\xi_{3}$ . The nodal positions at the back of the shell ($\xi_{3}=-1$ ) are denoted by $\mathbf{y}_{a}$ and those on the front of the shell ($\xi_{3}=+1$ ) are denoted by $\mathbf{x}_{a}$ , thus $$\begin{equation} \mathbf{x}\left(\xi_{i}\right)=\sum_{a}M_{a}\left(\xi_{1},\xi_{2}\right)\left(\frac{1+\xi_{3}}{2}\mathbf{x}_{a}+\frac{1-\xi_{3}}{2}\mathbf{y}_{a}\right)=\sum_{a}M_{a}\left(\xi_{1},\xi_{2}\right)\left(\mathbf{x}_{a}-\frac{1-\xi_{3}}{2}\mathbf{d}_{a}\right)\label{eq:solid-element-interpolation} \end{equation}$$ The vector from $\mathbf{y}_{a}$ to $\mathbf{x}_{a}$ is the director, $\mathbf{d}_{a}$ , $$\mathbf{d}_{a}=\mathbf{x}_{a}-\mathbf{y}_{a}$$ From this relation we can get the shell covariant basis vectors, $$\begin{equation} \begin{aligned}\mathbf{g}_{\alpha}\left(\xi_{i}\right) & =\frac{\partial\mathbf{x}}{\partial\xi_{\alpha}}=\sum_{a}\frac{\partial M_{a}}{\partial\xi_{\alpha}}\left(\frac{1+\xi_{3}}{2}\mathbf{x}_{a}+\frac{1-\xi_{3}}{2}\mathbf{y}_{a}\right)=\sum_{a}\frac{\partial M_{a}}{\partial\xi_{\alpha}}\left(\mathbf{x}_{a}-\frac{1-\xi_{3}}{2}\mathbf{d}_{a}\right)\\ \mathbf{g}_{3}\left(\xi_{i}\right) & =\frac{\partial\mathbf{x}}{\partial\xi_{3}}=\sum_{a}\frac{1}{2}M_{a}\left(\xi_{1},\xi_{2}\right)\left(\mathbf{x}_{a}-\mathbf{y}_{a}\right)=\sum_{a}\frac{1}{2}M_{a}\left(\xi_{1},\xi_{2}\right)\mathbf{d}_{a} \end{aligned} \label{eq:shell-covariant-basis} \end{equation}$$ from which we may evaluate the contravariant basis vectors $\mathbf{g}^{i}$ . Let the front-face and back-face displacements be denoted by $\mathbf{u}$ and $\mathbf{w}$ , respectively. It follows that $\mathbf{x}_{a}=\mathbf{X}_{a}+\mathbf{u}_{a}$ and $\mathbf{y}_{a}=\mathbf{Y}_{a}+\mathbf{w}_{a}$ , where $\mathbf{X}_{a}$ represents the shell nodal positions in the reference configuration, provided as nodal coordinates in the input file, and $\mathbf{Y}_{a}=\mathbf{X}_{a}-\mathbf{D}_{a}$ is evaluated from the user-defined referential shell thickness, and the surface surface normals evaluated at each node. If the shell surface is not planar in the reference configuration, users must be careful to select shell thicknesses that don't produce inverted elements (negative Jacobians) as a result of this extrapolation.