--- /dev/null
+\documentclass{article}
+\usepackage{amsmath}
+\begin{document}
+
+For ``normal'' elements, we have that
+\begin{gather*}
+ \varphi(\mathbf x)
+ =
+ \hat\varphi(\mathbf{\hat x})
+\end{gather*}
+where $\mathbf x = \Phi(\mathbf{\hat x})$. Thus, the gradient of a shape
+function in real coordinates is given by
+\begin{gather*}
+ \nabla\varphi(\mathbf x)
+ =
+ \frac{\partial \mathbf{\hat x}}{\partial \mathbf{x}}
+ \hat\nabla\hat\varphi(\mathbf{\hat x})
+\end{gather*}
+where the matrix $\frac{\partial \mathbf{\hat x}}{\partial \mathbf{x}}$ is the
+gradient of the inverse of the mapping for which it can be shown that
+\begin{gather*}
+ \frac{\partial \mathbf{\hat x}}{\partial \mathbf{x}}
+ =
+ \frac{\partial \Phi^{-1}(\mathbf{x})}{\partial \mathbf{x}}
+ =
+ \left(\frac{\partial \Phi(\mathbf{\hat x})}{\partial \mathbf{\hat x}}\right)^{-1}
+ =
+ J^{-1}.
+\end{gather*}
+Note that $\Phi$ is a polynomial mapping and so $J=\frac{\partial
+ \Phi(\mathbf{\hat x})}{\partial \mathbf{\hat x}}$ is easy to
+compute. Computing $J^{-1}$ then involves inverting a $d\times d$ matrix.
+
+As a consequence, we have just shown that we can write
+\begin{gather*}
+ \nabla\phi(\mathbf x) =
+ J^{-1} \hat\nabla \hat\varphi(\mathbf{\hat x})
+\end{gather*}
+and we can generalize this formula for second derivatives:
+\begin{gather*}
+ \nabla^2\phi(\mathbf x) =
+ [J^{-1} \hat\nabla][J^{-1} \hat\nabla] \hat\varphi(\mathbf{\hat x}).
+\end{gather*}
+Note, however, that for non-affine mappings, $J$ is not a constant matrix and
+so we need to apply the chain rule to obtain
+\begin{align*}
+ [\nabla^2\phi(\mathbf x)]_{ij}
+ &=
+ [J^{-1} \hat\nabla]_{i}[J^{-1} \hat\nabla]_j \hat\varphi(\mathbf{\hat x})
+ \\
+ &=
+ [J^{-1}_{ik} \hat\partial_k][J^{-1}_{jl} \hat\partial_l] \hat\varphi(\mathbf{\hat x})
+ \\
+ &=
+ J^{-1}_{ik} J^{-1}_{jl} \hat\partial_k \hat\partial_l
+ \hat\varphi(\mathbf{\hat x})
+ +
+ J^{-1}_{ik} (\hat\partial_k J^{-1}_{jl}) \hat\partial_l \hat\varphi(\mathbf{\hat x}).
+\end{align*}
+The difficulty is the derivative of the inverse of the Jacobian,
+$J^{-1}$. Since $J$ is a polynomial itself, the entries of $J^{-1}$ are
+rational functions and computing their derivatives is, at least,
+awkward. However, we can use that
+\begin{gather*}
+ 0 = \hat\partial_k \mathbf I_{pl} = \hat\partial_k (J_{pm}J^{-1}_{ml})
+ =
+ (\hat\partial_k J_{pm})J^{-1}_{ml}
+ +
+ J_{pm}(\hat\partial_k J^{-1}_{ml}).
+\end{gather*}
+Multiplying this equation with $J^{-1}_{jp}$ from the left and summing over
+$p$ yields the formula
+\begin{gather*}
+ \delta_{jm}
+ \hat\partial_k J^{-1}_{ml}
+ =
+ -
+ J^{-1}_{jp}(\hat\partial_k J_{pm})J^{-1}_{ml},
+\end{gather*}
+that is
+\begin{gather*}
+ \hat\partial_k J^{-1}_{jl}
+ =
+ -
+ J^{-1}_{jp}(\hat\partial_k J_{pm})J^{-1}_{ml}.
+\end{gather*}
+Inserting this into the formula for the second derivative yields
+\begin{align*}
+ [\nabla^2\phi(\mathbf x)]_{ij}
+ &=
+ J^{-1}_{ik} J^{-1}_{jl} \hat\partial_k \hat\partial_l
+ \hat\varphi(\mathbf{\hat x})
+ -
+ J^{-1}_{ik}
+ J^{-1}_{jp}(\hat\partial_k J_{pm})J^{-1}_{ml}
+ \hat\partial_l \hat\varphi(\mathbf{\hat x}).
+\end{align*}
+All of these terms are easy to compute since they only involve the already
+existing inverse Jacobian matrices $J^{-1}$ and the second derivatives
+$\hat\partial_k J_{pm} = \hat\partial_k \hat\partial_p \Phi_{m}$
+\marginpar{Is indeed $J_{pm} = \hat\partial_p \Phi_{m}$ or is it $J_{pm} = \hat\partial_m \Phi_{p}$?}
+of the polynomial forward mapping $\Phi$. Note that as expected, if the
+mapping is affine then $\hat\partial_k J_{pm}=0$ and the formula above reduces to
+the one expected.
+
+\end{document}