22. Analysis of the generalized and lifted shape operator

In this notebook we will provide a proof roadmap on how to rigorously analyze the generalized and lifted Weingarten tensor using numerical analysis and derive convergence rates. We follow [Gopalakrishnan, Neunteufel. Analysis of the generalized shape operator for surfaces (in preparation)].

22.1. The comparison problem

The exact shape operator and the discrete lifted shape operator do not naturally live on the same surface. Also, there are several difficulties we have to take care of:

• \(W\) lives on the smooth surface \(\mathcal S\),

• \(\kappa_h\) lives on the discrete surface \(\mathcal S_h\),

• the generalized operator is a distributional functional,

• the involved quantities (e.g. the signed dihedral angle) are highly nonlinear.

During the analysis we need to compare functionals across changing geometries.

22.2. Step 1: connect the two surfaces by interpolation of embeddings

Introduce a family of embeddings

\[ \Phi_h^t = \bar{\Phi}_h + t(\Phi_h - \bar{\Phi}_h), \qquad 0\le t\le 1, \]

where \(\bar{\Phi}_h:\hat{\Omega}\to\mathcal{S}\) and \(\Phi_h:\hat{\Omega}\to\mathcal{S}_h\) are the embeddings of the exact and approximated surface living on the same reference/parameter domain \(\hat{\Omega}\subset \mathbb{R}^2\).

This connects the exact triangulated surface to the discrete surface. The error is then written with the fundamental theorem of calculus into an integral formulation

\[ \widetilde W_h(\sigma_h) - \tilde{\bar{W}}_h(\sigma) = \int_0^1 \frac{d}{dt}\widetilde W(\Phi_h^t)(\sigma_t)\,dt . \]

Note, that \(\tilde{\bar{W}}_h(\sigma) = \int_{\mathcal{S}}\bar{W}:\sigma\).

It is important to note, that we have to pull back all quantities back to the common parameter domain \(\hat{\Omega}\), including the test functions \(\sigma\) in the Hellan-Herrmann-Johnson space.

Thus, in a first intermediate step, we need to investigate how the generalized Weingarten tensor behaves, if the surface (more precisely its embedding) varies.

22.3. Step 2: Linearization

Denote with \(\hat X=\frac{d}{dt}\Phi_h^t\) the perturbation vector of the surface and define \(X_t\circ\Phi_h^t=\hat{X}\). Further, define \(F=\hat{\nabla} \hat X\), \(J=\sqrt{\det(F^TF)}\).

Similar to shape optimization techniques, we can compute how the involved quantities change:

\[\begin{split} \begin{align*} &\frac{d}{dt} J = (\mathrm{div}_{\mathcal{S}} X)\circ\Phi\,J && \text{surface determinant}\\ &\frac{d}{dt}(\nu\circ\Phi) = -((\nabla_{\mathcal{S}}X)^T\nu)\circ\Phi && \text{surface normal}\\ &\frac{d}{dt}({\sigma}\circ\Phi) = (-2\mathrm{div}_{\mathcal{S}}(X){\sigma}+2\mathrm{sym}(\nabla_{\mathcal{S}}X{\sigma}))\circ\Phi && \text{HHJ test function}\\ &\frac{d}{dt}({\sigma}_{\mu\mu}\circ\Phi) = (-2(\nabla_{\mathcal{S}}X)_{\tau\tau}{\sigma}_{\mu\mu})\circ\Phi\\ &\frac{d}{dt}(\nabla_{\mathcal{S}}\nu\circ\Phi) = \left(\sum_{i=1}^3\nu_i\nabla_{\mathcal{S}}\nabla_{\mathcal{S}}X_i-(\nabla_{\mathcal{S}}X)^T\nabla_{\mathcal{S}}\nu + \nabla_{\mathcal{S}}\nu((\nabla_{\mathcal{S}}X^T\nu)\otimes\nu -\nabla_{\mathcal{S}}X)\right)\circ\Phi && \text{Weingarten tensor}\\ &\frac{d}{dt}(\sphericalangle_{\mathrm{sgn}}(\nu_L,\nu_R)\circ\Phi) = [\![ (\nabla_{\mathcal{S}}X)_{\nu\mu}]\!]\circ\Phi && \text{dihedral angle} \end{align*} \end{split}\]

Thus, we derive at the following linearization of the generalized shape operator

\[ \widetilde{W}({\sigma})= \sum_{T\in\mathcal T}\int_T W_T:\sigma\,ds+\sum_{E\in\mathcal E_{\rm int}}\int_E \sphericalangle_{\mathrm{sgn}}(\nu_L,\nu_R)\,\sigma_{\mu\mu}\,dl \]

with \(X=\dot{X}\):

\[ \begin{align*} \frac{d}{dt}\widetilde{W}({\sigma}) = a(\Phi;{\sigma},X)+b(\Phi;{\sigma},X), \end{align*} \]

where with \(\mathcal{H}_{\nu}(X) = \sum_{i=1}^3\mathrm{hesse}(X_i)\nu_i\)

\[\begin{split} \begin{align*} a(\Phi;{\sigma},X) &= \sum_{T\in \mathcal{T}}\int_T-\mathrm{div}(X)\nabla\nu:{\sigma} -\sum_{E\in\mathring{\mathcal{E}}}\int_E(\nabla X)_{\tau\tau}\sphericalangle_{\mathrm{sgn}}(\nu_L,\nu_R){\sigma}_{\mu\mu},\\ b(\Phi;{\sigma},X) &= \sum_{T\in\mathcal{T}}\int_T-\mathcal{H}_{\nu}(X):{\sigma} + \sum_{E\in\mathring{\mathcal{E}}}\int_E[\![ (\nabla X)_{\nu\mu} ]\!]_E{\sigma}_{\mu\mu}. \end{align*} \end{split}\]

The first bilinear form, \(a(\cdot;\cdot,\cdot)\) only contains first derivatives of \(X\) and the boundary term involving \(X\) is single-valued. The second bilinear for \(b(\cdot;\cdot,\cdot)\) contains bending-type terms involving second derivatives of \(X\) and is closely related to the HHJ method for Kirchhoff-Love shell bending. That is the analytic echo of the finite element choice and will be exploited in the numerical analysis.

22.4. Step 3: pull everything back

As motivated above we need to transform both bilinear forms back to the parameter domain \(\hat \Omega\) to perform the analysis.

\[\begin{split} \begin{align*} a(\Phi;{\sigma},X) &= \sum_{\hat{T}\in \hat{\mathcal{T}}}\int_{\hat{T}}-J^{-1}\mathrm{tr}(\hat{\nabla}\hat{X} F^{\dagger})\hat{S}:{\hat{\sigma}} -\sum_{\hat{E}\in\mathring{\hat{\mathcal{E}}}}\int_{\hat{E}}J_{\mathrm{bnd}}^{-3}(\hat{\nabla}_{\hat{\tau}} \hat{X})\cdot (F\hat{\tau})\sphericalangle_{\mathrm{sgn}}(\nu_L,\nu_R)\circ\Phi{\hat{\sigma}}_{\hat{\mu}\hat{\mu}},\\ b(\Phi;{\sigma},X) &= \sum_{\hat{T}\in \hat{\mathcal{T}}}\int_{\hat{T}}-\sum_{i=1}^3J^{-1}(F_1\times F_2)\left(\hat{\nabla}^2\hat{X}_i-\sum_{\alpha=1}^2(\hat{\nabla}_{\alpha}\hat{X}_i){\Gamma}^{\alpha}\right):{\hat{\sigma}} \\ &\quad+ \sum_{\hat{E}\in\mathring{\hat{\mathcal{E}}}}\int_{\hat{E}}[\![ \sum_{i=1}^3J^{-1}(F_1\times F_2)\cdot\left(\hat{\nabla}_{\hat{\mu}}\hat{X}_i-J_{\mathrm{bnd}}^2(F^TF)_{\hat{\tau}\hat{\mu}}\hat{\nabla}_{\hat{\tau}}\hat{X}_i\right) ]\!]_E{\hat{\sigma}}_{\hat{\mu}\hat{\mu}}. \end{align*} \end{split}\]

Above, we used the notation of the Christoffel symbols of second kind: \({\Gamma}^{\alpha}_{\beta\gamma} = \sum_{\ell=1}^3F^{\dagger}_{\alpha\ell}\hat{\nabla}_{\beta}F_{\ell \gamma}\). Further, \(\hat{S}=\sum_{i=1}^3\nu_i\circ\Phi\hat{\nabla}^2 \Phi_i\) arises from the pull-back of the Weingarten tensor \(\nabla_{\mathcal{S}}\nu\circ\Phi = F^{\dagger^T}\hat{S} F^{\dagger}\). The HHJ test function \(\sigma\) is transformed by the double Piola transformation \({\sigma}\circ\Phi = J^{-2}F{\hat{\sigma}} F^T\) and the normal vector reads \(\nu\circ\Phi = J^{-1}F_1\times F_2\).

22.5. Analysis of generalized shape operator

The \(H^{-1}\)-norm is well-suited for the analysis of the generalized shape operator as it consists of volume and codimension 1 facet (edge) contributions. Thus we will estimate

\[ \widetilde W_h(\sigma_h) - \tilde{\bar{W}}_h(\sigma) = \int_0^1 \left(a(\Phi_t;X_t,\sigma_t) + b(\Phi_t;X_t,\sigma_t)\right)\,dt, \]

by estimating the integrand, which are sufficiently linear to perform numerical analysis. Goal: estimate to get \(|a(\Phi_t;X_t,\sigma_t)|\le C \|\hat X\|\|\hat \sigma\|_{H^1}\) and \(|b(\Phi_t;X_t,\sigma_t)|\le C \|\hat X\|\|\hat \sigma\|_{H^1}\) with the weakest norm for \(\hat X\) as possible as \(\hat X=\Phi_h-\bar{\Phi}_h\) gives the convergence.

It is crucial to carefully estimate the involved nonlinear terms. For example one can proof for the signed dihedral angle that

\[ \begin{align*} \|\sphericalangle_{\mathrm{sgn}}(\nu_L,\nu_R)\|_{W^{1,\infty}}\le h\|\bar{\Phi}_h\|_{W^{\min\{k,2\},\infty}}. \end{align*} \]

We arrive at the following result Let \((\Phi_h)_{h>0}\in\mathrm{Lag}^k_h\) be a family of embeddings such that \(\|\Phi_h-\bar{\Phi}_h\|_{W^{1,\infty}}\!\to\! 0\). Then there holds

\[ \begin{align*} \|\widetilde{W}-\bar{W}\|_{H^{-1}} \le C(1+\max\limits_{\hat{T}\in\hat{\mathcal{T}}}h^{-1}_{\hat{T}}\|\Phi_h-\bar{\Phi}_h\|_{W^{\min\{k,2\},\infty}(\hat{T})})\|\Phi_h-\bar{\Phi}_h\|_{H^1}. \end{align*} \]

If the embedding \(\Phi_h\) has optimal approximation properties, we obtain the following convergence result

\[ \|\widetilde{W}-\bar{W}\|_{H^{-1}} \le C\,h^k. \]

Thus, we have shown that the dihedral angle used in discrete differential geometry always converges in the \(H^{-1}\)-norm.

22.6. Analysis of lifted generalized Weingarten tensor

Next, we investigate the convergence of the lifted generalized Weingarten tensor

\[ (\kappa_h,\tau_h)_{L^2(\mathcal T_h)} = \widetilde W_h(\tau_h). \]

For an optimal order interpolant \(\Phi_h\) of \(\bar{\Phi}_h\) we get the same convergence result: Let \((\Phi_h)_{h>0}\in\mathrm{Lag}^k_h\) be a family of embeddings and \({\kappa}_h\in M_h^{k-1}\) be the lifted Weingarten tensor. Then

\[ \|{\kappa}_h-\bar{W}\|_{H^{-1}} \le C\,h^{k+1}. \]

Note, that we need the technical assumption of a quasi-uniform triangulation as we use an \(H^1\)-like stability of an \(L^2\)-like projection operator including a global inverse estimate.

However, we can improve the convergence rates if the approximated surface comes from the canonical Lagrange interpolation operator, \(\Phi_h= \mathcal{I}^{\mathrm{Lag}^k}_h\bar{\Phi}_h\). This includes the important case of densely inscribed surfaces, i.e., the vertices of the triangulation lie on the exact surface. The key ingredient for the improved convergence is the magic orthogonality condition of the HHJ bilinear form

\[ b(\sigma_h, u-\mathcal{I}^{\mathrm{Lag}^k}_hu)=0,\qquad \forall \sigma_h \in M_h^{k-1}, \]

where

\[ b(\sigma, u) = \sum_T\int_T\sigma:\nabla^2 u\,dx -\int_{\partial T} \sigma_{nn}\partial_nu\,ds = \sum_T\int_T\mathrm{div}\mathrm{div}\sigma\, u\,dx -\int_{\partial T} (\mathrm{div}\sigma + \partial_t\sigma_{nt}\,)u\,ds +\sum_{V\in T} [\![\sigma_{nt}]\!]_V^T\,u(V). \]

Let \((\Phi_h)_{h>0}\in\mathrm{Lag}^k_h\) be a family of embeddings such that \(\Phi_h= \mathcal{I}^{\mathrm{Lag}^k}_h\bar{\Phi}_h\) for \(k\ge 1\). Let \({\kappa}_h\in M_h^{k-1}\) be the lifted Weingarten tensor. Then

\[ \|{\kappa}_h-\bar{W}\|_{H^{-1}} \le C\,h^{\textcolor{red}{k+1}}. \]

By using standard techniques, we can derive estimates in stronger norms, loosing the expected \(h\) powers. Thus, we have shown that the lifted Weingarten tensor always converges in the \(L^2\) norm, if the approximation stems from the canonical Lagrange interpolation operator.