Finite Element Methods for Geometric Modeling ... - Semantic Scholar
Finite Element Methods for Geometric Modeling and Processing Using General Fourth Order Geometric Flows Guoliang Xu State Key Laboratory of Scientific and Engineering Computing Institute of Computational Mathematics, Academy of Mathematics and System Sciences, Chinese Academy of Sciences, Beijing 100080, China [email protected]
Abstract. A variational formulation of a general form fourth order geometric partial differential equation is derived, and based on which a mixed finite element method is developed. Several surface modeling problems, including surface blending, hole filling and surface mesh refinement with the G1 continuity, are taken into account. The used geometric partial differential equation is universal, containing several well-known geometric partial differential equations as its special cases. The proposed method is general which can be used to construct surfaces for geometric design as well as simulate the behaviors of various geometric PDEs. Experimental results show that it is simple, efficient and gives very desirable results. Keywords: Geometric PDE, Surface blending, hole filling, Surface mesh refinement, Mixed finite element method.
1
Introduction
In this paper, we use a general form geometric partial differential equation (GPDE) to solve several surface modeling problems. The GPDE we use includes many well known equations as its special cases. It is nonlinear and geometrically intrinsic, which is solved using a mixed finite element method. The problems we solve include surface blending, hole filling and mesh refinement with the G1 continuity. For the problems of surface blending and hole filling, we are given triangular surface meshes of the surrounding area. Triangular surface patches need to be constructed to fill the openings enclosed by the surrounding surface mesh and interpolate its boundary with the G1 continuity (see Fig. 2 and 3). For mesh refinement problem, we are given an initial surface mesh. We construct a sequence of surface meshes which interpolates the vertices of the previous meshes with the specified G1 continuity (see Fig. 4). There are basically two approaches for solving a GPDE numerically on surface. One is the generalized finite divided difference method (see [25,28]), the
Project support in part by NSFC grant 60773165 and National Key Basic Research Project of China (2004CB318000).
F. Chen and B. J¨ uttler (Eds.): GMP 2008, LNCS 4975, pp. 164–177, 2008. c Springer-Verlag Berlin Heidelberg 2008
Finite Element Methods for Geometric Modeling
165
other is the finite element method (see [1,10,16,18,26]). The approach we adopt in this paper is a mixed finite element method consisting of using the Loop’s subdivision basis for representing the surfaces and the piecewise linear basis for representing the curvatures. Loop’s basis have been used in finite element analysis and design for thin-shell structures (see [6,7]). The generalized finite divided difference method has been used in [25,28]. It is well known that the finite divided difference method is simpler and easier to implement, but it lacks the convergence analysis. The finite element method is not as simple as the finite divided difference method, but is based on a well developed mathematical foundation. Previous Work. The well-known biharmonic equation and its various variations have been repeatedly used to solve the problems of geometry modeling, such as blending and the hole filling (see [3,4]), interactive surface design (see [15], [24]) and interactive sculpting (see [14]). A PDE method based on the biharmonic equation is proposed [19] with certain functional constraints, such as geometric constraints, physical criteria and engineering requirements, which can be incorporated into the shape design process. The main advantageous of using the biharmonic equation is that the equation is linear, hence it is easy to solve. However, the biharmonic equation is not geometrically intrinsic and its solution depends on the concrete surface parametrization used. Furthermore, the biharmonic equation is defined on a rectangular domain in general, hence it is inappropriate for modeling surfaces with arbitrary shaped boundaries. In contrast, the equations we use in this paper are geometrically intrinsic without the limitation of the shape of the surface boundaries. These equations are called geometric partial differential equations (GPDE). The simplest GPDE may be the mean curvature flow (MCF). MCF and its variations have been intensively used to smooth or fair noisy surfaces (see [1,9,11] for references). These second order equations have been shown to be the most important and efficient flows for fairing or denoising surface meshes. However, for solving the surface modeling and design problems, MCF cannot achieve the G1 continuity at the boundary. To achieve smooth joining of different surface patches in solving the surface blending problem (see [8,21,22,25]), the fourth order geometric flows, such as Willmore flow and surface diffusion flow, have been employed. For instance, the surface diffusion flow is used in [21] to fair surface meshes with the G1 conditions in a special case where the meshes are assumed to have the subdivision connectivity. This equation is also used in [22] for smoothing meshes while satisfying the G1 boundary conditions and the requirement of the outer fairness (the smoothness in the classical sense) and the inner fairness (the regularity of the vertex distribution). Willmore flow has been used for the surface restoration (see [8]) and surface modeling (see [26]) using the finite element method. Three fourth order flows (surface diffusion flow, Willmore flow and quasi surface flow) are treated separately using finite element method in [26]. Main Results. We derive a weak form formulation for a general form fourth order geometric flow. Based on this weak form a mixed form formulation is
166
G. Xu
proposed. The equation is solved by a newly proposed mixed finite element method over surfaces. The proposed approach is simple, efficient and general as well. It handles a class of geometric PDEs in the same way and solves several surface modeling problems in the same manner. The approach is applied to solving each of the surface modeling problems and gives very desirable results for a range of the complicated surface models. The rest of the paper is organized as follows: Section 2 introduces some background material from differential geometry and our earlier publications. Section 3 describes the used nonlinear GPDE and its variational form. In section 4, we give the details of the discretization of the variational form in both the spatial and the temporal directions. Comparative examples to illustrate the different effects achievable from the solutions of these GPDEs are given in section 5.
2
Differential Geometry Preliminaries
In this section, we introduce the used notations, curvatures and several geometric differential operators, including Laplace-Beltrami operator and GiaquintaHildebrandt operator etc. Some results used in this paper are also presented. Let S = {x(u, v) ∈ R3 ; (u, v) ∈ Ω ⊂ R2 } be a regular and sufficiently smooth parametric surface. To simplify the notation we sometimes write w = (u, v) = (u1 , u2 ). Let gαβ = xuα , xuβ and bαβ = n, xuα uβ be the coefficients of the first and the second fundamental forms of S with ∂x ∂2x , xuα uβ = , α, β = 1, 2, n = (xu × xv )/xu × xv , α ∂u ∂uα ∂uβ where ·, · and × denote the usual scalar and cross products of two vectors, respectively, in Euclidean space R3 . Set xuα =
g = det[ gαβ ],
[ g αβ ] = [ gαβ ]−1 , b = det[ bαβ ],
[ bαβ ] = [ bαβ ]−1 .
Curvatures. To introduce the notions of mean curvature and Gaussian curvature, let us first introduce the concept of Weingarten map or shape operator. The shape operator of surface S is a self-adjoint linear map on the tangent space Tx S := span{xu , xv } defined by W : Tx S → Tx S, such that W(xu ) = −nu , W(xv ) = −nv .
(1)
We can easily represent this linear map by a 2 × 2 matrix S = [ bαβ ][ g αβ ]. In particular, [nu , nv ] = −[xu , xv ]S T is valid. The eigenvalues k1 and k2 of S are principal curvatures of S and their arithmetic average and product are the mean curvature and the Gaussian curvature, respectively. That is H=
tr(S) k1 + k2 = , 2 2
K = k1 k2 = det(S).
Tangential gradient operator. Let f be a C 1 smooth function on S. Then the tangential gradient operator ∇ acting on f is given by (see [12], page 102) ∇f = [xu , xv ][ g αβ ][fu , fv ]T ∈ R3 ,
(2)
Finite Element Methods for Geometric Modeling
167
where fu and fv denote the first order partial derivatives of f with respect to arguments u and v. For a vector-valued function f = [f1 , · · · , fk ]T ∈ C 1 (S)k , we define its gradient by ∇f = [∇f1 , · · · , ∇fk ] ∈ R3×k . It is easy to see that ∇x = [xu , xv ][g αβ ][xu , xv ]T , ∇n = −[xu , xv ][g αβ ]S[xu , xv ]T ,
(3)
and both ∇x and ∇n are symmetric 3 × 3 matrices. Second tangent operator. Let f be a C 1 smooth function on S. Then we define the second tangent operator 3 acting on f by (see [27]) 1 b22 −b12 . (4) 3f = [xu , xv ][hαβ ][fu , fv ]T ∈ R3 , with [hαβ ] = g −b12 b11 Divergence operator. Let v be a C 1 smooth vector field on S. Then the divergence of v is defined by ∂ ∂ √ 1 , g [ g αβ ] [xu , xv ]T v . (5) div(v) = √ g ∂u ∂v For a matrix-valued function V = [v1 , · · · , vk ] ∈ C 1 (S)3×k , we define div(V) = [div(v1 ), · · · , div(vk )]T ∈ Rk . Laplace-Beltrami operator. Let f be a C 2 smooth function on S. Then ∇f is a smooth vector field on S. The Laplace-Beltrami operator (LBO) Δ applying to f is defined by (see [13], page 83) Δf = div(∇f ). Giaquinta-Hildebrandt operator. Let f be a C 2 smooth function on S. Then the Giaquinta-Hildebrandt operator acting on f is given by (see [17], page 84) 2f = div(3f ). For a vector-valued function f = [f1 , · · · , fk ]T ∈ C 1 (S)k , we define Δf = div(∇f ),
2f = div(3f ).
(6)
For the operators introduced above, we can prove the following theorems that are used in the sequel. Detailed proofs can be found in our technical report [29]. Theorem 1. Let x ∈ S. Then Δx = 2Hn, Δn = −2∇H − 2HΔx + 2x.
(7)
Theorem 2. Let x ∈ S. Then 2x = 2Kn, 2n = −∇K − KΔx = −∇K − H2x.
(8)
To derive the weak form of the considered geometric flows, we introduce two Green’s formulas.
168
G. Xu
Theorem 3 (Green’s formula I for LB operator). ([5], page 142) Let f ∈ C 2 (S), h ∈ C 1 (S), with at least one of them compactly supported. Then hΔf + ∇f, ∇h dA = 0. M
Theorem 4 (Green’s formula I for GH operator). Let f ∈ C 2 (S), h ∈ C 1 (S), with at least one of them compactly supported. Then hf + ∇f, 3h dA = 0. M
3
Weak and Mixed Forms of a General Fourth Order GPDE
Euler-Lagrange equation. Let f (H, K) ∈ C 1 (R × R) be a given Lagrange function. Define energy functional f (K, H)dA. (9) E (S) = S
In this paper we use a vector-valued Euler-Lagrange equation derived from the complete-variation x(w, ε) = x(w) + εΦ(w), w ∈ Ω, |ε| 1, Φ ∈ Cc∞ (Ω)3 ,
(10)
for the functional E (S). We refer a result from [28] as the following theorem: Theorem 5. Let f (H, K) ∈ C 1 (R × R). Then the Euler-Lagrange equation of the integral E (S) from the complete-variation (10) is 1 2(fK n) + Δ(fH n) − div(fH ∇n) − div[(f − 2KfK )∇x]) = 0 ∈ R3 . (11) 2 Theorem 6. The weak form of (11) can be written as
1 fH ∇x3φ + (f − 2HfH − 2KfK )∇x∇φ + fH nΔφ + fK n2φ dA = 0, 2 S (12) ∀φ ∈ Cc∞ (Ω). Proof. By direct calculation, it is not difficult to derive that [g αβ ]S = 2H[g αβ ] − K[bαβ ]. Using this equality and (3), we have ∇n = 3x − 2H∇x.
(13)
Multiplying both sides of (11) with φ ∈ Cc∞ (Ω), then taking integral over S, using Green formulas, and finally using (13), we have
1 2(fK n) + Δ(fH n) − div(fH ∇n) − div[(f − 2KfK )∇x]) φ dA 0= 2 S
1 fK n2φ + fH nΔφ + [fH (3x−2H∇x)]∇φ+[(f −2KfK )∇x]∇φ) dA. = 2 S Noticing that 3x∇φ = ∇x3φ, we obtain (12) in the end.
Finite Element Methods for Geometric Modeling
169
A General form Geometric Flow. For a given function f (H, K), let S0 be a given initial surface with boundary Γ . Then the geometric flow consists of finding a family {S(t); t ≥ 0} of smooth orientable surfaces in R3 which evolve according to the following equation ∂x 1 ∂t + 2(fK n) + 2 Δ(fH n) − div(fH ∇n) − div (f − 2KfK )∇x = 0, (14) S(0) = S0 , ∂S(t) = Γ. Weak form formulation. Using (12), we can obtain the weak form of equation (14) as follows: ⎧ ∂x ⎪ ⎪ fH ∇x3φ + (f − 2HfH − 2KfK )∇x∇φ φ dA + ⎪ ⎨ S ∂t S
+ 12 fH nΔφ + fK n2φ dA = 0, ∀φ ∈ Cc∞ (Ω), ⎪ ⎪ ⎪ ⎩ S(0) = S0 , ∂S(t) = Γ.
(15)
Assumptions. In order to make the generated flow meaningful and well-defined. We first require that f is a smooth function about its arguments. We also assume that f is an algebraic function, meaning it does not involve differential and integral operations about its variables. Suppose fH and fK could be represented as fH = 2αH + 2βK + μ, fK = 2γH + 2δK + ν. (16) Let h = 2Hn, k = 2Kn. Then using (15) and hψ dA = Δxψ dA = − ∇x∇ψ dA, ∀ψ ∈ Cc∞ (Ω), S S S kϕ dA = 2xϕ dA = − ∇x3ϕ dA, ∀ϕ ∈ Cc∞ (Ω), S
S
S
we obtain a mixed form of equation (14) as follows: Find [x, h, k] ∈ R3×3 , such that ⎧ ∂x ⎪ ⎪ φ dA + fH ∇x3φ + (f − 2HfH − 2KfK )∇x∇φ ⎪ ⎪ ∂t ⎪ S S ⎪ ⎪ 1 ⎪ ⎪ + (αh + βk + νn)Δφ+ (γh + δk + μn)2φ dA = 0, ∀φ ∈ Cc∞ (Ω), ⎨ 2 ⎪ hψ dA + ∇x∇ψ dA = 0, ∀ψ ∈ Cc∞ (Ω), ⎪ ⎪ ⎪ S S ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ kϕ dA + ∇x3ϕ dA = 0, ∀ϕ ∈ Cc∞ (Ω). S
S
(17) From the definitions of tangential operators, it is easy to derive that ∇x3φ = 3φ, ∇x∇φ = ∇φ. Further notice that 2H = nT h, 2K = nT k, the weak form (17) can be rewritten as: Find [x, h, k] ∈ R3×3 , such that
170
G. Xu
⎧ ∂x ⎪ ⎪ μ∇x3φ + η∇x∇φ + 3φnT (αh + βk) φ dA + ⎪ ⎪ ⎪ S ∂t S ⎪ ⎪ ⎪ +∇φnT [(ω − fH )h + ( − fK )k] ⎪ ⎪
⎪ 1 ⎨ + (αh + βk + νn)Δφ+ (γh + δk + μn)2φ dA = 0, 2 ⎪ ⎪ ⎪ hψ dA + ∇x∇ψ dA = 0, ∀ψ ∈ Cc∞ (Ω), ⎪ ⎪ ⎪ ⎪ S S ⎪ ⎪ ⎪ ⎪ ⎩ kϕ dA + ∇x3ϕ dA = 0, ∀ϕ ∈ Cc∞ (Ω), S
∀φ ∈ Cc∞ (Ω),
S
(18) where ω, and η are defined by f = 2ωH + 2K + η. System (18) is the starting point of the finite element discretization.
4
Solving the GPDE by Mixed Finite Element Methods
Consider a triangular surface mesh Ω with vertices x1 , x2 , · · · , xn . For each vertex xi , we associate it with a basis function φi . Then the surface S is approximately represented as x=
n
φj (q)xj ∈ R3 , q ∈ Ω,
(19)
j=1
where x1 , x2 , · · · , xn are regarded as the control vertices of S. Now we classify these control vertices into several categories. The first category consists of interior vertices, denoted as x1 , x2 , · · · , xn0 . The positions of these vertices are subject to change (unknown). For the problems of hole filling and blending, the interior vertices are those in the openings (see Fig. 1(a), the solid dots). For the problem of mesh refinement, all the non-interpolating vertices are interior (see Fig. 1(b), the solid dots). Apart from the interior vertices, the remaining vertices are classified as 1-ring neighbors (see Fig. 1(a), the empty
Ω0
Ω1
1
1
(a)
(b)
Fig. 1. Classification of vertices. (a) The solid dots are interior control vertices. The empty dots are 1-ring neighbor (boundary) vertices. The solid squares are 2-ring neighbor vertices. The shaded region is Ω0 . The region enclosed by the 2-ring neighbor vertices is Ω1 . (b) The solid dots are the interior control vertices. The empty ones are the interpolated vertices.
Finite Element Methods for Geometric Modeling
171
dots), 2-ring neighbors (see Fig. 1(a), the solid squares), · · · , according to their distance to the interior vertices. For instance, the 1-ring neighbors are all the vertices adjacent to the interior ones, the 2-ring neighbors are all the remaining vertices adjacent to the 1-ring neighbors, · · · . The k-ring neighbor control vertices are denoted as xnk−1 +1 , · · · , xnk . All the 1-ring, 2-ring, · · · , neighbor vertices are fixed. In our mixed form formulations of GPDE, the mean and Gaussian curvature normals are also treated as unknown vector-valued functions. Hence, for each control vertex xi , we also associate it with a continuous basis function ψi for representing the curvatures. For instance, mean and Gaussian curvature normals are represented as h(x) =
n
ψj (q)hj ∈ R3 ,
k(x) =
j=1
n
ψj (q)kj ∈ R3 , q ∈ Ω,
(20)
j=1
respectively. Here hj and kj denote the mean and Gaussian curvature normals at xj , respectively. Suppose that φi (q) and ψi (q) have compact supports and the sizes of their supports are 2 and 1, respectively. At a k-ring neighbor vertex, the surface point does not depend upon the interior vertices if k ≥ 2, so the mean and Gaussian curvature normals there are well defined from the surface. However, at the interior and the 1-ring neighbor vertices, the corresponding surface points relate to unknown vertices. Then the mean and Gaussian curvatures there cannot be determined from the surface. Hence we treat these curvatures as unknowns. In a word, the unknown vertices are x1 , x2 , · · · , xn0 , the unknown mean curvature normals are h1 , h2 , · · · , hn0 , hn0 +1 , · · · , hn1 , and the unknown Gaussian curvature normals are k1 , k2 , · · · , kn0 , kn0 +1 , · · · , kn1 . In the next two sub-sections, we will obtain matrix forms of the used GPDE by substituting (19)–(20) into (18). Note that the proposed scheme needs the values of the mean and Gaussian curvature normals at each vertex. For an interior or 1-ring neighbor vertex, the mean and Gaussian curvature normals are unknown in the equation to be solved, but initial values at t = 0 is required. For a k-ring neighbor vertex (k ≥ 2), these vectors need to be computed previously and keep fixed afterwards (see Remark 4.3 for how these vectors are computed). A default choice of these vectors is to set them as zeros. 4.1
Spatial Direction Discretizations
We construct one equation for each unknown vertex xj (j = 1, · · · , n0 ) and other two equations for the unknown vector-valued curvatures hj and kj (j = 1, · · · , n1 ). This is achieved by the following operations: (i) Substituting (19)–(20) into (18). (ii) Taking the test functions φ = φi (i = 1, · · · , n0 ) and
172
G. Xu ∂x (t)
j ψ = ϕ = ψi (i = 1, · · · , n1 ). (iii) Using the fact that ∂t = 0 for j > n0 (since xj is fixed). We then obtain the following matrix form of (18): ⎧ ∂Xn0 (t) (p) (m) (g) ⎪ + Ln3 Xn3 (t) + Ln2 Hn2 (t) + Ln2 Kn2 (t) + Bn0 = 0, ⎪ ⎨ Mn0 ∂t
⎪ ⎪ ⎩
(1)
(1)
Mn2 Hn2 (t) + Ln3 Xn3 (t) = 0,
(1) Mn2 Kn2 (t)
where
+
(2) Ln3 Xn3 (t)
(21)
= 0,
Xj (t) = [x1 (t), · · · , xj (t)]T ∈ Rj×3 , Hj (t) = [h1 (t), · · · , hj (t)]T ∈ Rj×3 , Kj (t) = [k1 (t), · · · , kj (t)]T ∈ Rj×3 ,
are unknowns and n ,K
0 MK = (mij )ij=1 , n0 ,K (p) (p) , LK = lij nij=1 1 ,K (1) (1) , LK = lij
ij=1
n1 ,K (1) (1) MK = mij , nij=1 n0 ,K 0 ,K (m) (m) (g) (g) LK = lij , LK = lij , ij=1 ij=1 n1 ,K n0 (2) (2) LK = lij , Bn0 = bTi i=1 ij=1
are the coefficient matrices. The elements of these matrices are defined as follows: (1) mij = φi φj dA, mij = ψi ψj dA, S S (p) lij = μ(∇ψj )T 3φi + η(∇ψj )T ∇φi dA ∈ R, S
1 (m) ψj [α3φi + (ω − fH )∇φi ]nT + ( αψj Δφi + γψj 2φi )I3 dA ∈ R3×3 , lij = 2 S
1 (g) ψj [β3φi + ( − fK )∇φi ]nT + ( βψj Δφi + δψj 2φi )I3 dA ∈ R3×3 . lij = 2 S (1) (2) T lij = (∇φj ) ∇ψi dA , lij = (∇φj )T 3ψi dA, S S 1 bi = [n( νΔφi + μ2φi )]dA. 2 S The integrations in these matrix elements are computed by a 12-point Gaussian quadrature rule. That is, each domain triangle is subdivided into four subtriangles and a three-point Gaussian quadrature rule is employed on each of the sub-triangles (see [2] for a set of Gaussian quadrature rules). In system (21), moving the terms relating to the known vertices xn0 +1 , · · · , xn3 , the known mean curvatures hn1 +1 , · · · , hn2 and the known Gaussian curvatures kn1 +1 , · · · , kn2 to the right-handed side, we can rewrite the system as follows ⎧ ∂Xn0 (t) (p) (m) (g) (0) ⎪ ⎪ ⎨ Mn0 ∂t +Ln0 Xn0 (t) + Ln1 Hn1 (t) + Ln1 Kn1 (t) = B , ⎪ ⎪ ⎩
(1)
(1)
(1)
(2)
Mn1 Hn1 (t)+Ln0 Xn0 (t) = B (1) ,
Mn1 Kn1 (t)+Ln0 Xn0 (t) = B (2) . (1)
(2)
Note that Mn0 and Mn1 are symmetric and positive definite matrices.
(22)
Finite Element Methods for Geometric Modeling
4.2
173
Finite Element Spaces
Now let us introduce two finite element spaces Eh and Fh for representing the surfaces and the curvatures, respectively. Let Eh = span[φ1 , · · · , φn ], Fh = span[ψ1 , · · · , ψn ], where the basis functions φi and ψi are defined for the vertex xi . In our implementation, φi is defined as the limit function of the Loop’s subdivision scheme for the control values one at xi and zero at other vertices. It is known that φi is a quartic box spline. If vertex xi is irregular, local subdivision is needed around xi until the parameter values of interest are interior to a regular patch (see [1] for details). The size of the support of φi is two. Our ψi is defined as a piecewise linear function which has the values one at xi and zero at other vertices. Hence the size of the support of ψi is one. Since ψi is a piecewise linear function, its evaluation is straightforward. For evaluating φi , we employ the fast computation scheme proposed by Stam [23]. 4.3
Temporal Direction Discretization (k)
(k)
Suppose we have approximate solutions Xn0 = Xn0 (tk ), Kn1 = Kn1 (tk ) and (k) Hn1 = Hn1 (tk ) of system (22) at tk . We want to construct approximate solutions (k+1) (k+1) (k+1) Xn0 , Kn1 and Hn1 for the next time step tk+1 = tk + τ (k) using a ∂X − semi-implicit Euler scheme. That is, the derivative ∂tn0 is replaced by [Xnk+1 0 Xnk0 ]/τ (k) , and all the M , L and B matrices in (22) are computed using the surface data at tk . Such a treatment yields a linear system of the equations with (k+1) (k+1) (k+1) Xn0 , Hn1 and Kn1 as unknowns: ⎡ ⎤ ⎡ (k+1)⎤ ⎡ (p) (m) (g) (1) (k)⎤ Mn0 + τ (k) Ln0 τ (k) Ln1 τ (k) Ln1 Xn0 τ (k) B (0)+ Mn0 Xn0 ⎢ ⎥ ⎢ (k+1)⎥ ⎣ (1) (1) ⎦. Ln0 Mn1 0 ⎦ ⎣Hn1 ⎦= B (1) ⎣ (2) (1) (k+1) (2) B 0 Mn1 Kn1 Ln0
(23)
The coefficient matrix of this system is highly sparse, hence an iterative method for solving it is desirable. We use Saad’s iterative method [20], named GMRES, to solve the system. The experiment shows that this iterative method works very well.
5
Examples
In this section, we give several examples to show the strength of the proposed approach. Considering the approach is general, there are infinitely many possibilities for choosing f as well as the geometric models. The examples provided here are just a few of them. In the illustrative figures in this section, τ and T stand for the temporal step-size and the number of iterations used. We present several application examples, including surface blending, hole filling and surface mesh refinement.
174
G. Xu
(a) Input
(b) Initial
(c) f = 1
(d) f = H 2
(e) f = K 2
(f) f = H 2 + K 2
(g) f = H 2 K 2
(h) f = k14 + k24
Fig. 2. Surface Blending: The input is four cylinders to be blended as shown in figure (a). The initial blending surface is shown in figure (b). We evolve the initial surface using the geometric flows with different Lagrange functions f (H, K). Figures (c)-(h) are the evolution results of the corresponding flows, where figure (c) is the result after 135 iterations with temporal step-size 0.01 and figures (d)–(h) show the results after 1000 iterations with temporal step-size 0.00001.
1. Surface Blending Given a collection of surface meshes with boundaries, we construct a fair surface to blend the meshes with the G0 or G1 continuity at the boundaries. The aim is to observe the smoothness at the blending boundaries. Fig. 2 shows some the blending results, where (a) is the input surface model to be blended. (b) shows the initial blending surface construction with noise added. (c)–(h) show the results using different f . It could be observed that all the used fourth order flows yield smooth joining bending surfaces at the boundaries. But the result yielded by f = K 2 is not good at all. The curves on these surfaces are isophotes. The smooth isophotes imply that the surfaces are at least G1 smooth. 2. Hole Filling Given a surface mesh with holes, we construct fair surfaces to fill the holes with the G0 or G1 continuity on the boundaries. Fig 3 shows such an example, where a bunny mesh with a lot of complex shaped holes is given (figure (a)). An initial G0 minimal surface filler construction of the holes is shown in (b), which is generated by the geometric flow with Lagrange function f = 1. The fair blending surface (figure (c)) is generated by taking f = H 2 . The mesh of (d) is the same as (c) but shaded in one color. It is easy to see that the fourth order flow used generate G1 smooth filling surface, while the second order flow yields G0 smooth surface.
Finite Element Methods for Geometric Modeling
(a)
(b)
(c)
175
(d)
Fig. 3. Hole filing: (a) shows a bunny mesh riddled with all complex shaped holes. (b) shows a minimal surface filler construction. (c) (shaded in different colors) and (d) (shaded in one color) are the faired filler surfaces, after 50 iterations, generated by taking the Lagrange function f = H 2 . The time step-length is chosen to be 10−6 .
(a) input
(b) refine once
(e) input
(f) subdivide once
(c) refine twice
(d) refine 3 times
(g) subdivide twice (h) subdivide 3 times
Fig. 4. Mesh refinement: (a) shows an input initial triangular mesh. (b)-(d) are the iterative refinements, where (b), (c) and (d) are generated using τ = 0.005, 0.0001, 0.000001, respectively. The iteration numbers are 15, 50, 50 respectively. (e) shows the same mesh as (a) but with isophotes. (f)–(h) are the repeatedly subdivision results using Butterfly subdivision scheme.
3. Surface Mesh Refinement Mesh refinement is a process of subdividing each triangles into several subtriangles by inserting new vertices at certain places. For simplicity, we subdivide each triangle into 4 subtriangles by inserting a new vertex for each edge. The
176
G. Xu
position of the newly inserted vertices are determined by the geometric PDEs. The original vertices are fixed. The first row of Fig. 4 shows a sequence of meshes, where (a) shows the input triangular surface mesh. Figures (b)-(d) show the iteratively subdivision results using a geometric PDE by taking f = H 2 . For comparing with the stationary interpolatory subdivision scheme, we present in the second row the subdivision results of Butterfly scheme. The isophotes on the surfaces illustrate that the geometric PDE approach yields high quality surfaces.
6
Conclusions
We have developed a mixed finite method for solving a general form fourth order geometric partial differential equation. We applied the proposed method to solve several surface modeling problems, including surface blending, hole filling and mesh refinement with the G1 continuity. The used geometric partial differential equation is universal, containing several well-known geometric partial differential equations as its special cases. The proposed method is general which can be used to construct surfaces for geometric design as well as simulate the behaviors of various geometric PDEs. Experimental results have shown that the proposed approach can be used to handle a large number of geometric PDEs and the numerical algorithm is simple, efficient and gives very desirable results.
References 1. Bajaj, C., Xu, G.: Anisotropic diffusion of surface and functions on surfaces. ACM Transaction on Graphics 22(1), 4–32 (2003) 2. Bajaj, C., Xu, G., Warren, J.: Acoustics Scattering on Arbitrary Manifold Surfaces. In: Proceedings of Geometric Modeling and Processing, Theory and Application, Japan, pp. 73–82 (2002) 3. Bloor, M.I.G., Wilson, M.J.: Generating blend surfaces using partial differential equations. Computer Aided Design 21(3), 165–171 (1989) 4. Bloor, M.I.G., Wilson, M.J.: Using partial differential equations to generate freeform surfaces. Computer Aided Design 22(4), 221–234 (1990) 5. Chavel, I.: Riemannian Geometry – a Modern Introduction. Cambridge University Press, Cambridge (1993) 6. Cirak, F., Ortiz, M.: C 1 -conforming subdivision elements for finite deformation thin-shell analysis. Internat. J. Numer. Methods Engrg. 51(7), 813–833 (2001) 7. Cirak, F., Ortiz, M., Schroder, P.: Subdivision Surfaces: A New Paradigm for ThinShell Finite-Element Analysis. Internat. J. Numer. Methods Engrg. 47, 2039–2072 (2000) 8. Clarenz, U., Diewald, U., Dziuk, G., Rumpf, M., Rusu, R.: A finite element method for surface restoration with boundary conditions. Computer Aided Geometric Design 21(5), 427–445 (2004) 9. Clarenz, U., Diewald, U., Rumpf, M.: Anisotropic geometric diffusion in surface processing. In: Proceedings of Viz2000, IEEE Visualization, Salt Lake City, Utah, pp. 397–405 (2000)
Finite Element Methods for Geometric Modeling
177
10. Deckelnick, K., Dziuk, G.: A fully discrete numerical scheme for weighted mean curvature flow. Numerische Mathematik 91, 423–452 (2002) 11. Desbrun, M., Meyer, M., Schr¨ oder, P., Barr, A.H.: Implicit fairing of irregular meshes using diffusion and curvature flow. In: SIGGRAPH 1999, Los Angeles, USA, pp. 317–324 (1999) 12. do Carmo, M.P.: Differential Geometry of Curves and Surfaces. Prentice-Hall, Englewood Cliffs (1976) 13. do Carmo, M.P.: Riemannian Geometry. Birkh¨ auser, Boston, Basel, Berlin (1992) 14. Du, H., Qin, H.: Direct manipulation and interactive sculpting of PDE surfaces. 19(3), 261–270 (2000) 15. Du, H., Qin, H.: Dynamic PDE-based surface design using geometric and physical constraint. Graphical Models 67(1), 43–71 (2005) 16. Dziuk, G.: An algorithm for evolutionary surfaces. Numerische Mathematik 58, 603–611 (1991) 17. Giaquinta, M., Hildebrandt, S.: Calculus of Variations. A Series of Comprehensive Studies in Mathematics, vol. I(310). Springer, Berlin (1996) 18. Kobbelt, L., Hesse, T., Prautzsch, H., Schweizerhof, K.: Iterative Mesh Generation for FE-computation on Free Form Surfaces. Engng. Comput. 14, 806–820 (1997) 19. Lowe, T., Bloor, M., Wilson, M.: Functionality in blend design. Computer-Aided Design 22(10), 655–665 (1990) 20. Saad, Y.: Iterative Methods for Sparse Linear Systems. Second Edition with corrections (2000) 21. Schneider, R., Kobbelt, L.: Generating Fair Meshes with G1 Boundary conditions. In: Geometric Modeling and Processing, Hong Kong, China, pp. 251–261 (2000) 22. Schneider, R., Kobbelt, L.: Geometric fairing of irregular meshes for free-form surface design. Computer Aided Geometric Design 18(4), 359–379 (2001) 23. Stam, J.: Fast Evaluation of Loop Triangular Subdivision Surfaces at Arbitrary Parameter Values. In: SIGGRAPH 1998 Proceedings (1998), CD-ROM supplement 24. Ugail, H., Bloor, M., Wilson, M.: Techniques for interactive design using the PDE method. ACM Transaction on Graphics 18(2), 195–212 (1999) 25. Xu, G., Pan, Q., Bajaj, C.L.: Discrete surface modelling using partial differential equations. Computer Aided Geometric Design 23(2), 125–145 (2006) 26. Xu, G., Pan, Q.: G1 Surface Modelling Using Fourth Order Geometric Flows. Computer-Aided Design 38(4), 392–403 (2006) 27. Xu, G., Zhang, Q.: Construction of Geometric Partial Differential Equations in Computational Geometry. Mathematica Numerica Sinica 28(4), 337–356 (2006) 28. Xu, G., Zhang, Q.: A General Framework for Surface Modeling Using Geometric Partial Differential Equations. In: Computer Aided Geometric Design (to appear, 2007) 29. Zhang, Q., Xu, G.: Geometric partial differential equations for minimal curvature variation surfaces. In: Research Report No. ICM-06-03. Institute of Computational Mathematics, Chinese Academy of Sciences (2006)
Abstract. A variational formulation of a general form fourth order geometric partial differential equation is derived, and based on which a mixed finite element method is developed. Several surface modeling problems, including surface blending, hole filling and surface mesh refinement with the G1 continuity, are taken into account. The used geometric partial differential equation is universal, containing several well-known geometric partial differential equations as its special cases. The proposed method is general which can be used to construct surfaces for geometric design as well as simulate the behaviors of various geometric PDEs. Experimental results show that it is simple, efficient and gives very desirable results. Keywords: Geometric PDE, Surface blending, hole filling, Surface mesh refinement, Mixed finite element method.
1
Introduction
In this paper, we use a general form geometric partial differential equation (GPDE) to solve several surface modeling problems. The GPDE we use includes many well known equations as its special cases. It is nonlinear and geometrically intrinsic, which is solved using a mixed finite element method. The problems we solve include surface blending, hole filling and mesh refinement with the G1 continuity. For the problems of surface blending and hole filling, we are given triangular surface meshes of the surrounding area. Triangular surface patches need to be constructed to fill the openings enclosed by the surrounding surface mesh and interpolate its boundary with the G1 continuity (see Fig. 2 and 3). For mesh refinement problem, we are given an initial surface mesh. We construct a sequence of surface meshes which interpolates the vertices of the previous meshes with the specified G1 continuity (see Fig. 4). There are basically two approaches for solving a GPDE numerically on surface. One is the generalized finite divided difference method (see [25,28]), the
Project support in part by NSFC grant 60773165 and National Key Basic Research Project of China (2004CB318000).
F. Chen and B. J¨ uttler (Eds.): GMP 2008, LNCS 4975, pp. 164–177, 2008. c Springer-Verlag Berlin Heidelberg 2008
Finite Element Methods for Geometric Modeling
165
other is the finite element method (see [1,10,16,18,26]). The approach we adopt in this paper is a mixed finite element method consisting of using the Loop’s subdivision basis for representing the surfaces and the piecewise linear basis for representing the curvatures. Loop’s basis have been used in finite element analysis and design for thin-shell structures (see [6,7]). The generalized finite divided difference method has been used in [25,28]. It is well known that the finite divided difference method is simpler and easier to implement, but it lacks the convergence analysis. The finite element method is not as simple as the finite divided difference method, but is based on a well developed mathematical foundation. Previous Work. The well-known biharmonic equation and its various variations have been repeatedly used to solve the problems of geometry modeling, such as blending and the hole filling (see [3,4]), interactive surface design (see [15], [24]) and interactive sculpting (see [14]). A PDE method based on the biharmonic equation is proposed [19] with certain functional constraints, such as geometric constraints, physical criteria and engineering requirements, which can be incorporated into the shape design process. The main advantageous of using the biharmonic equation is that the equation is linear, hence it is easy to solve. However, the biharmonic equation is not geometrically intrinsic and its solution depends on the concrete surface parametrization used. Furthermore, the biharmonic equation is defined on a rectangular domain in general, hence it is inappropriate for modeling surfaces with arbitrary shaped boundaries. In contrast, the equations we use in this paper are geometrically intrinsic without the limitation of the shape of the surface boundaries. These equations are called geometric partial differential equations (GPDE). The simplest GPDE may be the mean curvature flow (MCF). MCF and its variations have been intensively used to smooth or fair noisy surfaces (see [1,9,11] for references). These second order equations have been shown to be the most important and efficient flows for fairing or denoising surface meshes. However, for solving the surface modeling and design problems, MCF cannot achieve the G1 continuity at the boundary. To achieve smooth joining of different surface patches in solving the surface blending problem (see [8,21,22,25]), the fourth order geometric flows, such as Willmore flow and surface diffusion flow, have been employed. For instance, the surface diffusion flow is used in [21] to fair surface meshes with the G1 conditions in a special case where the meshes are assumed to have the subdivision connectivity. This equation is also used in [22] for smoothing meshes while satisfying the G1 boundary conditions and the requirement of the outer fairness (the smoothness in the classical sense) and the inner fairness (the regularity of the vertex distribution). Willmore flow has been used for the surface restoration (see [8]) and surface modeling (see [26]) using the finite element method. Three fourth order flows (surface diffusion flow, Willmore flow and quasi surface flow) are treated separately using finite element method in [26]. Main Results. We derive a weak form formulation for a general form fourth order geometric flow. Based on this weak form a mixed form formulation is
166
G. Xu
proposed. The equation is solved by a newly proposed mixed finite element method over surfaces. The proposed approach is simple, efficient and general as well. It handles a class of geometric PDEs in the same way and solves several surface modeling problems in the same manner. The approach is applied to solving each of the surface modeling problems and gives very desirable results for a range of the complicated surface models. The rest of the paper is organized as follows: Section 2 introduces some background material from differential geometry and our earlier publications. Section 3 describes the used nonlinear GPDE and its variational form. In section 4, we give the details of the discretization of the variational form in both the spatial and the temporal directions. Comparative examples to illustrate the different effects achievable from the solutions of these GPDEs are given in section 5.
2
Differential Geometry Preliminaries
In this section, we introduce the used notations, curvatures and several geometric differential operators, including Laplace-Beltrami operator and GiaquintaHildebrandt operator etc. Some results used in this paper are also presented. Let S = {x(u, v) ∈ R3 ; (u, v) ∈ Ω ⊂ R2 } be a regular and sufficiently smooth parametric surface. To simplify the notation we sometimes write w = (u, v) = (u1 , u2 ). Let gαβ = xuα , xuβ and bαβ = n, xuα uβ be the coefficients of the first and the second fundamental forms of S with ∂x ∂2x , xuα uβ = , α, β = 1, 2, n = (xu × xv )/xu × xv , α ∂u ∂uα ∂uβ where ·, · and × denote the usual scalar and cross products of two vectors, respectively, in Euclidean space R3 . Set xuα =
g = det[ gαβ ],
[ g αβ ] = [ gαβ ]−1 , b = det[ bαβ ],
[ bαβ ] = [ bαβ ]−1 .
Curvatures. To introduce the notions of mean curvature and Gaussian curvature, let us first introduce the concept of Weingarten map or shape operator. The shape operator of surface S is a self-adjoint linear map on the tangent space Tx S := span{xu , xv } defined by W : Tx S → Tx S, such that W(xu ) = −nu , W(xv ) = −nv .
(1)
We can easily represent this linear map by a 2 × 2 matrix S = [ bαβ ][ g αβ ]. In particular, [nu , nv ] = −[xu , xv ]S T is valid. The eigenvalues k1 and k2 of S are principal curvatures of S and their arithmetic average and product are the mean curvature and the Gaussian curvature, respectively. That is H=
tr(S) k1 + k2 = , 2 2
K = k1 k2 = det(S).
Tangential gradient operator. Let f be a C 1 smooth function on S. Then the tangential gradient operator ∇ acting on f is given by (see [12], page 102) ∇f = [xu , xv ][ g αβ ][fu , fv ]T ∈ R3 ,
(2)
Finite Element Methods for Geometric Modeling
167
where fu and fv denote the first order partial derivatives of f with respect to arguments u and v. For a vector-valued function f = [f1 , · · · , fk ]T ∈ C 1 (S)k , we define its gradient by ∇f = [∇f1 , · · · , ∇fk ] ∈ R3×k . It is easy to see that ∇x = [xu , xv ][g αβ ][xu , xv ]T , ∇n = −[xu , xv ][g αβ ]S[xu , xv ]T ,
(3)
and both ∇x and ∇n are symmetric 3 × 3 matrices. Second tangent operator. Let f be a C 1 smooth function on S. Then we define the second tangent operator 3 acting on f by (see [27]) 1 b22 −b12 . (4) 3f = [xu , xv ][hαβ ][fu , fv ]T ∈ R3 , with [hαβ ] = g −b12 b11 Divergence operator. Let v be a C 1 smooth vector field on S. Then the divergence of v is defined by ∂ ∂ √ 1 , g [ g αβ ] [xu , xv ]T v . (5) div(v) = √ g ∂u ∂v For a matrix-valued function V = [v1 , · · · , vk ] ∈ C 1 (S)3×k , we define div(V) = [div(v1 ), · · · , div(vk )]T ∈ Rk . Laplace-Beltrami operator. Let f be a C 2 smooth function on S. Then ∇f is a smooth vector field on S. The Laplace-Beltrami operator (LBO) Δ applying to f is defined by (see [13], page 83) Δf = div(∇f ). Giaquinta-Hildebrandt operator. Let f be a C 2 smooth function on S. Then the Giaquinta-Hildebrandt operator acting on f is given by (see [17], page 84) 2f = div(3f ). For a vector-valued function f = [f1 , · · · , fk ]T ∈ C 1 (S)k , we define Δf = div(∇f ),
2f = div(3f ).
(6)
For the operators introduced above, we can prove the following theorems that are used in the sequel. Detailed proofs can be found in our technical report [29]. Theorem 1. Let x ∈ S. Then Δx = 2Hn, Δn = −2∇H − 2HΔx + 2x.
(7)
Theorem 2. Let x ∈ S. Then 2x = 2Kn, 2n = −∇K − KΔx = −∇K − H2x.
(8)
To derive the weak form of the considered geometric flows, we introduce two Green’s formulas.
168
G. Xu
Theorem 3 (Green’s formula I for LB operator). ([5], page 142) Let f ∈ C 2 (S), h ∈ C 1 (S), with at least one of them compactly supported. Then hΔf + ∇f, ∇h dA = 0. M
Theorem 4 (Green’s formula I for GH operator). Let f ∈ C 2 (S), h ∈ C 1 (S), with at least one of them compactly supported. Then hf + ∇f, 3h dA = 0. M
3
Weak and Mixed Forms of a General Fourth Order GPDE
Euler-Lagrange equation. Let f (H, K) ∈ C 1 (R × R) be a given Lagrange function. Define energy functional f (K, H)dA. (9) E (S) = S
In this paper we use a vector-valued Euler-Lagrange equation derived from the complete-variation x(w, ε) = x(w) + εΦ(w), w ∈ Ω, |ε| 1, Φ ∈ Cc∞ (Ω)3 ,
(10)
for the functional E (S). We refer a result from [28] as the following theorem: Theorem 5. Let f (H, K) ∈ C 1 (R × R). Then the Euler-Lagrange equation of the integral E (S) from the complete-variation (10) is 1 2(fK n) + Δ(fH n) − div(fH ∇n) − div[(f − 2KfK )∇x]) = 0 ∈ R3 . (11) 2 Theorem 6. The weak form of (11) can be written as
1 fH ∇x3φ + (f − 2HfH − 2KfK )∇x∇φ + fH nΔφ + fK n2φ dA = 0, 2 S (12) ∀φ ∈ Cc∞ (Ω). Proof. By direct calculation, it is not difficult to derive that [g αβ ]S = 2H[g αβ ] − K[bαβ ]. Using this equality and (3), we have ∇n = 3x − 2H∇x.
(13)
Multiplying both sides of (11) with φ ∈ Cc∞ (Ω), then taking integral over S, using Green formulas, and finally using (13), we have
1 2(fK n) + Δ(fH n) − div(fH ∇n) − div[(f − 2KfK )∇x]) φ dA 0= 2 S
1 fK n2φ + fH nΔφ + [fH (3x−2H∇x)]∇φ+[(f −2KfK )∇x]∇φ) dA. = 2 S Noticing that 3x∇φ = ∇x3φ, we obtain (12) in the end.
Finite Element Methods for Geometric Modeling
169
A General form Geometric Flow. For a given function f (H, K), let S0 be a given initial surface with boundary Γ . Then the geometric flow consists of finding a family {S(t); t ≥ 0} of smooth orientable surfaces in R3 which evolve according to the following equation ∂x 1 ∂t + 2(fK n) + 2 Δ(fH n) − div(fH ∇n) − div (f − 2KfK )∇x = 0, (14) S(0) = S0 , ∂S(t) = Γ. Weak form formulation. Using (12), we can obtain the weak form of equation (14) as follows: ⎧ ∂x ⎪ ⎪ fH ∇x3φ + (f − 2HfH − 2KfK )∇x∇φ φ dA + ⎪ ⎨ S ∂t S
+ 12 fH nΔφ + fK n2φ dA = 0, ∀φ ∈ Cc∞ (Ω), ⎪ ⎪ ⎪ ⎩ S(0) = S0 , ∂S(t) = Γ.
(15)
Assumptions. In order to make the generated flow meaningful and well-defined. We first require that f is a smooth function about its arguments. We also assume that f is an algebraic function, meaning it does not involve differential and integral operations about its variables. Suppose fH and fK could be represented as fH = 2αH + 2βK + μ, fK = 2γH + 2δK + ν. (16) Let h = 2Hn, k = 2Kn. Then using (15) and hψ dA = Δxψ dA = − ∇x∇ψ dA, ∀ψ ∈ Cc∞ (Ω), S S S kϕ dA = 2xϕ dA = − ∇x3ϕ dA, ∀ϕ ∈ Cc∞ (Ω), S
S
S
we obtain a mixed form of equation (14) as follows: Find [x, h, k] ∈ R3×3 , such that ⎧ ∂x ⎪ ⎪ φ dA + fH ∇x3φ + (f − 2HfH − 2KfK )∇x∇φ ⎪ ⎪ ∂t ⎪ S S ⎪ ⎪ 1 ⎪ ⎪ + (αh + βk + νn)Δφ+ (γh + δk + μn)2φ dA = 0, ∀φ ∈ Cc∞ (Ω), ⎨ 2 ⎪ hψ dA + ∇x∇ψ dA = 0, ∀ψ ∈ Cc∞ (Ω), ⎪ ⎪ ⎪ S S ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ kϕ dA + ∇x3ϕ dA = 0, ∀ϕ ∈ Cc∞ (Ω). S
S
(17) From the definitions of tangential operators, it is easy to derive that ∇x3φ = 3φ, ∇x∇φ = ∇φ. Further notice that 2H = nT h, 2K = nT k, the weak form (17) can be rewritten as: Find [x, h, k] ∈ R3×3 , such that
170
G. Xu
⎧ ∂x ⎪ ⎪ μ∇x3φ + η∇x∇φ + 3φnT (αh + βk) φ dA + ⎪ ⎪ ⎪ S ∂t S ⎪ ⎪ ⎪ +∇φnT [(ω − fH )h + ( − fK )k] ⎪ ⎪
⎪ 1 ⎨ + (αh + βk + νn)Δφ+ (γh + δk + μn)2φ dA = 0, 2 ⎪ ⎪ ⎪ hψ dA + ∇x∇ψ dA = 0, ∀ψ ∈ Cc∞ (Ω), ⎪ ⎪ ⎪ ⎪ S S ⎪ ⎪ ⎪ ⎪ ⎩ kϕ dA + ∇x3ϕ dA = 0, ∀ϕ ∈ Cc∞ (Ω), S
∀φ ∈ Cc∞ (Ω),
S
(18) where ω, and η are defined by f = 2ωH + 2K + η. System (18) is the starting point of the finite element discretization.
4
Solving the GPDE by Mixed Finite Element Methods
Consider a triangular surface mesh Ω with vertices x1 , x2 , · · · , xn . For each vertex xi , we associate it with a basis function φi . Then the surface S is approximately represented as x=
n
φj (q)xj ∈ R3 , q ∈ Ω,
(19)
j=1
where x1 , x2 , · · · , xn are regarded as the control vertices of S. Now we classify these control vertices into several categories. The first category consists of interior vertices, denoted as x1 , x2 , · · · , xn0 . The positions of these vertices are subject to change (unknown). For the problems of hole filling and blending, the interior vertices are those in the openings (see Fig. 1(a), the solid dots). For the problem of mesh refinement, all the non-interpolating vertices are interior (see Fig. 1(b), the solid dots). Apart from the interior vertices, the remaining vertices are classified as 1-ring neighbors (see Fig. 1(a), the empty
Ω0
Ω1
1
1
(a)
(b)
Fig. 1. Classification of vertices. (a) The solid dots are interior control vertices. The empty dots are 1-ring neighbor (boundary) vertices. The solid squares are 2-ring neighbor vertices. The shaded region is Ω0 . The region enclosed by the 2-ring neighbor vertices is Ω1 . (b) The solid dots are the interior control vertices. The empty ones are the interpolated vertices.
Finite Element Methods for Geometric Modeling
171
dots), 2-ring neighbors (see Fig. 1(a), the solid squares), · · · , according to their distance to the interior vertices. For instance, the 1-ring neighbors are all the vertices adjacent to the interior ones, the 2-ring neighbors are all the remaining vertices adjacent to the 1-ring neighbors, · · · . The k-ring neighbor control vertices are denoted as xnk−1 +1 , · · · , xnk . All the 1-ring, 2-ring, · · · , neighbor vertices are fixed. In our mixed form formulations of GPDE, the mean and Gaussian curvature normals are also treated as unknown vector-valued functions. Hence, for each control vertex xi , we also associate it with a continuous basis function ψi for representing the curvatures. For instance, mean and Gaussian curvature normals are represented as h(x) =
n
ψj (q)hj ∈ R3 ,
k(x) =
j=1
n
ψj (q)kj ∈ R3 , q ∈ Ω,
(20)
j=1
respectively. Here hj and kj denote the mean and Gaussian curvature normals at xj , respectively. Suppose that φi (q) and ψi (q) have compact supports and the sizes of their supports are 2 and 1, respectively. At a k-ring neighbor vertex, the surface point does not depend upon the interior vertices if k ≥ 2, so the mean and Gaussian curvature normals there are well defined from the surface. However, at the interior and the 1-ring neighbor vertices, the corresponding surface points relate to unknown vertices. Then the mean and Gaussian curvatures there cannot be determined from the surface. Hence we treat these curvatures as unknowns. In a word, the unknown vertices are x1 , x2 , · · · , xn0 , the unknown mean curvature normals are h1 , h2 , · · · , hn0 , hn0 +1 , · · · , hn1 , and the unknown Gaussian curvature normals are k1 , k2 , · · · , kn0 , kn0 +1 , · · · , kn1 . In the next two sub-sections, we will obtain matrix forms of the used GPDE by substituting (19)–(20) into (18). Note that the proposed scheme needs the values of the mean and Gaussian curvature normals at each vertex. For an interior or 1-ring neighbor vertex, the mean and Gaussian curvature normals are unknown in the equation to be solved, but initial values at t = 0 is required. For a k-ring neighbor vertex (k ≥ 2), these vectors need to be computed previously and keep fixed afterwards (see Remark 4.3 for how these vectors are computed). A default choice of these vectors is to set them as zeros. 4.1
Spatial Direction Discretizations
We construct one equation for each unknown vertex xj (j = 1, · · · , n0 ) and other two equations for the unknown vector-valued curvatures hj and kj (j = 1, · · · , n1 ). This is achieved by the following operations: (i) Substituting (19)–(20) into (18). (ii) Taking the test functions φ = φi (i = 1, · · · , n0 ) and
172
G. Xu ∂x (t)
j ψ = ϕ = ψi (i = 1, · · · , n1 ). (iii) Using the fact that ∂t = 0 for j > n0 (since xj is fixed). We then obtain the following matrix form of (18): ⎧ ∂Xn0 (t) (p) (m) (g) ⎪ + Ln3 Xn3 (t) + Ln2 Hn2 (t) + Ln2 Kn2 (t) + Bn0 = 0, ⎪ ⎨ Mn0 ∂t
⎪ ⎪ ⎩
(1)
(1)
Mn2 Hn2 (t) + Ln3 Xn3 (t) = 0,
(1) Mn2 Kn2 (t)
where
+
(2) Ln3 Xn3 (t)
(21)
= 0,
Xj (t) = [x1 (t), · · · , xj (t)]T ∈ Rj×3 , Hj (t) = [h1 (t), · · · , hj (t)]T ∈ Rj×3 , Kj (t) = [k1 (t), · · · , kj (t)]T ∈ Rj×3 ,
are unknowns and n ,K
0 MK = (mij )ij=1 , n0 ,K (p) (p) , LK = lij nij=1 1 ,K (1) (1) , LK = lij
ij=1
n1 ,K (1) (1) MK = mij , nij=1 n0 ,K 0 ,K (m) (m) (g) (g) LK = lij , LK = lij , ij=1 ij=1 n1 ,K n0 (2) (2) LK = lij , Bn0 = bTi i=1 ij=1
are the coefficient matrices. The elements of these matrices are defined as follows: (1) mij = φi φj dA, mij = ψi ψj dA, S S (p) lij = μ(∇ψj )T 3φi + η(∇ψj )T ∇φi dA ∈ R, S
1 (m) ψj [α3φi + (ω − fH )∇φi ]nT + ( αψj Δφi + γψj 2φi )I3 dA ∈ R3×3 , lij = 2 S
1 (g) ψj [β3φi + ( − fK )∇φi ]nT + ( βψj Δφi + δψj 2φi )I3 dA ∈ R3×3 . lij = 2 S (1) (2) T lij = (∇φj ) ∇ψi dA , lij = (∇φj )T 3ψi dA, S S 1 bi = [n( νΔφi + μ2φi )]dA. 2 S The integrations in these matrix elements are computed by a 12-point Gaussian quadrature rule. That is, each domain triangle is subdivided into four subtriangles and a three-point Gaussian quadrature rule is employed on each of the sub-triangles (see [2] for a set of Gaussian quadrature rules). In system (21), moving the terms relating to the known vertices xn0 +1 , · · · , xn3 , the known mean curvatures hn1 +1 , · · · , hn2 and the known Gaussian curvatures kn1 +1 , · · · , kn2 to the right-handed side, we can rewrite the system as follows ⎧ ∂Xn0 (t) (p) (m) (g) (0) ⎪ ⎪ ⎨ Mn0 ∂t +Ln0 Xn0 (t) + Ln1 Hn1 (t) + Ln1 Kn1 (t) = B , ⎪ ⎪ ⎩
(1)
(1)
(1)
(2)
Mn1 Hn1 (t)+Ln0 Xn0 (t) = B (1) ,
Mn1 Kn1 (t)+Ln0 Xn0 (t) = B (2) . (1)
(2)
Note that Mn0 and Mn1 are symmetric and positive definite matrices.
(22)
Finite Element Methods for Geometric Modeling
4.2
173
Finite Element Spaces
Now let us introduce two finite element spaces Eh and Fh for representing the surfaces and the curvatures, respectively. Let Eh = span[φ1 , · · · , φn ], Fh = span[ψ1 , · · · , ψn ], where the basis functions φi and ψi are defined for the vertex xi . In our implementation, φi is defined as the limit function of the Loop’s subdivision scheme for the control values one at xi and zero at other vertices. It is known that φi is a quartic box spline. If vertex xi is irregular, local subdivision is needed around xi until the parameter values of interest are interior to a regular patch (see [1] for details). The size of the support of φi is two. Our ψi is defined as a piecewise linear function which has the values one at xi and zero at other vertices. Hence the size of the support of ψi is one. Since ψi is a piecewise linear function, its evaluation is straightforward. For evaluating φi , we employ the fast computation scheme proposed by Stam [23]. 4.3
Temporal Direction Discretization (k)
(k)
Suppose we have approximate solutions Xn0 = Xn0 (tk ), Kn1 = Kn1 (tk ) and (k) Hn1 = Hn1 (tk ) of system (22) at tk . We want to construct approximate solutions (k+1) (k+1) (k+1) Xn0 , Kn1 and Hn1 for the next time step tk+1 = tk + τ (k) using a ∂X − semi-implicit Euler scheme. That is, the derivative ∂tn0 is replaced by [Xnk+1 0 Xnk0 ]/τ (k) , and all the M , L and B matrices in (22) are computed using the surface data at tk . Such a treatment yields a linear system of the equations with (k+1) (k+1) (k+1) Xn0 , Hn1 and Kn1 as unknowns: ⎡ ⎤ ⎡ (k+1)⎤ ⎡ (p) (m) (g) (1) (k)⎤ Mn0 + τ (k) Ln0 τ (k) Ln1 τ (k) Ln1 Xn0 τ (k) B (0)+ Mn0 Xn0 ⎢ ⎥ ⎢ (k+1)⎥ ⎣ (1) (1) ⎦. Ln0 Mn1 0 ⎦ ⎣Hn1 ⎦= B (1) ⎣ (2) (1) (k+1) (2) B 0 Mn1 Kn1 Ln0
(23)
The coefficient matrix of this system is highly sparse, hence an iterative method for solving it is desirable. We use Saad’s iterative method [20], named GMRES, to solve the system. The experiment shows that this iterative method works very well.
5
Examples
In this section, we give several examples to show the strength of the proposed approach. Considering the approach is general, there are infinitely many possibilities for choosing f as well as the geometric models. The examples provided here are just a few of them. In the illustrative figures in this section, τ and T stand for the temporal step-size and the number of iterations used. We present several application examples, including surface blending, hole filling and surface mesh refinement.
174
G. Xu
(a) Input
(b) Initial
(c) f = 1
(d) f = H 2
(e) f = K 2
(f) f = H 2 + K 2
(g) f = H 2 K 2
(h) f = k14 + k24
Fig. 2. Surface Blending: The input is four cylinders to be blended as shown in figure (a). The initial blending surface is shown in figure (b). We evolve the initial surface using the geometric flows with different Lagrange functions f (H, K). Figures (c)-(h) are the evolution results of the corresponding flows, where figure (c) is the result after 135 iterations with temporal step-size 0.01 and figures (d)–(h) show the results after 1000 iterations with temporal step-size 0.00001.
1. Surface Blending Given a collection of surface meshes with boundaries, we construct a fair surface to blend the meshes with the G0 or G1 continuity at the boundaries. The aim is to observe the smoothness at the blending boundaries. Fig. 2 shows some the blending results, where (a) is the input surface model to be blended. (b) shows the initial blending surface construction with noise added. (c)–(h) show the results using different f . It could be observed that all the used fourth order flows yield smooth joining bending surfaces at the boundaries. But the result yielded by f = K 2 is not good at all. The curves on these surfaces are isophotes. The smooth isophotes imply that the surfaces are at least G1 smooth. 2. Hole Filling Given a surface mesh with holes, we construct fair surfaces to fill the holes with the G0 or G1 continuity on the boundaries. Fig 3 shows such an example, where a bunny mesh with a lot of complex shaped holes is given (figure (a)). An initial G0 minimal surface filler construction of the holes is shown in (b), which is generated by the geometric flow with Lagrange function f = 1. The fair blending surface (figure (c)) is generated by taking f = H 2 . The mesh of (d) is the same as (c) but shaded in one color. It is easy to see that the fourth order flow used generate G1 smooth filling surface, while the second order flow yields G0 smooth surface.
Finite Element Methods for Geometric Modeling
(a)
(b)
(c)
175
(d)
Fig. 3. Hole filing: (a) shows a bunny mesh riddled with all complex shaped holes. (b) shows a minimal surface filler construction. (c) (shaded in different colors) and (d) (shaded in one color) are the faired filler surfaces, after 50 iterations, generated by taking the Lagrange function f = H 2 . The time step-length is chosen to be 10−6 .
(a) input
(b) refine once
(e) input
(f) subdivide once
(c) refine twice
(d) refine 3 times
(g) subdivide twice (h) subdivide 3 times
Fig. 4. Mesh refinement: (a) shows an input initial triangular mesh. (b)-(d) are the iterative refinements, where (b), (c) and (d) are generated using τ = 0.005, 0.0001, 0.000001, respectively. The iteration numbers are 15, 50, 50 respectively. (e) shows the same mesh as (a) but with isophotes. (f)–(h) are the repeatedly subdivision results using Butterfly subdivision scheme.
3. Surface Mesh Refinement Mesh refinement is a process of subdividing each triangles into several subtriangles by inserting new vertices at certain places. For simplicity, we subdivide each triangle into 4 subtriangles by inserting a new vertex for each edge. The
176
G. Xu
position of the newly inserted vertices are determined by the geometric PDEs. The original vertices are fixed. The first row of Fig. 4 shows a sequence of meshes, where (a) shows the input triangular surface mesh. Figures (b)-(d) show the iteratively subdivision results using a geometric PDE by taking f = H 2 . For comparing with the stationary interpolatory subdivision scheme, we present in the second row the subdivision results of Butterfly scheme. The isophotes on the surfaces illustrate that the geometric PDE approach yields high quality surfaces.
6
Conclusions
We have developed a mixed finite method for solving a general form fourth order geometric partial differential equation. We applied the proposed method to solve several surface modeling problems, including surface blending, hole filling and mesh refinement with the G1 continuity. The used geometric partial differential equation is universal, containing several well-known geometric partial differential equations as its special cases. The proposed method is general which can be used to construct surfaces for geometric design as well as simulate the behaviors of various geometric PDEs. Experimental results have shown that the proposed approach can be used to handle a large number of geometric PDEs and the numerical algorithm is simple, efficient and gives very desirable results.
References 1. Bajaj, C., Xu, G.: Anisotropic diffusion of surface and functions on surfaces. ACM Transaction on Graphics 22(1), 4–32 (2003) 2. Bajaj, C., Xu, G., Warren, J.: Acoustics Scattering on Arbitrary Manifold Surfaces. In: Proceedings of Geometric Modeling and Processing, Theory and Application, Japan, pp. 73–82 (2002) 3. Bloor, M.I.G., Wilson, M.J.: Generating blend surfaces using partial differential equations. Computer Aided Design 21(3), 165–171 (1989) 4. Bloor, M.I.G., Wilson, M.J.: Using partial differential equations to generate freeform surfaces. Computer Aided Design 22(4), 221–234 (1990) 5. Chavel, I.: Riemannian Geometry – a Modern Introduction. Cambridge University Press, Cambridge (1993) 6. Cirak, F., Ortiz, M.: C 1 -conforming subdivision elements for finite deformation thin-shell analysis. Internat. J. Numer. Methods Engrg. 51(7), 813–833 (2001) 7. Cirak, F., Ortiz, M., Schroder, P.: Subdivision Surfaces: A New Paradigm for ThinShell Finite-Element Analysis. Internat. J. Numer. Methods Engrg. 47, 2039–2072 (2000) 8. Clarenz, U., Diewald, U., Dziuk, G., Rumpf, M., Rusu, R.: A finite element method for surface restoration with boundary conditions. Computer Aided Geometric Design 21(5), 427–445 (2004) 9. Clarenz, U., Diewald, U., Rumpf, M.: Anisotropic geometric diffusion in surface processing. In: Proceedings of Viz2000, IEEE Visualization, Salt Lake City, Utah, pp. 397–405 (2000)
Finite Element Methods for Geometric Modeling
177
10. Deckelnick, K., Dziuk, G.: A fully discrete numerical scheme for weighted mean curvature flow. Numerische Mathematik 91, 423–452 (2002) 11. Desbrun, M., Meyer, M., Schr¨ oder, P., Barr, A.H.: Implicit fairing of irregular meshes using diffusion and curvature flow. In: SIGGRAPH 1999, Los Angeles, USA, pp. 317–324 (1999) 12. do Carmo, M.P.: Differential Geometry of Curves and Surfaces. Prentice-Hall, Englewood Cliffs (1976) 13. do Carmo, M.P.: Riemannian Geometry. Birkh¨ auser, Boston, Basel, Berlin (1992) 14. Du, H., Qin, H.: Direct manipulation and interactive sculpting of PDE surfaces. 19(3), 261–270 (2000) 15. Du, H., Qin, H.: Dynamic PDE-based surface design using geometric and physical constraint. Graphical Models 67(1), 43–71 (2005) 16. Dziuk, G.: An algorithm for evolutionary surfaces. Numerische Mathematik 58, 603–611 (1991) 17. Giaquinta, M., Hildebrandt, S.: Calculus of Variations. A Series of Comprehensive Studies in Mathematics, vol. I(310). Springer, Berlin (1996) 18. Kobbelt, L., Hesse, T., Prautzsch, H., Schweizerhof, K.: Iterative Mesh Generation for FE-computation on Free Form Surfaces. Engng. Comput. 14, 806–820 (1997) 19. Lowe, T., Bloor, M., Wilson, M.: Functionality in blend design. Computer-Aided Design 22(10), 655–665 (1990) 20. Saad, Y.: Iterative Methods for Sparse Linear Systems. Second Edition with corrections (2000) 21. Schneider, R., Kobbelt, L.: Generating Fair Meshes with G1 Boundary conditions. In: Geometric Modeling and Processing, Hong Kong, China, pp. 251–261 (2000) 22. Schneider, R., Kobbelt, L.: Geometric fairing of irregular meshes for free-form surface design. Computer Aided Geometric Design 18(4), 359–379 (2001) 23. Stam, J.: Fast Evaluation of Loop Triangular Subdivision Surfaces at Arbitrary Parameter Values. In: SIGGRAPH 1998 Proceedings (1998), CD-ROM supplement 24. Ugail, H., Bloor, M., Wilson, M.: Techniques for interactive design using the PDE method. ACM Transaction on Graphics 18(2), 195–212 (1999) 25. Xu, G., Pan, Q., Bajaj, C.L.: Discrete surface modelling using partial differential equations. Computer Aided Geometric Design 23(2), 125–145 (2006) 26. Xu, G., Pan, Q.: G1 Surface Modelling Using Fourth Order Geometric Flows. Computer-Aided Design 38(4), 392–403 (2006) 27. Xu, G., Zhang, Q.: Construction of Geometric Partial Differential Equations in Computational Geometry. Mathematica Numerica Sinica 28(4), 337–356 (2006) 28. Xu, G., Zhang, Q.: A General Framework for Surface Modeling Using Geometric Partial Differential Equations. In: Computer Aided Geometric Design (to appear, 2007) 29. Zhang, Q., Xu, G.: Geometric partial differential equations for minimal curvature variation surfaces. In: Research Report No. ICM-06-03. Institute of Computational Mathematics, Chinese Academy of Sciences (2006)
