Construction of minimal subdivision surface with a given boundary

Computer-Aided Design 43 (2011) 374–380

Contents lists available at ScienceDirect

Computer-Aided Design journal homepage: www.elsevier.com/locate/cad

Construction of minimal subdivision surface with a given boundary Qing Pan a , Guoliang Xu b,∗ a

College of Mathematics and Computer Science, Hunan Normal University, Changsha, 410081, China

b

LSEC, Institute of Computational Mathematics, Academy of Mathematics and System Sciences, Chinese Academy of Sciences, Beijing 100190, China

article

info

Article history: Received 9 June 2009 Accepted 27 December 2010 Keywords: Minimal subdivision surface Loop subdivision Mean curvature flow

abstract The fascinating characters of minimal surface make it to be widely used in shape design. While the flexibility and high quality of subdivision surface make it a powerful mathematical tool for shape representation. In this paper, we construct minimal subdivision surfaces with given boundaries using the mean curvature flow, a second order geometric partial differential equation. This equation is solved by a finite element method where the finite element space is spanned by the limit functions of an extended Loop’s subdivision scheme proposed by Biermann et al. Using this extended Loop’s subdivision scheme we can treat a surface with boundary, thereby construct the perfect minimal subdivision surfaces with any topology of the control mesh and any shaped boundaries. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction Minimal surfaces have several desirable properties because of which they are often used as models in architecture. First of all, minimal surfaces have the least surface area, which make them to be widely used in large scale and light roof constructions. Secondly, minimal surfaces have separable property, meaning: any sub-patch, no matter how small, cut from a minimal surface still has the least area of all surface sub-patches with the same boundary. Thirdly, minimal surfaces have balanced surface tension, which stabilizes the whole construction since the tension is in equilibrium at each point on a roof, as on a soap film. Finally, there are no umbilicus points on a minimal surface, hence no water could stay on a minimal surface roof. Architecture inspired from minimal surfaces embodies the unite of economy and beauty. The most spectacular buildings of that architectural style are the roofs of the Munich Olympic stadium, the former Kongreßhalle in Berlin, and many of the smaller tent roofs on pavilions in parks and other public places. In art world we can see plenty of ingenious sculpture works playing the ultimate of minimal surfaces. Scientists and engineers have anticipated the nanotechnology applications of minimal surfaces in areas of molecular engineering and materials science. In the area of computer aided design (CAD), minimal surfaces are also studied. To construct a continuous surface with a given boundary, Monterde (see [1]) used a four-sided Bézier surface to approximate a minimal surface. The problem which is called Plateau–Bézier problem was solved by replacing the area functional with the Dirichlet functional. Cosin and Monterde [2] studied

∗

Corresponding author. Tel.: +86 10 62624352; fax: +86 10 62542285. E-mail address: [email protected] (G. Xu).

0010-4485/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.cad.2010.12.013

the properties that the control vertices must satisfy and showed that in the bicubical case all minimal surfaces are, up to an affine transformation, pieces of the Enneper surfaces. Triangular Bézier surface was constructed by Arnal et al. [3] using a variational approach. Excellent work has been done on the use of minimal surfaces in geometric modeling and shape design (see [4–6]). Discrete minimal surfaces were studied by Polthier in [7]. Minimal surfaces were also produced as the steady solution of the mean curvature flow (see [8]) for both continuous and discrete cases. For the continuous case, the constructed surfaces are usually Bézier surfaces, or B-spline surfaces. Obviously, Bézier surfaces, B-spline or NURB surfaces have to be three or four sided. This is a serious limitation for designing minimal surfaces with any shaped boundaries. In this paper, our intention is to construct minimal subdivision surfaces with piecewise B-spline curve boundaries. There is no limitation on the number of spline pieces. B-spline has been widely accepted as a representation tool for curves and surfaces in industrial design. Using B-spline to represent a surface boundary is preferable and acceptable. To represent a surface patch with any topology, perhaps subdivision surfaces are the best candidates, since there is no limitation on the topology of the control mesh. However, subdivision surfaces, such as Loop’s subdivision surface and the Catmull–Clark subdivision surface are traditionally closed, which cannot be used directly for serving our purpose. For many surface modeling problems, such as the construction of bodies of cars, aircrafts, machine parts and roofs, surfaces are usually constructed in a piecewise manner with fixed boundaries. In such a case Loop’s subdivision scheme could not be applied near the boundaries of the control mesh. Therefore, extension of Loop’s subdivision scheme for the control mesh with boundaries is definitely required. On this aspect, an excellent work

Q. Pan, G. Xu / Computer-Aided Design 43 (2011) 374–380

a

b

375

c

Fig. 1.1. (a) is the boundary curve. (b) is the initial construction. (c) is its minimal subdivision surface.

has been done by Biermann et al. [9] and that is just sufficient for our use. Loop’s subdivision has been shown to be a very powerful tool in several areas. For instance, it is used successfully in smooth surface reconstruction from scattered data (see [10,11]) and the thinshell finite element analysis (see [12]) for describing the geometry and the displacement fields. Subdivision techniques can handle arbitrary topology surfaces. The variational subdivision method combines the two approaches for effective designing of freeform surfaces. Kobbelt et al. in [13,14] proposed a mesh refinement method based on variational methods that try to overcome certain drawbacks of uniform subdivision. A new approach of subdivision based on the evolution of surfaces under curvature motion is presented in [15]. We intend to construct the minimal subdivision surfaces (see Fig. 1.1) using the mean curvature flow, a second order geometric partial differential equation. The equation is solved by a finite element method. The contribution of this paper includes the proposal of a method for constructing minimal subdivision surfaces with given boundaries. Several schemes, such as the extended Loop’s subdivision scheme, the fast evaluation of the basis functions and the finite element method for the initialboundary problem of the mean curvature flow with the Dirichlet boundary condition, are combined together to form an efficient and mathematically sound approach. The rest of this paper is organized as follows: In Section 2 we describe some related results of Loop’s subdivision scheme including interior surface patches and boundary surface patches. Section 3 shows the mean curvature flow used in this paper, the details of its discretization and the numerical computation for the solution of this flow. In Section 4 we give several examples of minimal subdivision surface construction and some comparison results for approximation error to illustrate the effects of our method. Section 5 is the conclusion. 2. Loop’s subdivision scheme and its extension We shall discretize the proposed diffusion problem in a function space which is defined by the limit of the extended Loop’s subdivision. In this section we summarize Loop’s subdivision scheme and its extension by Biermann et al. [9] for completeness. 2.1. Loop’s subdivision scheme Beginning with an initial triangular control mesh with arbitrary topology structure, we can achieve a sequence of increasingly refined meshes which consist of only triangles. The vertex with valence 6 is named ordinary point, and the vertex whose valence is not 6 is named extraordinary point. In the refinement step of Loop’s subdivision scheme, each triangle is subdivided into

four sub-triangles. The control vertices of the refined mesh are generated from the control vertices of the previous step by a portfolio of weight coefficients. There are two sets of refined rules: vertex recomputation rule and edge refinement rule which can be referred to [16] for details. It is obvious that the valence of all newly generated vertices is 6, and the vertices inherited from the initial mesh keep their original valences. Finally the obtained sequence of meshes converges to a limit surface composed of unlimited number of surface patches. The limit surface of Loop’s subdivision is C 2 everywhere except at the extraordinary points where it is C 1 . 2.2. Extension of Loop’s subdivision scheme Loop’s subdivision scheme is usable only for closed (without boundaries) triangular control meshes and the subdivision scheme is applied uniformly to all vertices and edges. However, for a lot of surface modeling problems, surfaces are usually constructed in a piecewise manner with fixed boundaries. In such a case, Loop’s subdivision scheme cannot be applied near the boundary of the control mesh. Therefore, an extension of Loop’s subdivision scheme for control mesh with boundaries is required. The first of this type of extension perhaps is the work of Hoppe et al. [10]. Their subdivision scheme results in C 1 continuous surfaces. However, this work could not be used to solve our modeling problem because the limit boundary curves depend on the interior vertices which may lead to gaps between two surface patches. Another extension proposed by Biermann et al. [9] is perfectly suitable for our purpose. The extended subdivision scheme proposed by Biermann et al. is applied to a tagged triangular mesh that could have boundaries. Vertices are classified as corner vertices, boundary vertices and interior vertices whose recomputation rules are reset with new weight coefficients. Edges are tagged as boundary edges, subboundary edges and interior edges. Boundary edges are at the boundary of the control mesh in general, but they could also lie in the interior of the mesh and are called creases (see [9]) (in this paper, we name creases as boundary edges). The sub-boundary edges are the ones that are not boundary edges but incident to boundary vertices. The remainders are called interior ones. This extension also assigns the refinement rules for the boundary edges and sub-boundary edges. In describing the evaluation process for the subdivision surfaces, we classify triangles as boundary triangles, sub-boundary triangles and interior triangles. The triangles containing boundary vertices are named as boundary triangles. The ones adjacent to the boundary triangles are called sub-boundary triangles. All the others are called interior ones. 2.3. Evaluation of Loop’s subdivision surfaces The limit surface of a subdivision scheme is defined by an infinite iteration procedure. There is no closed form for the limit

376

Q. Pan, G. Xu / Computer-Aided Design 43 (2011) 374–380

surface in general. This makes the exact evaluation of the surface at any point a non-trivial task. It is well known that Stam [17] proposed a fast method for evaluating the limit surfaces of Loop’s subdivision scheme at any parameter value. This advanced technique was successively implemented and applied in several other work (see [18,19]). 2.3.1. Evaluating patches for interior triangles Each triangle of the control mesh corresponds to one triangular patch of the limit surface. The triangle is regarded as the parameter domain of the surface patch. If all vertices of the triangle have valence 6 and none of its 2-ring neighbor vertices is a boundary vertex, the resulting limit surface patch is exactly described as a single quartic box spline, for which an explicit closed form exists [17]. If a triangle is irregular, i.e., at least one of its vertices has a valence other than 6 or there is at least one non-interior neighbor vertex, the resulting surface patch is not a quartic box spline. In order to evaluate the surface at any parametric value, the mesh needs to be subdivided repeatedly until the parameter values of interest are interior to a regular triangle. For the evaluation of irregular patches, we use the scheme proposed by Stam [17]. 2.3.2. Evaluating patches for boundary and sub-boundary triangles Evaluating sub-boundary triangular patches. Subdividing a subboundary triangle once using the extended Loop’s subdivision scheme results in four interior triangles. Each of the corresponding patches can be evaluated using the technique described in the previous sub-section. Evaluating boundary triangular patches. Subdividing a boundary triangle once using the extended Loop’s subdivision scheme results in three (respectively one) boundary triangles and one (respectively three) sub-boundary triangle. The patches for subboundary triangles are evaluated using the above mentioned method. The boundary triangles may need to be further subdivided if the parameter values, where the surface patch needs to be evaluated, are in this domain. This process is carried out repeatedly till the parameter values to be evaluated are within a sub-boundary triangle. k subdivisions are needed to compute the kth power of the subdivision matrix. This can be computed using the method presented in [20]. 3. Minimal surface construction using finite element method We use the mean curvature flow to construct the minimal surface. In this section we describe the flow and its numerical solving method. For ease of description, we first introduce some notations and preliminaries. 3.1. Notations and preliminaries Let M = {x(u1 , u2 ) ∈ R3 , (u1 , u2 ) ∈ Ω ⊂ R2 } be a piece of regular smooth parametric surface. The tangent plane and the unit normal vector at x of M are respectively defined as Tx M = span{xu1 , xu2 } and

n=

xu1 × xu2 , ‖xu1 × xu2 ‖

where xu1 and xu2 are the coordinate tangent vectors, and × denotes the usual cross product operation of two vectors in Euclidean space R3 . The first fundamental form of M is I = ⟨dx, dx⟩ = g11 du1 du1 + 2g12 du1 du2 + g22 du2 du2 with the coefficients gαβ = ⟨xuα , xuβ ⟩, α, β = 1, 2, and ⟨·, ·⟩ denotes the usual scalar product of two vectors in R3 . The second fundamental form of M is II = −⟨dx, dn⟩ = b11 du1 du1 + 2b12 du1 du2 + b22 du2 du2

with the coefficients bαβ = ⟨xuα uβ , n⟩, α, β = 1, 2. Set

[g αβ ] = [gαβ ]−1 .

g = det[gαβ ],

Then the mean curvature of surface M is given as H =

b11 g22 − 2b12 g12 + b22 g11 2g

.

Tangential gradient operator. Let f be a function in the space C 1 (M ) consisting of all C 1 smooth functions on M. Then the tangential gradient operator ∇M acting on f is defined as

∇M f = [xu1 , xu2 ][g αβ ][fu1 , fu2 ]T ∈ R3 , where fu1 = ∂ f (x(u1 , u2 ))/∂ u1 and fu2 = ∂ f (x(u1 , u2 ))/∂ u2 . Divergence operator. Let v be a C 1 smooth vector field on the parametric surface M. Then the divergence of v is defined as 1 divM (v) = √ g

[

]  ∂ ∂ √ αβ T , g [g ][xu1 v, xTu2 v]T . 1 2 ∂u ∂u

Laplace–Beltrami operator. Let f be a function in the space C 2 (M ) consisting of all C 2 smooth functions on M. Then the Laplace–Beltrami operator acting on f is defined as ∆M f = divM (∇M f ). It is well known that ∆M x = 2Hn. Green’s formula. Let M be an orientable surface, and v be a smooth vector field on M. Suppose f ∈ C01 (M ), then

∫

divM (f v)dA = M

∫

[⟨v, ∇M f ⟩ + f divM (v)]dA = 0,

(3.1)

M

where C01 (M ) is a space consisting of all C 1 smooth functions on M with compact support. 3.2. Mean curvature flow In three-dimensional space R3 , starting from a compact immersed orientable surface M0 with its general surface point x ∈ M0 , a family {M (t ) : t ≥ 0} of smooth orientable surfaces evolves according to the mean curvature flow

∂x = 2Hn, ∂t

M (0) = M0 ,

(3.2)

where t is the time-evolving parameter. The surface satisfying H = 0 is called aminimal surface which arises as the critical of the area functional M dA. Since ∆M x = 2Hn, the mean curvature flow can be written as

∂x = ∆M x, ∂t

M (0) = M0 .

(3.3)

The steady solution of this flow is a minimal surface. Since the equation is nonlinear, it is difficult to obtain its analytical solution. To solve the equation numerically, we employ a finite element method based on the extended Loop’s subdivision scheme described above. 3.3. Finite element method for the mean curvature flow Let M (t ) be the limit surface of the extended Loop’s subdivision scheme for the control mesh Md (t ) where t is the time parameter. VM (t ) ⊂ C 1 (M (t )) is a finite-dimensional function space defined by

Q. Pan, G. Xu / Computer-Aided Design 43 (2011) 374–380

the extended Loop’s subdivision scheme for the discrete function values on the vertices. Multiplying a trial function θ ∈ VM (t ) ∩ C01 (M (t )) to both sides of (3.3) and applying Green’s formula (3.1), we obtain the following weak form equation of (3.3).

 Find x(t ) ∈ VM3 (t ) , such that   ] [ ∫  ∂ x(t )  T θ + (∇M x(t )) ∇M θ dA = 0, ∂t M (t )    ∀θ ∈ VM (t ) ∩ C01 (M (t )),   M (0) = M .

(3.4)

Suppose φi is a basis function of VM (t ) corresponding to vertex xi (t ), i = 1, . . . , m. Assume that x1 (t ), . . . , xm0 (t ) are the interior vertices, and xm0 +1 (t ) = xm0 +1 (0), . . . , xm (t ) = xm (0) are the fixed boundary vertices. Then x(t ) could be represented as x(t ) =

m0 −

m −

xi (t )φi +

i =1

xi (t )φi ,

xi (t ) ∈ R3 .

i=m0 +1

Take θ = φj (j = 1, . . . , m0 ), (3.4) can be rewritten as

∫ ∫ − m0 m0 −  ′  (∇M φi )T ∇M φj dA xi (t ) xi (t ) φi φj dA +    M (t ) M (t )  i=1 i =1 ∫ m − (3.5) = − x ( t ) (∇M φi )T ∇M φj dA, j = 1, . . . , m0 ,  i   M ( t )  i=m0 +1   (0) xj (0) = xj , j = 1, . . . , m, where x0j is the jth vertex of the initial mesh Md (0). System (3.5) is a set of nonlinear ordinary differential equations for the unknown vertices xi (t ), i = 1, . . . , m0 . The system is nonlinear because the domain M (t ), over which the integrations are taken, is also unknown. It should be noted that φi depends upon the topology structure of the mesh around xi (t ), but not on the positions of the control vertices. Hence φi is kept the same for any time t. To solve (3.5), we use the forward Euler scheme to discrete x′i (t ) as x′i (t ) ≈

xik+1 − xki

τ

for a given temporal step-size τ , and use a semi-implicit scheme to discretize the remaining terms. We obtain a linear system m0 −

xki +1

i=1

=

m0 − i =1

∫ Mk

xki

φi φj dA + τ

m0 −

xki +1

i=1

∫ Mk

φi φj dA − τ

j = 1 , . . . , m0 ,

m − i=m0 +1

∫ Mk

x0i

(∇M φi )T ∇M φj dA

∫ Mk

377

the basis function associated with the vertex xi , including the corner vertices and the boundary vertices, of a control mesh Md (t ). φi is defined by the limit of the extended Loop’s subdivision scheme with zero function values at every vertex, except for the vertex xi , where it is one. The support of φi is local and covers the 2-ring neighborhoods of vertex xi . We need to evaluate φi and its partial derivatives for forming linear system (3.6). Let ej , j = 1, . . . , mi , be the 2-ring neighborhood elements of xi . If ej is regular, the explicit box spline expression exists for φi on ej . If ei is irregular, local subdivision, as described in Section 2.3, is needed around ei until the parameter values of interest are interior to a regular patch. In the real computation, parameter values are taken to be the Gaussian quadrature knots in the interior of the unit triangle. Only a few subdivision steps are required to get the Gaussian quadrature knots in the interior of a regular triangle. of  The 3limit surface  Loop’s subdivision is expressed as M = ∑ x∈R :x= xi φi . Each triangular surface patch of M is defined locally by only a few related basis functions since the supports of the basis functions are compact. For an interior triangle [xi xj xk ], the related basis functions that define the surface patch over the triangle are uniquely determined by the valences ni , nj and nk of vertices xi , xj and xk , respectively. Hence, interior triangles that have the same valence for each of the three vertices will have the same set related basis functions. To reduce the computation costs of evaluating these functions in the numerical integration, interior triangles are classified into categories according to their vertex valences. All members in one category will have the same vertex valences, hence the same set of related basis functions. The surface evolution in time t does not change the vertex valences of the mesh. Hence, for one category of patches, we only need to evaluate the basis functions once. Using tree structures, the classification could be conducted within a linear time. If a triangle is not interior one, its related basis functions depend on its neighboring extended Loop’s subdivision scheme. Hence, the classification strategy could not be used. In this case, each triangle is classified as a single category. However, since the evolution process does not change the topology of the control mesh, only its vertex positions, the basis functions of the finite element space are fixed. This means that we only need to calculate the basis functions φi once. In fact, what we need is to compute the basis function values and its partial derivatives at the Gaussian quadrature knots. 3.5. Parameterization of subdivision surface and functions on the surface

(∇M φi )T ∇M φj dA, (3.6)

for the unknowns xik+1 , where M k is the limit surface for the control vertices {xki }. The coefficient matrix of system (3.6) is sparse since the basis functions φi have local support (see Section 3.4). System (3.6) is iteratively solved for k = 0, 1, . . . , using the GREMS method till a termination condition is satisfied. The termination condition is specified as max ‖xik+1 − xki ‖ ≤ ϵ i

for a given small value ϵ (we take it as 0.001). There are some previous work on the convergence analysis of the finite element method for the mean curvature flow (see [21–23]). All these works are for graph surfaces (function form surfaces). 3.4. Basis functions The basis functions of the finite element function space VM (t ) are the quartic box splines as mentioned above. We use φi to represent

In Riemannian geometry, differentiable functions are smooth and C ∞ . However, our discretized version of the diffusion problem is in class C 1 . As we mentioned earlier, the functions are defined by the limit of Loop’s subdivision. Such a function is C 2 smooth everywhere except at the extraordinary vertices, where it is C 1 . The function is locally parameterized as the image of the unit triangle defined by T = {(v, w) ∈ R2 : v ≥ 0, w ≥ 0, v + w ≤ 1}, where (1−v−w, v, w) is the barycentric coordinate of the triangle. Using this parameterization, our discretized representation of M is k ˚ ˚ ˚ M = α=1 Tα , Tα ∩ Tβ = ∅ for α ̸= β , where Tα is the interior of the triangular function patch Tα . Each triangular surface patch is assumed to be parameterized locally as xα : T → Tα ; (v, w) → xα (v, w),

(3.7)

where xα (v, w) is defined as an explicit box spline, whether it is a regular or an irregular surface patch, in [16]. Unlike the differentiable structure of a manifold, our parameterization has no overlap. Each point p ∈ M has unique parameter coordinates,

378

Q. Pan, G. Xu / Computer-Aided Design 43 (2011) 374–380

a

a′

a′′

b

b′

b′′

c

c′

c′′

Fig. 4.1. (a), (b) and (c) are the boundary curves of the three models respectively. (a′ ), (b′ ) and (c′ ) are their corresponding initial constructed surfaces. (a′′ ), (b′′ ) and (c′′ ) are their minimal subdivision surfaces.

except at the boundary of the patches. Under this parameterization, the function itself and its partial derivatives can be computed directly. The integration of a function on the surface M is calculated as

∫ f dA := M

−∫ α

 f (xα (v, w)) det(gij ) dv dw,

(3.8)

T

where gij are the coefficients of the first fundamental form of the surface M. The integration on the triangle T is computed adaptively by numerical methods.

as (see [25]) Vt (x)T (x) = I − n(xki )n(xki )T U0 (xki ).





(3.10)

Here I is a 3 × 3 unit matrix, n is the surface normal of the limit surface and

U0 (xki ) =

1

−

(xkj − xki ),

card(N1 (i)) j∈N (i) 1

where N1 (i) is the index set of the 1-ring neighboring vertices of vertex xi , card(N1 (i)) denotes the center of the vertices in N1 (i).

3.6. Mesh regularization 4. Illustrative examples Introducing a mesh regularization mechanism in our evolution process is a necessary step because the surface movement driven by the mean curvature flow may cause a very irregular (nonuniform) distribution of the mesh vertices. Since the tangential displacement does not change the shape of surfaces, but its parameterization (see [24]), we add a tangential movement for the numerical stability of our PDE evolution. The tangential movement can be written as

∂x = Vt (x)T (x), (3.9) ∂t where T (x) is a tangent direction at the control vertex x and Vt (x) is its tangential velocity. The right-hand side Vt (x)T (x) is chosen

In this section, we present a few examples of minimal subdivision surfaces constructed by the proposed method. Some numerical results showing the effectiveness of our construction approach are also given. 4.1. Graphical examples Fig. 4.1 shows three examples of minimal subdivision surfaces with fixed boundaries. We construct the initial surfaces only from the boundary information, where the boundary curves could have discontinuity on its tangent direction, as shown in Fig. 4.1(a).

Q. Pan, G. Xu / Computer-Aided Design 43 (2011) 374–380

379

a

a′

a′′

b

b′

b′′

c

c′

c′′

Fig. 4.2. The first row is the Helicoid surface model: (a) is the boundary curve. (a′ ) is the initial construction. (a′′ ) is its minimal subdivision surface. The second row is the Catenoid surface model: (b) is the boundary curve. (b′ ) is the initial construction. (b′′ ) is its minimal subdivision surface. The third row is the Enneper surface model: (c) is the boundary curve. (c′ ) is the initial construction. (c′′ ) is its minimal subdivision surface.

Table 4.1 k describes the subdivision times, where we subdivide the nine models at six different levels respectively. The data from the second to the tenth rows are the maximum approximate errors of the mean curvature computed from the discrete solutions of the PDE evolution. Asymptotic maximum errors Models

k=0

k=1

k=2

k=3

k=4

k=5

Fig. 4.1(a) Fig. 4.1(b) Fig. 4.1(c) Fig. 4.2(a) Fig. 4.2(b) Fig. 4.2(c) Fig. 4.3(a) Fig. 4.3(b) Fig. 4.3(c)

1.225E−1 5.882E−1 1.361E−1 9.981E−2 2.218E−1 3.702E−2 1.670E−1 1.140E−1 7.260E−2

5.941E−2 3.700E−1 7.227E−2 5.345E−2 5.501E−2 1.553E−2 8.159E−2 5.601E−2 3.385E−2

2.945E−2 2.344E−1 4.238E−2 2.934E−2 1.373E−2 6.441E−3 4.030E−2 2.782E−2 1.664E−2

1.586E−2 1.547E−1 2.546E−2 1.773E−2 3.181E−3 2.584E−3 2.001E−2 1.359E−2 8.057E−3

9.616E−3 1.082E−1 1.546E−2 1.163E−2 6.394E−4 1.018E−3 1.006E−2 6.844E−3 4.054E−3

6.587E−3 9.678E−2 9.837E−3 7.916E−3 1.580E−4 4.089E−4 5.161E−3 3.690E−3 2.202E−3

Helicoid, Catenoid and Enneper are the minimal surfaces with analytic forms. In Fig. 4.2, we discretize these three analytic surfaces at a rough level, linearly refine them several times, and make them as their initial constructions. We achieve the minimal subdivision surfaces as the steady solutions of (3.4). It is clearly seen that the initial surfaces are very different from the final minimal subdivision surfaces. We also construct three practical models in Fig. 4.3, such as umbrella, roof models of pavilion and Chinese palace. Their initial surfaces are also constructed only from the boundary information at a rough level, and they are linearly refined several times. The resulting surfaces fully reflect the beauty of minimal surfaces.

4.2. Refinement and convergence In order to show that the proposed method is effective, we list the maximum errors of the mean curvature H computed from the discrete solutions of our numerical method for the nine models. The asymptotic errors are presented in Table 4.1. In this table, we construct their initial surfaces as the initial value of the PDE evolution through subdividing the nine models at gradually more and more dense levels respectively according to the extended Loop’s subdivision scheme. From the numerical results, we can see that the maximum of the approximation errors of the mean curvature monotonously decline with the increasing of subdivision times k. Hence, our numerical method is convergent.

380

Q. Pan, G. Xu / Computer-Aided Design 43 (2011) 374–380

a

a′

a′′

b

b′

b′′

c

c′

c′′

Fig. 4.3. The first row is a roof model of the pavilion with four corners: (a) is the boundary curve. (a′ ) is the initial construction. (a′′ ) is the resulting minimal subdivision surface. The second row is an umbrella model: (b) is the boundary curve. (b′ ) is the initial construction. (b′′ ) is the resulting minimal subdivision surface. The third row is a roof model of the Chinese palace: (c) is the boundary curve. (c′ ) is the initial construction. (c′′ ) is the resulting minimal subdivision surface.

5. Conclusions There has been much rich work on minimal surfaces. Subdivision technology can provide a simple and efficient method to describe freeform surfaces with arbitrary topology, at the same time satisfy some smoothness requirement. In this paper, we constructed efficiently minimal subdivision surfaces with given boundaries through the mean curvature flow. For the numerical solution of the flow, we have adopted the finite element method based on the extended Loop’s subdivision scheme. The experiment results showed that our method is convergent and yields high quality minimal surface. Acknowledgements The first author was supported in part by NSFC grant 10701071 and Program for Excellent Talents in Hunan Normal University (No. ET10901). The second author was supported in part by NSFC grant 60773165 and NSFC key project under the grant 10990013.

[8] [9] [10]

[11]

[12]

[13] [14] [15] [16] [17] [18]

References [1] Monterde J. Bézier surface of minimal area: the Dirichlet approach. Computer Aided Geometric Design 2004;21:117–36. [2] Cosin C, Monterde J. Bézier surfaces of minimal area. In: Proceedings of the int. conf. of computational science. Lecture notes in computer science, vol. 2330. Amsterdam: Springer-Verlag; 2002. p. 72–81. [3] Arnal A, Lluch A, Monterde J. Triangular Bézier surfaces of minimal area. Lecture notes in computer science, vol. 2669. 2003. [4] Jin W, Wang G. Geometric modeling using minimal surfaces. Chinese Journal of Computers 1999;22(12):1276–80. [5] Man J, Wang G. Approximating to nonparameterzied minimal surface with Bspline surface. Journal of Software 2003;14(4):824–9. [6] Man J, Wang G. Minimal surface modeling using finite element method. Chinese Journal of Computers 2003;26(4):507–10. [7] Polthier K. Computational aspects of discrete minimal surfaces. In: J. Hass, D. Hoffman, A. Jaffe, H. Rosenberg, R. Schoen and M. Wolf (Eds.), Proc.

[19] [20] [21] [22]

[23] [24]

[25]

of the clay summer school on global theory of minimal surfaces, 2002. citeseer.ist.psu.edu/polthier02computational.html. Xu G. Geometric partial differential equation methods in computational geometry. Beijing (China): Science Press; 2008. Biermann H, Levin A, Zorin D. Piecewise-smooth subdivision surfaces with normal control. In: SIGGRAPH, 2000, pp. 113–120. Hoppe H, DeRose T, Duchamp T, Halstend M, Jin H, McDonald J, Schweitzer J, Stuetzle W. Piecewise smooth surfaces reconstruction. In: Computer graphics proceedings, annual conference series, ACM SIGGRAPH94, 1994, pp. 295–302. Hoppe H, DeRose T, Duchampi T, McDonald J, Stuetzle W. Mesh optimization. In: Computer graphics proceedings, annual conference series, SIGGRAPH93, 1993, pp. 19–26. Cirak F, Ortiz M, Schröder P. Subdivision surfaces: a new paradigm for thinshell finite-element analysis. International Journal for Numerical Methods in Engineering 2000;47:2039–72. Kobbelt L. A variational approach to subdivision. Computer Aided Geometric Design 1996;13(8):743–61. Kobbelt L, Schröder P. A mutiresolution framework for variational subdivision. ACM Transactions on Computer Graphics 1998;17(4):209–37. Diewald U, Morigi S, Rumpf M. A cascadic geometric filtering approach to subdivision. Computer Aided Geometric Design 2002;19(9):675–94. Loop C. Smooth subdivision surfaces based on triangles. Master’s thesis. Technical report. Department of Mathematices, University of Utah. 1978. Stam J. Fast evaluation of Loop triangular subdivision surfaces at arbitrary parameter values. In: SIGGRAPH ’98 proceedings. 1998 [CD-ROM supplement]. Bajaj C, Xu G. Anisotropic diffusion of subdivision surfaces and functions on surfaces. ACM Transactions on Graphics 2003;22(1):4–32. Xu G, Pan Q. G1 surface modelling using fourth order geometric flows. Computer-Aided Design 2006;38(4):392–403. Zorin D, Kristjansson D. Evaluation of piecewise smooth subdivision surfaces. The Visual Computer 2002;18(5–6):299–315. Deckelnick K, Dziuk G. A fully discrete numerical scheme for weighted mean curvature flow. Numerische Mathematik 2002;91:423–52. Deckelnick K, Dziuk G. Mean curvature flow and related topics. In: Blowey JF, et al., editors. Frontiers in numerical analysis. Berlin: Springer; 2003. p. 63–108. 10th LMS-EPSRC numerical analysis summer school, Durham, UK, July 7–19, 2002. Huisken G. Flow by mean curvature of convex surfaces into spheres. Journal of Differential Geometry 1984;20:237–66. Epstein CL, Gage M. The curve shortening flow. In: Chorin A, Majda A, editors. Wave motion: theory, modeling, and computation. New York: SpringerVerlag; 1987. p. 15–59. Xu G, Pan Q, Bajaj C. Discrete surface modelling using partial differential equations. Computer Aided Geometric Design 2006;23(2):125–45.