Finite element procedures for computing normals and mean curvature on triangulated surfaces and their use for mesh refinement

TL;DR

This paper introduces finite element methods with stabilization techniques for accurately computing normals and mean curvature on triangulated surfaces, enabling adaptive mesh refinement and surface reconstruction.

Contribution

It presents a stabilized finite element approach for mean curvature and normals, improving convergence and enabling adaptive mesh refinement on triangulated surfaces.

Findings

Achieves first order L^2-convergence for mean curvature vector.

Provides a method for continuous normal recovery using L^2-projections.

Demonstrates improved surface reconstruction and mesh refinement results.

Abstract

In this paper we consider finite element approaches to computing the mean curvature vector and normal at the vertices of piecewise linear triangulated surfaces. In particular, we adopt a stabilization technique which allows for first order $L^2$-convergence of the mean curvature vector and apply this stabilization technique also to the computation of continuous, recovered, normals using $L^2$-projections of the piecewise constant face normals. Finally, we use our projected normals to define an adaptive mesh refinement approach to geometry resolution where we also employ spline techniques to reconstruct the surface before refinement. We compare or results to previously proposed approaches.

Figures

Tables

Table 1. (a) Torus with a = 1 𝑎 1 a=1

$h$	$L^2$	$L_{\mathrm{stab}}^2$	MWE	MWA	MWAAT	MWELR	MWRELR	MWSERL	DLLB
0.0527	0.0762	0.0360	0.1535	0.0723	0.1950	0.0820	0.0874	0.0722	0.5623
0.0245	0.0369	0.0077	0.0452	0.0222	0.0836	0.0292	0.0320	0.0224	0.6175
0.0113	0.0152	0.0018	0.0154	0.0080	0.0349	0.0115	0.0128	0.0081	0.3338
0.0057	0.0066	0.0005	0.0060	0.0032	0.0147	0.0047	0.0052	0.0033	0.2672
0.0029	0.0029	0.0002	0.0026	0.0014	0.0066	0.0021	0.0023	0.0014	0.1728

Table 2. (a) Torus with a = 1 𝑎 1 a=1

$h$	$L^2$	$L_{\mathrm{stab}}^2$	MWE	MWA	MWAAT	MWELR	MWRELR	MWSERL	DLLB
0.0527	0.0762	0.0360	0.1535	0.0723	0.1950	0.0820	0.0874	0.0722	0.5623
0.0245	0.0369	0.0077	0.0452	0.0222	0.0836	0.0292	0.0320	0.0224	0.6175
0.0113	0.0152	0.0018	0.0154	0.0080	0.0349	0.0115	0.0128	0.0081	0.3338
0.0057	0.0066	0.0005	0.0060	0.0032	0.0147	0.0047	0.0052	0.0033	0.2672
0.0029	0.0029	0.0002	0.0026	0.0014	0.0066	0.0021	0.0023	0.0014	0.1728

Table 3. (b) Torus with a = 4 𝑎 4 a=4

$h$	$L^2$	$L_{\mathrm{stab}}^2$	MWE	MWA	MWAAT	MWELR	MWRELR	MWSERL	DLLB
0.0657	0.2557	0.2557	0.5310	0.3618	0.4477	0.3828	0.3701	0.3884	0.3359
0.0294	0.0670	0.0666	0.1638	0.1190	0.1664	0.1231	0.1234	0.1209	0.1097
0.0132	0.0228	0.0150	0.0413	0.0287	0.0572	0.0314	0.0325	0.0288	0.0308
0.0067	0.0103	0.0038	0.0130	0.0085	0.0242	0.0100	0.0107	0.0085	0.0114
0.0034	0.0044	0.0009	0.0049	0.0028	0.0116	0.0036	0.0041	0.0027	0.0044

Table 4. (a) Torus with a = 1 𝑎 1 a=1

$h$	${ϵ}_{\mathrm{MWA}}-{ϵ}_{L^2}$	relative change	${ϵ}_{\mathrm{MWA}}-{ϵ}_{L_{\mathrm{stab}}^2}$	relative change
0.0527	-0.0039	-0.0542	0.0363	0.5022
0.0245	-0.0147	-0.6625	0.0145	0.6547
0.0113	-0.0072	-0.8959	0.0062	0.7764
0.0057	-0.0034	-1.0657	0.0027	0.8477
0.0029	-0.0015	-1.0784	0.0012	0.8876

Table 5. (a) Torus with a = 1 𝑎 1 a=1

$h$	${ϵ}_{\mathrm{MWA}}-{ϵ}_{L^2}$	relative change	${ϵ}_{\mathrm{MWA}}-{ϵ}_{L_{\mathrm{stab}}^2}$	relative change
0.0527	-0.0039	-0.0542	0.0363	0.5022
0.0245	-0.0147	-0.6625	0.0145	0.6547
0.0113	-0.0072	-0.8959	0.0062	0.7764
0.0057	-0.0034	-1.0657	0.0027	0.8477
0.0029	-0.0015	-1.0784	0.0012	0.8876

Table 6. (b) Torus with a = 4 𝑎 4 a=4

$h$	${ϵ}_{\mathrm{MWA}}-{ϵ}_{L^2}$	relative change	${ϵ}_{\mathrm{MWA}}-{ϵ}_{L_{\mathrm{stab}}^2}$	relative change
0.0657	0.1060	0.2931	0.1060	0.2931
0.0294	0.0521	0.4375	0.0524	0.4402
0.0132	0.0059	0.2064	0.0137	0.4781
0.0067	-0.0018	-0.2078	0.0048	0.5597
0.0034	-0.0016	-0.5827	0.0019	0.6837

Table 7. (a) Torus with a = 1 𝑎 1 a=1 .

$h$	$L^2$	$L_{\mathrm{stab}}^2$	MWE	MWA	MWAAT	MWELR	MWRELR	MWSERL	DLLB
0.0527	-	-	-	-	-	-	-	-	-
0.0245	0.9452	2.0151	1.5909	1.5386	1.1032	1.3449	1.3087	1.5214	-0.1219
0.0113	1.1524	1.8875	1.3999	1.3230	1.1332	1.2153	1.1924	1.3205	0.7987
0.0057	1.2124	1.8968	1.3641	1.3375	1.2636	1.3078	1.3005	1.3333	0.3249
0.0029	1.1844	1.6348	1.2098	1.1932	1.1687	1.1897	1.1841	1.2009	0.6328

Table 8. (a) Torus with a = 1 𝑎 1 a=1 .

$h$	$L^2$	$L_{\mathrm{stab}}^2$	MWE	MWA	MWAAT	MWELR	MWRELR	MWSERL	DLLB
0.0527	-	-	-	-	-	-	-	-	-
0.0245	0.9452	2.0151	1.5909	1.5386	1.1032	1.3449	1.3087	1.5214	-0.1219
0.0113	1.1524	1.8875	1.3999	1.3230	1.1332	1.2153	1.1924	1.3205	0.7987
0.0057	1.2124	1.8968	1.3641	1.3375	1.2636	1.3078	1.3005	1.3333	0.3249
0.0029	1.1844	1.6348	1.2098	1.1932	1.1687	1.1897	1.1841	1.2009	0.6328

Table 9. (b) Torus with a = 4 𝑎 4 a=4 .

$h$	$L^2$	$L_{\mathrm{stab}}^2$	MWE	MWA	MWAAT	MWELR	MWRELR	MWSERL	DLLB
0.0657	-	-	-	-	-	-	-	-	-
0.0294	1.6651	1.6711	1.4612	1.3812	1.2299	1.4099	1.3647	1.4506	1.3909
0.0132	1.3512	1.8705	1.7267	1.7826	1.3383	1.7122	1.6712	1.7999	1.5934
0.0067	1.1654	2.0310	1.6914	1.7818	1.2655	1.6789	1.6342	1.7863	1.4589
0.0034	1.2520	2.1265	1.4343	1.6452	1.0709	1.4678	1.3909	1.6451	1.3711

Table 10. (a) Torus with a = 1 𝑎 1 a=1

$h$	${\gamma }^{\ast }$	${ϵ}_{L^2}-{{ϵ}_{L_{\mathrm{stab}}^2}|}_{{\gamma }^{\ast }}$	relative change
0.0527	0.0147	0.0402	0.5277
0.0245	0.0371	0.0292	0.7923
0.0113	0.0693	0.0134	0.8821
0.0057	0.1057	0.0061	0.9262
0.0029	0.1623	0.0028	0.9459

Table 11. (a) Torus with a = 1 𝑎 1 a=1

$h$	${\gamma }^{\ast }$	${ϵ}_{L^2}-{{ϵ}_{L_{\mathrm{stab}}^2}|}_{{\gamma }^{\ast }}$	relative change
0.0527	0.0147	0.0402	0.5277
0.0245	0.0371	0.0292	0.7923
0.0113	0.0693	0.0134	0.8821
0.0057	0.1057	0.0061	0.9262
0.0029	0.1623	0.0028	0.9459

Table 12. (b) Torus with a = 4 𝑎 4 a=4

$h$	${\gamma }^{\ast }$	${ϵ}_{L^2}-{{ϵ}_{L_{\mathrm{stab}}^2}|}_{{\gamma }^{\ast }}$	relative change
0.0657	0	0	0
0.0294	0.0000	0.0003	0.0048
0.0132	0.0013	0.0078	0.3424
0.0067	0.0030	0.0065	0.6354
0.0034	0.0061	0.0035	0.8001

Table 13. Table 5: PN regular refinement of a torus with a = 4 𝑎 4 a=4 , initial mesh-size of 0.0321 and initial geometrical error of 0.0863.

$N_e$	$h$	${{ϵ}_{\mathrm{geom}}|}_{{𝒏}_{\mathrm{L2stab}}}$	${\Vert {𝒏}_f-{𝒏}_{\mathrm{L2stab}}\Vert }_{{\Sigma }_h}$	${{ϵ}_{\mathrm{geom}}|}_{{𝒏}_{\mathrm{MWA}}}$	${\Vert {𝒏}_f-{𝒏}_{\mathrm{MWA}}\Vert }_{{\Sigma }_h}$
1856	0.0900	0.0351	1.3934	0.0357	1.4351
7424	0.0451	0.0253	0.6976	0.0262	0.6981
29696	0.0226	0.0240	0.3453	0.0250	0.3354
118784	0.0113	0.0239	0.1766	0.0249	0.1619

Table 14. Table 6: Local refinement of a torus with a = 4 𝑎 4 a=4 , initial mesh-size of 0.0321 and initial geometrical error of 0.0863.

$N_e$	$h$	${{ϵ}_{\mathrm{geom}}|}_{{𝒏}_{\mathrm{L2stab}}}$	${\Vert {𝒏}_f-{𝒏}_{\mathrm{L2stab}}\Vert }_{{\Sigma }_h}$	$N_e$	$h$	${{ϵ}_{\mathrm{geom}}|}_{{𝒏}_{\mathrm{MWA}}}$	${\Vert {𝒏}_f-{𝒏}_{\mathrm{MWA}}\Vert }_{{\Sigma }_h}$
1278	0.1085	0.0470	1.6393	1278	0.1085	0.0472	1.6831
3470	0.0659	0.0286	0.9966	3390	0.0667	0.0296	1.0123
9420	0.0400	0.0247	0.5945	9068	0.0408	0.0258	0.5952
23862	0.0252	0.0237	0.3675	23060	0.0256	0.0248	0.3643
62068	0.0156	0.0234	0.2285	60294	0.0158	0.0246	0.2205

Table 15. Table 7: Local refinement with projection to the exact surface of a torus with a = 4 𝑎 4 a=4 , initial mesh-size of 0.1618.

$N_e$	$h$	${{ϵ}_{\mathrm{geom}}|}_{{𝒏}_{\mathrm{L2stab}}}$	Rate	${\Vert {𝒏}_f-{𝒏}_{\mathrm{L2stab}}\Vert }_{{\Sigma }_h}$	Effectivity Index, $E$
1278	0.0635	0.0393	-	1.5308	1.0055
3474	0.0307	0.0143	1.3900	0.8872	1.0147
8736	0.0151	0.0054	1.3608	0.5671	1.0024
21149	0.0095	0.0025	1.7165	0.3527	1.0013
53047	0.0056	0.0010	1.7790	0.2225	1.0009
131740	0.0028	0.0004	1.3741	0.1429	0.9999

Table 16. Table 8: Mean curvature errors for torus with a = 1 𝑎 1 a=1 .

$h$	${\gamma }_H$	$H_h$	$H_h$ rates	SSF	SSF rates	DLLB	DLLB rates
0.0527	0.05	0.1643	-	0.2183	-	0.4377	-
0.0245	0.05	0.0553	1.4178	0.0668	1.5412	0.3703	0.2179
0.0113	0.05	0.0250	1.0302	0.0404	0.6547	0.3126	0.2200
0.0057	0.05	0.0137	0.8746	0.0314	0.3670	0.2477	0.3394

Table 17. Table 9: Stabilization factor γ H subscript 𝛾 𝐻 \gamma_{H} as a function of mesh-size h ℎ h .

$h$	${\gamma }^{\ast }$	${ϵ}_H^{\ast }$	$\Delta {ϵ}_H$	relative change	$ϵ$ rate
0.0527	0.0576	0.1636	0.3715	2.2702	-
0.0245	0.1157	0.0454	0.4081	8.9971	1.6736
0.0113	0.2081	0.0124	0.3951	31.8512	1.6770

Equations

$$U_{\delta }(\Sigma )=\{\mathbf{x}\in {\mathbb{R}}^3:|\varphi (\mathbf{x})|<\delta \}$$

$${\Sigma }_h={\cup }_{K\in {\mathcal{K}}_h}K\subset U_{{\delta }_0}(\Sigma ).$$

$$v^e=v\circ \mathbf{p}\text{on Uδ0(Σ)}$$

$$(w^l{)}^e=w^l\circ \mathbf{p}=w\text{on Σh}$$

$$\mathbf{H}=-{\Delta }_{\Sigma }{\mathbf{x}}_{\Sigma }$$

$$\mathbf{H}=({\kappa }_1+{\kappa }_2)\mathbf{n}$$

$$(\mathbf{H},\mathbf{v}{)}_{\Sigma }=({\nabla }_{\Sigma }{\mathbf{x}}_{\Sigma },{\nabla }_{\Sigma }\mathbf{v}{)}_{\Sigma }\forall \mathbf{v}\in W$$

$$(\mathbf{v},\mathbf{w}{)}_{\omega }={\int }_{\omega }\mathbf{v}\cdot \mathbf{w}dx$$

$$\Vert \mathbf{v}{\Vert }_{\omega }^2={\int }_{\omega }\mathbf{v}\cdot \mathbf{v}dx.$$

$$({\mathbf{H}}_h,\mathbf{v}{)}_{{\Sigma }_h}+J_h({\mathbf{H}}_h,\mathbf{v})=({\nabla }_{{\Sigma }_h}{\mathbf{x}}_{{\Sigma }_h},{\nabla }_{{\Sigma }_h}\mathbf{v}{)}_{{\Sigma }_h}\forall \mathbf{v}\in W_h$$

$$J_h(u,v)$$

$$[{\mathbf{t}}_E\cdot {\nabla }_{{\Sigma }_h}u]={\mathbf{t}}_{E,K_1}\cdot {\nabla }_{{\Sigma }_h}{u}_1+{\mathbf{t}}_{E,K_2}\cdot {\nabla }_{{\Sigma }_h}{u}_2$$

$$u\approx \sum _i{\phi }_i{U}_i,$$

$${\nabla }_{{\Sigma }_h}u\approx \sum _i{\nabla }_{{\Sigma }_h}{\phi }_i{U}_i,$$

$${\nabla }_{{\Sigma }_h}{\phi }_i=:\left[\begin{array}{c}\frac{\partial {\phi }_i}{\partial x_{{\Sigma }_h}}\\ \frac{\partial {\phi }_i}{\partial y_{{\Sigma }_h}}\\ \frac{\partial {\phi }_i}{\partial z_{{\Sigma }_h}}\end{array}\right]={\mathbf{P}}_{{\Sigma }_h}\left[\begin{array}{c}\frac{\partial {\phi }_i}{\partial x}\\ \frac{\partial {\phi }_i}{\partial y}\\ \frac{\partial {\phi }_i}{\partial z}\end{array}\right].$$

$$(\mathbf{t}\cdot {\nabla }_{{\Sigma }_h}){\phi }_i=t_x\frac{\partial {\phi }_i}{\partial x_{{\Sigma }_h}}+t_y\frac{\partial {\phi }_i}{\partial y_{{\Sigma }_h}}+t_z\frac{\partial {\phi }_i}{\partial z_{{\Sigma }_h}}.$$

$$\mathbf{\Phi }:=\left[\begin{array}{ccccccc}{\phi }^1& 0& 0& {\phi }^2& 0& 0& \cdots \\ 0& {\phi }^1& 0& 0& {\phi }^2& 0& \cdots \\ 0& 0& {\phi }^1& 0& 0& {\phi }^2& \cdots \end{array}\right],$$

$${\mathbf{B}}_E:=\left[({\mathbf{t}}_1\cdot {\nabla }_{{\Sigma }_h})\mathbf{\Phi },({\mathbf{t}}_2\cdot {\nabla }_{{\Sigma }_h})\mathbf{\Phi }\right],$$

$$\mathbf{J}=\sum _{E\in \mathcal{E}}h{\int }_E{{\mathbf{B}}_E}^{⊺}{\mathbf{B}}_{E}dE.$$

$$\left(\mathbf{M}+{\gamma }_H\mathbf{J}\right)\mathbf{H}=\mathbf{S}\mathbf{x},$$

$$\mathbf{M}={\int }_{{\Sigma }_h}{\mathbf{\Phi }}^{⊺}\mathbf{\Phi }dx,$$

$$\mathbf{S}={\int }_{{\Sigma }_h}{\left({\nabla }_{{\Sigma }_h}\mathbf{\phi }\right)}^{⊺}{\nabla }_{{\Sigma }_h}\mathbf{\phi }dx={\int }_{{\Sigma }_h}{{\mathbf{B}}_S}^{⊺}{\mathbf{B}}_{S}dx,$$

$${\mathbf{B}}_S:=\left[\begin{array}{ccccccc}\frac{\partial {\phi }_1}{\partial x_{{\Sigma }_h}}& 0& 0& \frac{\partial {\phi }_2}{\partial x_{{\Sigma }_h}}& 0& 0& \cdots \\ \frac{\partial {\phi }_2}{\partial y_{{\Sigma }_h}}& 0& 0& \frac{\partial {\phi }_2}{\partial y_{{\Sigma }_h}}& 0& 0& \cdots \\ \frac{\partial {\phi }_1}{\partial z_{{\Sigma }_h}}& 0& 0& \frac{\partial {\phi }_2}{\partial z_{{\Sigma }_h}}& 0& 0& \cdots \\ 0& \frac{\partial {\phi }_1}{\partial x_{{\Sigma }_h}}& 0& 0& \frac{\partial {\phi }_2}{\partial x_{{\Sigma }_h}}& 0& \cdots \\ 0& \frac{\partial {\phi }_1}{\partial y_{{\Sigma }_h}}& 0& 0& \frac{\partial {\phi }_2}{\partial y_{{\Sigma }_h}}& 0& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right],$$

$$f(u,v)=\frac{1}{2}\left(a{u}^2+2buv+c{v}^2\right)$$

$${\int }_{A_M}\mathbf{K}(\mathbf{x})dA=\frac{1}{2}\sum _{j\in N}\left(\cot {\alpha }_{ij}+\cot {\beta }_{ij}\right)\left({\mathbf{x}}_i-{\mathbf{x}}_j\right),$$

$${\mathbf{K}}_h({\mathbf{x}}_i):=\frac{1}{2A_V}\sum _{j\in N}\left(\cot {\alpha }_{ij}+\cot {\beta }_{ij}\right)\left({\mathbf{x}}_i-{\mathbf{x}}_j\right).$$

$$H_h({\mathbf{x}}_i)=\frac{1}{2}|{\mathbf{K}}_h({\mathbf{x}}_i)|.$$

$$({\mathbf{n}}_h,\mathbf{v}{)}_{{\Sigma }_h}+J_h({\mathbf{n}}_h,\mathbf{v})=({\mathbf{n}}_K,\mathbf{v}{)}_{{\Sigma }_h}\forall \mathbf{v}\in W_h$$

$$\left(\mathbf{M}+{\gamma }_n\mathbf{J}\right){\mathbf{n}}_h=\mathbf{b},$$

$$\mathbf{b}:={\int }_{{\Sigma }_h}{\mathbf{\Phi }}^{⊺}{\mathbf{n}}_{K}d{\Sigma }_h.$$

Full text

Keywords: finite element method, discrete curvature, continuous interior penalty, projection method

1 Introduction

Our aim in this paper is to apply finite element techniques for computing geometrical quantities of interest in computer graphics applications, and to show that they can give accurate results, indeed more accurate that classical approaches. We restrict ourselves to closed surfaces approximated by piecewise linear simplices, and on such surfaces we consider three issues:

• •

accurate computation of the mean curvature vector;

• •

accurate computation of surface normals;

• •

adaptive refinement techniques for resolving the curvature.

We discretize the normal and curvature vectors using a piecewise linear finite element method based on tangential differential calculus, following the approach initiated by Dziuk [9]. This results in piecewise linear, continuous, vector fields on the discrete surface. In order to make comparisons with standard methods of computing curvature and normals, which are typically only represented at the vertices of the triangulated surface, we will focus mainly on the nodal values of the finite element fields.

Mean curvature. The mean curvature vector on a discrete surface plays an important role in computer graphics and computational geometry, as well as in certain surface evolution problems, see, e.g. [4, 5, 6, 8, 10, 11, 25]. It can be obtained by letting the Laplace–Beltrami operator act on the embedding of the surface in $\mathbb{R}^{3}$, and various formulas based on this fact have been suggested in the literature, see [21] and the references therein. It is known that the standard mean curvature vector based on the finite element discretization of the Laplace–Beltrami operator on a piecewise linear triangulated surface cannot be expected, in general, to give any order of convergence in the $L^{2}$ norm. More generally, for triangulated piecewise polynomial surfaces of order $k$ the expected convergence in $L^{2}$ norm is $k-1$, cf. [15, 7]. Convergence will also not occur in other standard discretization methods without restrictive assumptions on the mesh, see [29]. In this paper we employ a stabilized piecewise linear finite element method first suggested in [13] for approximation of the mean curvature vector, giving first order convergence in the $L^{2}$ norm for piecewise linear surfaces. The stabilization consists of adding suitably scaled terms involving the jumps in the tangent gradient of the discrete mean curvature vector in the direction of the outer co-normals at each edge in the surface mesh to the $L^{2}$–projection of the discrete Laplace–Beltrami operator used to compute the discrete mean curvature vector.

Normal vectors. Accurately determining the vertex normals on triangulated surfaces is of great importance in computer graphics for the computation of smooth shading [12, 24], and it is important in surface meshing/re-meshing [26, 23, 28] as well as smoothing (fairing) techniques [16]. We here extend the method suggested for computing the mean curvature vector, which can be seen as a general stabilization approach, to the problem of computing accurate vertex normals by stabilized $L^{2}$–projections.

Adaptive mesh improvement. Mesh improvement when the geometry is given by an analytical expression (or is otherwise known) can be obtained by local refinement of the simplices, putting new vertices on the known surface. The goal is then to resolve the curvature of the mesh in some predefined way. We suggest an approach based on the difference between the piecewise constant facet normals and the computed finite element normal field. This gives an estimate of the error in discrete facet normals which is closely related to the curvature of the geometry as will be discussed below. If the geometry is not a priori known but we are simply given a point cloud or a mesh, interpolation using vertex normals is standard, cf., e.g., Boschiroli et al. [3]. We combine one such approach, the PN triangle of Vlachos et al. [28], with our finite element normal fields and adaptive scheme in order to enhance the refined geometry.

The outline of the remainder of the paper is as follows: In Section 2 we introduce the discrete surface approximations, in Section 3 we define the stabilized mean curvature vector, in Section 4 we discuss a different schemes for computing vertex normals, including our stabilized projection method, in Section 5 we present an adaptive algorithm for resolving curvature, and in Section 6 we give some representative numerical results.

2 Meshed surfaces

Consider an embedded orientable closed surface $\mathbb{R}^{3}\supset\Sigma\in C^{2}$ with exterior unit normal $\bm{n}$. Let $\phi$ be the signed distance function such that $\nabla\phi=\bm{n}$ on $\Sigma$ and let $\bm{p}(\bm{x})=\bm{x}-\phi(\bm{x})\bm{n}(\bm{p}(\bm{x}))$ be the closest point mapping. Let $U_{\delta}(\Sigma)$ be the open tubular neighborhood

[TABLE]

for $\delta>0$ of $\Sigma$. Then there is $\delta_{0}>0$ such that the closest point mapping $\bm{p}(\bm{x})$ assigns precisely one point on $\Sigma$ to each $\bm{x}\in U_{\delta_{0}}(\Sigma)$.

We triangulate $\Sigma$ using a elementwise planar mesh $\mathcal{K}_{h}$ to obtain a quasiuniform triangulated surface

[TABLE]

Using the closest point mapping any function $v$ on $\Sigma$ can be extended to $U_{\delta_{0}}(\Sigma)$ using the pull back

[TABLE]

and the lifting $w^{l}$ of a function $w$ defined on $\Sigma_{h}$ to $\Sigma$ is defined as the push forward

[TABLE]

3 Approximation of the mean curvature vector

3.1 The continuous mean curvature vector

We define the tangential surface gradient $\nabla_{\Sigma}$ by $\nabla_{\Sigma}:=\bm{P}_{\Sigma}\nabla$, where $\nabla$ is the $\mathbb{R}^{3}$ gradient and $\bm{P}_{\Sigma}(\bm{x})=\bm{I}-\bm{n}(\bm{x})\otimes\bm{n}(\bm{x})$ is the projection onto the tangent plane $T_{\Sigma}(\bm{x})$ of $\Sigma$ at point $\bm{x}\in\Sigma$. The mean curvature vector $\bm{H}:\Sigma\rightarrow\mathbb{R}$ is then defined by

[TABLE]

where $\bm{x}_{\Sigma}:\Sigma\ni\bm{x}\mapsto\bm{x}\in\mathbb{R}^{3}$ is the coordinate map of $\Sigma$ into $\mathbb{R}^{3}$ and $\Delta_{\Sigma}=\nabla_{\Sigma}\cdot\nabla_{\Sigma}$ is the Laplace–Beltrami operator.

The relation between the mean curvature vector and mean curvature is given by the identity

[TABLE]

where $\kappa_{1}$ and $\kappa_{2}$ are the two principal curvatures and $(\kappa_{1}+\kappa_{2})/2=:H$ is the mean curvature, see [4].

The mean curvature vector satisfies the following weak problem: find $\bm{H}\in W=[H^{1}(\Sigma)]^{3}$ such that

[TABLE]

where $\nabla_{\Sigma}\bm{w}=\bm{w}\otimes\nabla_{\Sigma}$ for a vector valued function $\bm{w}$ and

[TABLE]

is the $L^{2}$–inner product on the set $\omega$ with associated norm

[TABLE]

Given the discrete coordinate map $\bm{x}_{\Sigma_{h}}:\Sigma_{h}\ni\bm{x}\mapsto\bm{x}\in\mathbb{R}^{3}$ and a discrete projection operator $\bm{P}_{\Sigma_{h}}=\bm{I}-\bm{n}_{h}\otimes\bm{n}_{n}$, where $\boldsymbol{n}_{h}$ denotes the piecewise constant facet normals, we define the stabilized discrete mean curvature vector $\bm{H}_{h}$ as follows. Let $V_{h}$ be the space of piecewise linear continuous functions defined on $\mathcal{K}_{h}$ and seek $\bm{H}_{h}\in W_{h}=[V_{h}]^{3}$ such that

[TABLE]

where $\nabla_{\Sigma_{h}}=\bm{P}_{\Sigma_{h}}\nabla$ and the stabilization term $J_{h}(\cdot,\cdot)$ is defined by

[TABLE]

Here $\gamma\geq 0$ is a stabilization parameter and $\mathcal{E}_{h}=\{E\}$ is the set of edges in the partition $\mathcal{K}_{h}$ of $\Sigma_{h}$. The jump of the tangential derivative in the direction of the outer co-normals at an edge $E\in\mathcal{E}_{h}$ shared by elements $K_{1}$ and $K_{2}$ in $\mathcal{K}_{h}$ is defined by

[TABLE]

where $u_{i}=u|_{K_{i}}$, $i=1,2,$ and $\bm{t}_{E,K_{i}}$ are the co-normals, i.e., the unit vectors orthogonal to $E$, tangent and exterior to $K_{i}$, $i=1,2$, see Figure 1. This stabilization method allows for proving first order convergence of the curvature vector $\|\bm{H}-\bm{H}_{h}^{l}\|_{\Sigma}\lesssim h,$ see [13].

3.2 Implementation issues

Using the standard Galerkin approximation,

[TABLE]

where $\varphi_{i}$ is the finite element basis functions and $U_{i}$ the nodal approximations of $u$ we have that

[TABLE]

were we define the tangential gradient of the basis function by

[TABLE]

The tangential derivative of the basis function is given by

[TABLE]

For vector–valued unknowns $\bm{u}$ we have $\bm{u}\approx\bm{\Phi}{\mathbf{u}}$ where ${\mathbf{u}}$ denotes nodal values and

[TABLE]

and using the notation $\bm{t}_{1}$ and $\bm{t}_{2}$ for the two co-normals on a given edge $E$, we define

[TABLE]

and the discrete stabilization matrix is given by

[TABLE]

The linear system corresponding to (6) becomes

[TABLE]

where ${\mathbf{M}}$ is the so called mass matrix, given by

[TABLE]

$\gamma_{H}$ is the mean curvature specific stabilization factor, ${\mathbf{S}}$ is the discrete Laplace-Beltrami operator defined by

[TABLE]

with

[TABLE]

${\mathbf{x}}$ is the coordinate vector of the nodal positions in the mesh, and ${\mathbf{H}}$ denotes the vector of vertex values of the approximate mean curvature vector.

3.3 Alternative approximations of the mean curvature vector

There exists several well known approaches to mean curvature estimation, for an extensive overview, see [18]. In the context of finite elements an alternative to ours is proposed by Heine in [14].

We shall compare our method to two types of such approaches:

1. fitting a surface locally to each vertex, see e.g., [1, Chap. 8.5] and 2) computing the discrete local Laplace-Beltrami operator, see e.g., [21, 8].

Smooth surface fit. Curvatures can be computed using a locally fitted quadratic function

[TABLE]

around a point ${\mathbf{x}_{i}}$ with $u$ and $v$ are local coordinates of the tangential plane to ${\mathbf{x}_{i}}$ such that $f(0,0)=0$. The tangential plane is determined using one of the edges connected to ${\mathbf{x}_{i}}$ and the normal at the same point. The idea is to compute the shape operator or Weingarten map of this function and subsequently the curvature. See [1, Chap. 8.5] for further details.

The discrete local Laplace-Beltrami operator. Let ${\mathbf{K}}$ denote the Laplace–Beltrami operator so that ${\mathbf{K}}(\bm{x})=2H(\bm{x})\bm{n}(\bm{x})$ at a given point $\bm{x}$ on the surface. On triangulated surfaces, one can use Gauss’ theorem to extract a discrete version of this operator in the nodes ${\mathbf{x}_{i}}$ of the mesh, cf. Meyer et al. [21]. The integral of ${\mathbf{K}}$ over the discrete 1-ring surface M on a triangulated surface is then given by

[TABLE]

where the angles $\alpha_{ij}$ and $\beta_{ij}$ are opposite to the edge $i\ j$ and $N$ is the set of neighbour vertices to ${\mathbf{x}_{i}}$, see Figure 2a. Given some definition $A_{V}$ of the local area surrounding a vertex ${\mathbf{x}_{i}}$ we can then define the discrete approximation ${\mathbf{K}}_{h}$ of ${\mathbf{K}}$ as

[TABLE]

In [21], it is proposed to use the Voroni regions as the definition for the local area, and an algorithm to improve the robustness for arbitrary meshes was provided. Similarly, Desbrun et. al [8] used the barycentric area to average the discrete Laplacian. In both cases, in order to compute the vertex normal, ${\mathbf{K}_{h}({\mathbf{x}_{i})}}$ is simply normalized and cases where the curvature is zero are treated by computing the mean face-normal of the 1-ring neighbourhood. The mean (discrete) curvature at the vertices is then given by

[TABLE]

4 Normal vector approximation

4.1 Stabilized projection of the normal field

In analogy with (6) we define the recovered discrete normal vector $\bm{n}_{h}$ as follows: find $\bm{n}_{h}\in W_{h}=[V_{h}]^{3}$ such that

[TABLE]

where $\bm{n}_{K}$ is the piecewise constant exterior normal to the facet elements $K$. The corresponding linear system becomes

[TABLE]

where $\gamma_{n}$ is the normal-specific stabilization factor, ${\mathbf{n}}_{h}$ the vector of vertex normals, and

[TABLE]

Note that (12) can be efficiently solved using a conjugate gradient method since ${\mathbf{M}}$ is symmetric, positive definite and sparse.

When translating the computed normal vector field to a set of discrete vertex normals, these will here be normalized (the nodal vectors contained in ${\mathbf{n}}_{h}$ are not in general of unit length).

4.2 Alternative approaches to computing vertex normals

Traditionally, vertex normals are estimated either from a local neighborhood of surrounding face normals using some type of local averaging, see e.g., [17, 23] and the references therein. Other estimation methodologies also exists such as local smooth surface fits, see, e. g., [20]. We use the notations for the local vertex normals introduced in [23] and give a brief description; see Figure 2b for an explanation of the notations used.

Mean weighted equally. Arguably, the most widespread estimation of the vertex normal was introduced by Gouraud [12] as

[TABLE]

where $\bm{n}_{i}$ is the face-normal of triangle $i$ , $n$ is the total number of triangles that share a common vertex for which the vertex normal is to be estimated and $|.|$ denotes the norm. Note that we shall subsequently omit making the normalization step of the vertex normal explicit and assume $\bm{n}=\hat{\bm{n}}:=\bm{n}/|\bm{n}|$.

Mean weighted by angle. A vertex normal approximation using angles between the inner edges was proposed by Thürrner and Wüthrich [27].

[TABLE]

where $\alpha_{i}$ is the angle between two edges $\bm{e}_{k}$ and $\bm{e}_{k+1}$ of a face $i$ sharing the vertex.

Mean weighted by sine and edge length reciprocals. Max [19] proposed several methods of weighting the face normals, one of which is to weight by the sine and edge length reciprocals to take into account the difference in lengths of surrounding edges.

[TABLE]

Mean weighted by areas of adjacent triangles. Another method proposed by Max [19] is to weight the normals by the area of the face.

[TABLE]

where the symbol $\times$ denotes the vector cross product.

Mean weighted by edge length reciprocals. Max [19] also proposed to just use the edge length reciprocals as weights.

[TABLE]

Mean weighted by square root of edge length reciprocals. Finally, Max [19] also suggested to use the square root of the length reciprocals.

[TABLE]

Normal from the discretized local Laplace-Beltrami operator. Another approach is to define the normal using the discretized local Laplace-Beltrami operator (DLLB) defined in Section 3.3. The normal is defined by normalizing the discrete mean curvature vector.

[TABLE]

In the numerical example below, Section 6, we compare the accuracy of these different approaches.

5 Adaptive algorithm

5.1 Error estimate

We base our adaptive algorithm on the Zienkiewicz–Zhu approach [30] which employs the difference between recovered derivatives and actual discrete piecewise derivatives of a finite element solution. By analogy we consider the piecewise constant normals to play the role of the piecewise derivatives, and compare these to the $L^{2}-$projected normals.

Since we are focusing on vertex normals, and since we will in the following compare methods that only produce such normals, we define a norm which is an approximation of the $L^{2}$–norm,

[TABLE]

where $\text{meas}(K)$ denotes the area of $K$ and $\boldsymbol{x}^{i}_{K}$ the vertex coordinates on $K$. This represents a Newton–Cotes numerical integration scheme for the $L^{2}(\Sigma_{h})$–norm using the vertices as integration points. The error in normals is thus approximated

[TABLE]

and we aim at achieving

[TABLE]

where TOL is a given tolerance. We note that we also have

[TABLE]

where $h$ is the local mesh size and $\nabla_{\Sigma}\bm{n}$ is the curvature tensor, which indicates that we counter large curvature by reduced mesh size for resolution of the geometry.

5.2 Triangle refinement

In cases where the exact geometry is not accessible, we consider triangle refinement approaches that utilise vertex normals for interpolation. An overview of such methods is given by Boschiroli et al. in [3]. Nagata [22] proposed a simple quadratic interpolation of triangles using vertex normals and positions at the end-nodes. The approach by Nagata depends on a curvature parameter that fixes a curvature coefficient in order to stabilize the method. The curvature coefficient is highly dependent on the vertex normal, and in cases where normals are near parallel, the method cannot capture inflections and without a stabilizing parameter, cusps will be introduced to the surface, see [23] where the authors point out this problem and suggest a possible solution. The solution suggested in [23] eliminates the problem of cusps in the interpolated surface but also eliminates the inflection, since the segment becomes linear. Another approach is to use higher order interpolation which are able to capture inflection points.

5.2.1 PN triangles

Vlachos et al. [28] proposed a cubic interpolation scheme that similarly to Nagata only depends on the positions and vertex normals of a triangular patch. We here write their algorithm in a vectorized manner. Let then $b:\,\mathbb{R}^{2}\rightarrow\mathbb{\mathbb{R}}^{3}$ denote a cubic triangular patch given by

[TABLE]

Here ${\mathbf{U}}$ is the matrix representation of the parameters defined by

[TABLE]

where $u=i/N$, $v=j/N$ for $i,j=\{0,1,\ldots,N\}$ such that $w:=1-u-v\geq 0$. Here $N$ gives a subtriangulation of the initial patch, see Figure 3. ${\mathbf{B}}$ denotes the cubic coefficients in matrix form and is given by

[TABLE]

where $\bm{b}$ denote the control points of the control grid for the PN triangle, see Figure (4), and are defined as follows:

[TABLE]

where $\bm{p}_{i}$ and $\bm{n}_{i}$ are the input corner points and normals. Finally the total set of interpolated points is given as a matrix product by

[TABLE]

Note that ${\mathbf{U}}$ can be evaluated for a certain number of refinements $N$ in a pre-processing step. In the local refinement section of this paper we use $N=1$ see Figure 5. As for the internal vertex normal computation, we do not interpolate the normals locally, instead we compute $\bm{n}_{h}$ using (12) for the total mesh in each iteration. The reason behind why we limit the tessellation step to 1 is the subsequent complexity of the local refinement procedure.

5.3 Local refinement procedure

Since the PN refinement with $N=1$ splits the face of a flat triangle into four child elements, we need a way of handling the hanging nodes. In this work we adapt the Red-Green refinement method proposed by Banks et al in [2]. This method preserves the aspect ratio of the initial mesh which is crucial in order to secure the accuracy of the associated finite element method.

6 Numerical examples

6.1 Geometry

We choose to analyze the errors on an implicitly defined torus which we can modify in order to generate slightly more complex features. The surface equation for the torus is given by

[TABLE]

where $R$ is the torus radius, $r$ the tube radius and $a$ is a “squish-factor” used to squish the torus in the z-direction in order to induce a higher curvature on the inside and outside, see Figure 6. In the following, the torus will be analyzed with $a=1$ and $a=4$, in order to compare errors with respect to strongly and smoothly varying curvature.

6.2 Vertex normal error

What follows is a comparison of different vertex normals with the exact normal. The measure for the mesh-size used in this context is defined as

[TABLE]

where $N_{\text{v}}$ denotes the number of vertices in the mesh.

Using an implicitly defined surface $\Sigma=\{\bm{x}:\phi(\bm{x})=0\}$, where $\phi$ is a signed distance function with the property $|\nabla\phi|=1,$ we have that $\bm{n}(\bm{x}_{\Sigma})=\nabla\phi(\bm{x}_{\Sigma})$. As discussed above, we will use (20) and define the error as

[TABLE]

where $\bm{n}_{a}$ is the approximate and $\bm{n}_{e}$ the exact normal defined by $\bm{n}_{e}=\nabla\phi$, computed at the vertex $i$ using $\bm{n}_{e}^{i}=\nabla\phi(\bm{x}^{i})$. The convergence rates are defined as

[TABLE]

6.3 Evaluation of the accuracy of computed vertex normals

The vertex normal error analysis was done on an unstructured mesh of a torus with $R=1$, and $r=1/2$ and $a=\{1,4\}$, see Figure 6.

The convergence of $L^{2}$ errors $\epsilon$ defined in (34) are shown in Figure 7 were it can be seen that the stabilized $L^{2}$–projection of the normals converges optimally. The raw data for this graph is available in Table 1. The relative difference between the stabilized $L^{2}$ normals and the next best traditional method $\bm{n}_{\mathrm{MWA}}$ can be seen in Table 2 where we can see a relative error decrease from $\epsilon_{\text{MWA}}$ of $\sim 29\%$ to $\sim 88\%$ depending on mesh-size and geometry. The convergence rates can be viewed in Table 3. In the next section we shall analyze the impact of the stabilization on the normal errors.

6.4 Effect of the stabilization on the accuracy of the computed normal

In this section we analyze the influence of the stabilization factor on the vertex normal error numerically by employing a golden search method to find the optimal stabilization factor $\gamma_{n}^{*}$ that minimizes the normal error $\epsilon$ defined in (34).

[TABLE]

where we use $\gamma_{n}^{0}=0$ and $\gamma_{n}^{1}=1$. This is done for several mesh-sizes and on a torus with $a=1$ and $a=4$, see Figure 8 and 9. Notice how the curves become more planar, i.e., choosing a “good” $\gamma_{n}$ becomes less sensitive with the decrease in $h$.

The error difference is shown in Figure 10 and Table 4 where $\epsilon_{L^{2}}$ is the $L^{2}$ error, defined in (34), of the $L^{2}$ vertex normals without stabilization and $\epsilon_{L_{\mathrm{stab}}^{2}}$ is the error of the stabilized $L^{2}$ vertex normal which is stabilized with a optimal stabilization factor $\gamma^{*}$.

6.5 Interpolation

In a 2D case we can see in Figure 11 how the choice of vertex normals affects the resulting cubic interpolation. The initial mesh is coarse and the (unstabilized) $L^{2}$–projected normals are not just depending on the nearest neighbors to each vertex but globally. The resulting difference is apparent.

We compare the impact of different vertex normals on the interpolation by measuring the geometrical error defined as

[TABLE]

where $\bm{x}_{\Sigma}(\bm{n})$ denotes the discrete surface interpolated with a particular normal approximation method. We measure the $L^{2}_{h}$–norm of the signed distance. The refinement algorithm employed is the PN triangles using 1 tessellation per face see Figure 5. The mesh-size in this section is defined as

[TABLE]

where $N_{\text{e}}$ denotes the number of elements and $A_{K}$ is the area of the $K$-th element. The initial mesh-size is $h=0.1618$ and the initial $L^{2}_{h}$–norm of the signed distance error is $\epsilon_{\mathrm{geom}}=0.0863$. See Figure 12 for the convergence comparison, Table 5 for the regular refinement data and Table 6 for the local refinement data.

Examples of interpolation using PN triangles with local refinement are shown in Figure 13 for the Torus, Figure 14 for the Utah teapot and Figure 15 for the Stanford bunny.

The local refinement method is compared to local refinement with projection to the exact surface, see Figure 16. We compare the approximite normal error with the exact normal error by computing the effectivity index, given by

[TABLE]

where $\bm{n}_{e}$ is the exact normal to the surface, $\bm{n}_{f}$ is the face normal and $\bm{n}_{L^{2}_{\mathrm{stab}}}$ is the recovered stabilized $L^{2}$–projected normal, see Table 7.

6.6 Mean curvature

The mean curvature is computed on a structured and unstructured torus with $R=1$, $r=1/2$ and $a=1$. We compare the mean curvature approximation to the exact mean curvature, the smooth surface fit approach (SSF) and the discrete local Laplace-Beltrami (DLLB) approach described in Section 3.3, and our stabilized discrete curvature vector solving (10). In the last case we compute the mean curvature $H_{h}$ through

[TABLE]

where $\bm{n}_{h}$ denotes the normal computed using the stabilized $L^{2}$–projection from (12). In our computational experience, this gives a more accurate result than the immediate $H_{h}=\frac{1}{2}|\bm{H}_{h}|$.

In Figure 17 we give iso-plots of the mean curvature. Figure 18 shows the convergence of mean curvature.

6.6.1 Stabilization sensitivity

We begin by analyzing how sensitive the mean curvature approximation is with respect to the stabilization factor $\gamma_{H}$. The mean curvature error is defined by

[TABLE]

where $H_{\mathrm{exact}}$ is the exact mean curvature computed on a torus with $R=1,$ $r=\frac{1}{2}$ and $a=1$. A golden search method is used to find the optimal stabilization factor and subsequently to numerically analyze the impact of the stabilization choice with regards to the error and mesh-size.

[TABLE]

where $\gamma_{H}^{0}=0$ and $\gamma_{H}^{1}=0.15$. The result of this optimization is a discrete function of mean curvature error with respect to the stabilization factor, $\epsilon_{H}(\gamma_{H})$, see Figure 19. In Table 9 we present the optimal stabilization factors and differences in mean curvature error defined as $\Delta\epsilon_{H}=\epsilon_{H}(\gamma_{H}^{0})-\epsilon_{H}(\gamma_{H}^{*})$, where $\gamma_{H}^{0}=0$ and $\gamma_{H}^{*}$ is the value of gamma that minimizes $\epsilon_{H}$. Notice how the curves become more planar, i.e., choosing a $\gamma_{H}$ that improves the solution becomes less sensitive with the decrease in $h$.

Acknowledgements

This research was supported in part by the Swedish Foundation for Strategic Research Grant No. AM13-0029, the Swedish Research Council Grants Nos. 2011-4992, 2013-4708, and the Swedish Research Programme Essence.