ICES Report 03-51

SURFACES RECONSTRUCTING OF WITH BIQUARTIC PATCHES AND OTHER TECHNIQUES Dong Xue, Leszek Demkowicz Texas Institute for Computational and Engineering Science The University of Texas at Austin Austin, TX 78712

Waldek Rachowicz Cracow University of Technology, Poland Abstract

We present in this paper an approach to reconstruct a surface that interpolates given surface is represented in a parametric form vertices with prescribed normal vectors. The implemented in two main steps: construction of a curve network using Hermite curves, followed by a construction of corresponding cubic normals, and surface fitting using biquartic rectangular patches interpolation. The proposed scheme is based on the idea of - and compatibility conditions for the curve net and the rectangular patches, respectively.

1

Key words: surface, rectangular patch, Hermite curve, interpolation Cross Boundary Derivatives(CDBs).

Acknowledgment The work has been done in collaboration with Zhang Yongjie from the Computational Visualization Center (CCV). All computations presented in the work were done through the National Science Foundation’s National Partnership for Advanced Computational Infrastructure (NPACI).

Contents 1 Introduction

4

1.1

Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.2

The Report Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2 Analysis

6

2.1

The

Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.2

Why cubic curve? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.3

Why cubic normal? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3 Biquartic 3.1

13

Curve Network Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.1.1

Hermite Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.1.2

Cubic Normals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3.2

3.3

Surface

surface fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.2.1

Vertex Nodes and Mid-Edge Nodes Contributions . . . . . . . . . . . . . . 19

3.2.2

Middle node Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Modeling of a Human Head . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Other

Interpolation Methods

32

4.1

Inverse Distance Weighting (IDW) Interpolation . . . . . . . . . . . . . . . . . . . 33

4.2

Bicubic

Surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2

5 Conclusions and Future Works

43

A Singular Value Decomposition (SVD)

44

B Subroutines

49

3

1 Introduction In a smooth surface reconstruction, two main forms are widely used to interpolate a parametric representations and implicit representations.

surface:

There are various successful approaches of using implicit representation in interpolating of a surface. Bajaj and Xu pioneered implicit algebraic surface patches, A-patches [1, 2, 3, 4] (see also [5],[6]). A-patches are defined by a manifold triangulation in . This technique is based on the idea of zero contour of trivariate function. Other implicit schemes include B-patches [7], and Spatches [8] which also provide elegant solutions to the smoothing problem at the cost of a currently non standard patch representation.

The proposed surface reconstruction scheme is based on an explicit parametric representation of rectangular patches focusing an arbitrary unstructured surface mesh. The scheme has been implemented within our Geometric Modeling Package (GMP), fitting into a general class of both explicit and implicit parameterizations. The GMP is used to support geometrical representation for Finite Element (FE) discretization. Compared with implicit parameterizations, parametric representations have the advantage of simplicity and better hold on the effect of the geometric data (in our case -an unstructured surface grid of bilinear quadrilaterals with normals prescribed at the grid points) on the ultimate shape of the reconstructed surface.

A number of methods for constructing smooth parametric surfaces has been proposed. The splines assembled from B-splines are widely used, because they have explicit formulas and built-in continuity. However the tensor product B-spline representation has a major shorting, the B-spline mesh must be a structured mesh. Using non uniform rational B-splines (NURBS) does not solve this problem since the trimming destroys one of the chief advatages of the B-spline representation, its built-in smoothness so that one ends up with the tricky task of smoothly joining the trimmed pieces [9, 10]. J. Peter pioneered the idea of surface splines [11, 12, 13, 14] (see also[15, 16]), which devises a representation that removes the regularity restrictions at the cost of creating a refined mesh of quardrilateral subcells. The mesh refinement can separate irregular vertices by implementing initial mesh refinement, edge cutting and quadratic meshing. The refined mesh guarantees that each orginal vertex is surrounded by vertices of degree four. The approach is then reduced to the B-spline paradigm by constructing surface parametrization using quadratic patches. The scheme is in the same conceptual frame work as tensor-product B-splines, and local evaluates by averaging and obeys the convex hull property [17, 18]. Other works include generalized subdivision scheme, [19, 20] and reparametrization scheme [21, 22].

4

1.1 Motivation The -adaptive FE method is one of the most powerful methodologies for simulating complex engineering problems, especially those involving singular solutions, complicated geometry, or multiple media with varying properties. The present work is motivated by an application of methods to simulate the absorption and diffraction of EM-waves in the human head. The EM waves are emitted by a dipole antenna which models a mobile phone. The issue of a precise geometric modeling is especially sensitive in high accuracy simulations. In our previous work, we interfaced GMP with the geometric data provided by J. Zhang at CCV, to obtain the connectivity information, and construct a piecewise trilinear model of the human head [23]. The data provided by Zhang is generated from MRI scans; the head model is represented in terms of the following nodal information: Nodal connectivities for each hexahedron. Coordinates of each node in a physical frame. Normal vectors

for each node on the surface of the head.

Obviously, a linear geometry model is insufficient for high accuracy simulations. The hp-adaptive method requires at least a continuous geometry reconstruction. If we only use continuous geometry reconstruction, the FE code adapts meshes to resolve the non-physical scattering of waves on edges resulting from the poor, low order geometry representation. The term interpolation is used here to describe ways of fitting a curve or surface to a set of data points or curves [24]. The objective of this project is to reconstruct a surface that interpolates given vertices with prescribed normal vectors. The proposed surface reconstruction scheme is based on a parametrized surface by implementing a local interpolation on rectangular patches, and it has the following properties: local control of shape, numerical stability, smoothness and continuity, ability to evaluate derivatives. The reconstruction of free-form geometry is the main modeling challenge. Generally speaking, there is no restriction on the number of cells meeting at a mesh point or the number of edges adjacent to a mesh cell. The surface Mesh cells need not be planar. 5

1.2 The Report Structure This paper is structured as follows: Section 1 explains the motivation for developing the surface reconstruction schemes. Section 2 lists the main steps for the surface reconstruction scheme, establishes the compatibility conditions, and explains the reason for using Hermite curve and cubic normals. Section 3 elaborates on the processes of constructing a curve network and biquartic rectangular patches. Several examples are given to demonstrate the effectiveness of the scheme. Section 4 presents details on implementing two other alternative surface interpolation schemes: Inverse Distance Weighting (IDW) 4.1 and Bicubic surface Interpolation 4.2. Section 5 summarizes and Appendix A illustrates the Singular Value Decomposition (SVD) used in section 3. Appendix B lists routines developed in the course of this project.

2 Analysis 2.1 The Constraints

Consider a rectangular patch with parameters and its neighbor patch with parameters as illustrated in Fig. 1. Concentrate on the common boundary where . If both rectangular patches are sufficiently smooth, and are compatible if and only if the surface normal to is well defined and agrees with the normal of at each point of the boundary:

(2.1) guarantees non-vanishing normals. For surfaces, continuity means oriented tangent plane continuity [25], smoothness for surfaces differs from smoothness for

bivariate functions. In the proposed methodology, the two steps,

!#"$ %'&)(+* ,.-0/ !213"

surface reconstruction will be done in

Step 1: Construction of a curve network. We construct admissible parametric curves , which are continuous at the vertices. Continuity is achieved by using an alternative sufficient constraint that forces the mesh curves to interpolate vertex specified at each data while having compatible normals point on the boundary of a rectangular patch, i.e., along the mesh curves, see Fig.2.

! "4$ %'&5(6* 78,.-0/.9 1:"

Step 2: surface fitting. We interpolate between curves of the network obtained in Step 1 using a smooth parametrized surface, by implementing an algorithm for local interpolation of rectangular patches, see Fig.3. A rectangular patch is the image of a bivariate polynomial 6

X1( u , v ) v

v

1

u

w

0

1

X 2 ( u, w )

u u 1

0

1

w

Figure 1: Parametrization of two abutting rectangular patches

N( ξ )

n1

n2

X’( ξ ) p2

X( ξ )

p1

, and its normal vector , ! ! 1 " $ % '&)( * parametrized with parameters 9 / restricted to a standard domain, 9 " Figure 2: A parametric curve

.

The compatibility conditions for curves should guarantee the mity to the boundaries. They consist of two parts,

smoothness as well as confor-

compatibility conditions for curves:

, / 11 , / 11 , / & , /0&

(2.2)

compatibility conditions for normals:

/ !21 , / ! 1 7, / ! 1 7

(2.3)

In a similar way, the compatibility conditions for rectangular patches also consist of two parts,

/0& 1 1 /

compatibility conditions for patches:

/ ! ! 1 1 , )/ ! ! 1 1

/ 5& , '/

(2.4)

compatibility conditions for patches:

! 1 ! 1 ! / 7 , )/

1/ ! / ! '& 1 7 , / ! 1

/0& 1

/ 1 / 1

7, / ! 1

/ 1 1 0/ & ! ! 1 1 , '/ ! ! 1 1

/ / % , /

where is the normal function along the curve; tives (CBDs).

! 1 ! 1 ! / & 7 , '/

! 1 ! 1 ! /% 7 , / % - & are the Cross Boundary Deriva

ξ2 1

X ( ξ 1, ξ 2 ) 1

0

ξ1

n3 n4

N c3 ( ξ 1 ) N c4 ( ξ 2 ) (ξ ) X n 2 p4 c3 1 n X c4 ( ξ 2 ) 1 p 3 N c2 ( ξ 2 ) Xc2 ( ξ 2 )

p1

N c1 ( ξ 1 )

X c1( ξ 1 )

p2

Figure 3: Interpolation of a rectangular patch

2.2 Why cubic curve? We usually want the reconstructed curve to be as smooth as possible, i.e., minimize the wiggles. avoid using high-degree polynomials.

8

Cubic curves are commonly used in graphics because curves of lower order commonly have too little flexibility, while curves of higher order are usually considered unnecessarily complex and can easily introduce undesired wiggles. The use of cubic curve can also be explained in an optimization sense. As the second derivative of curve approximates curvature, we determine an optimal curve reconstruction scheme by seeking the solution to the following variational problem:

, / !21

,

Given endpoints and normal vector , find a curve parametrization which satisfies compatibility conditions in equation (2.2) and minimizes the mean (linearized) curvature,

$ %'&)( *

, / ! 1 ! * (2.5) / %'& 1 , and it admits a unique solution which The problem can be stated formally in space

satisfies the following (equivalent to the minimization problem) variational statement,

(2.6) , , ! % for every parametrization, test function , satisfying homogeneous essential boundary condi-

tions,

, / 1 1 , / 11 , /0& , / & 6

Integration by parts leads to

,

,

!

(2.7)

, , ! , , , ! , , , , ,

(2.8)

Restricting ourselves first to test functions that vanish on the boundary along with their first order derivatives, we get

, , ! This implies that , and, therefore, ,

variational statement reduces to,

, "

/ %'& 1

must be a cubic polynomial

, / & 1 , / & 1 , / 1 , / 1 9

(2.9)

, " .

The

, / & 1 , , / 1 , we conclude that (2.10) , / & 1 , / & 1 !21 . The direction of unit binorThe unit tangent vector on the curve is defined as /

! 1 !21 / ! 1 / ! 1 . mal is determined by / . The unit principal normal / 1 1 Decomposing , / & and , /0& into their tangent and principal normal components, , /0& 11 / , /0& 1 1 1 1 / , /0& 1 1 1 1 , /0& / , /0& / , / & (2.11) Assuming

we reduce equation(2.10) to,

/ , /0& 1 1 / , /0& 1 1 / , /0& 1 1 / , / & 1 1 (2.12) This, together with essential conditions for test function in equation (2.8), implies the final natural ! boundary condition at & , , /0& 1 % 1 or equivalently, , / & 6 . In the same way, , / 1 , / 1 (2.13)

The cubic curve reconstruction has the following properties: it minimizes the mean (linearized) curvature, it provides the lowest order polynomials representation that interpolates two points and allows for the gradient at each point to be defined which makes the continuity possible it reduces automatically to linear polynomials ( a straight line segment), for appropriate vertex data. The algebraic form of a parametric cubic curve is given by the following vector equation:

! , / !21 ! !

(2.14)

where are vector algebraic coefficients. The coefficients are not the most convenient way of controlling the shape of the curve in typical modeling situations, nor do they contribute much to an intuitive understanding of the curve. A practical alternative is offered by the Hermite 10

n 2(ξ )

n 1(ξ )

t1 ( ξ )

t 2 (ξ ) p2

p1 Figure 4: A Hermite curve

interpolation, which allows us to define a curve segment in terms of its endpoints [26]. More specifically, we need to know the coordinates and tangent vectors at both endpoints, see , we obtain the algebraic coefficients Fig.4. Evaluating equation(2.14) and its derivative at in terms of the boundary conditions by solving a set of four equations in four unknowns,

! %'&

(2.15) Substituting the coefficients into equation (2.14) and rearranging the terms produces, , / !21 / ! ! ! ! & !21 1 /! ! ! 1 ! 1 / !21 ! 1 / !21 !21 (2.16) / / / /

where

!21 ! 1 ! 1 ! / ! 1 / ! / !21 / /

! / !21 !21 ! / / 2! 1

are the Hermite Basis Functions. Differentiating equation(2.16), we get,

, / 2! 1

(2.17)

Matrix algebra and its notation scheme offer a compact mathematical form for representing a curve. We can rewrite equation (2.16) as a product of two matrices. in matrix form, where

is a

, / ! 1

$ ! ! ! 5& ( and

& & &

Hermite basis transformation matrix,

&

(2.18)

&

$ (

(2.19)

is the Hermite control matrix,

Equation (2.18) is also known as the cubic Hermite spline equation. 11

(2.20)

2.3 Why cubic normal?

Figure 5: Five basis functions and their first derivatives To obtain the desired conditions in equation (2.5), we need to find the normal derivative along the mesh curves. The approach does away with position vectors and tangent vectors . In other words, the curve reconstruction and , and uses an independent normal function defined along involves not only defining the curve but also specifying the variation of normals the curve and matching the normals at the endpoints,

7 , /! 1

78, / 1

7 ,

(2.21) 7 ,/&1

! ! 1 The parametrization for a rectangular patch 4/ is then requested to satisfy compatibility conditions (2.4) and compatibility conditions (2.5). The use of bicubic polynomials " which in turn led to geometrically unacceptable situation. implied 7, ! ! Biquartic functions involve linear combinations of twenty five monomials in and . The scalar-valued interpolation involves tensor products of five one dimensional shape functions, the !1 !21 Hermite basis / + & '' and a fifth bubble function / , illustrated in Fig.5. The order of the surface interpolation introduces four additional edge shape functions which result in a total 1 scalar unknowns to satisfy condition. Thus, we have five coefficients per of / edge. As Condition (4.84) must be enforced along each edge, the maximum polynomial order of 7 , / !21 is three, 7, " .

12

3 Biquartic

&

Surface

3.1 Curve Network Construction

, / ! 1

78, / ! 1

The curve network construction includes Hermite curves interpolation along with its cubic normal interpolation . The curve network should satisfy compatibility conditions in compatibility conditions in (2.3) for normals. The curve net will later be (2.2) for curves and extended to a smooth surface. 3.1.1 Hermite Curves

of polynomial (2.16). A Hermite curve is constructed by determining the vector coefficients The two unknown vectors can be expressed in terms of their Euclidean components, and . The boundary conditions resulted from the variational principle are,

/

1

/

1

essential boundary conditions:

, / 11

, / & %

natural boundary conditions:

, / 11 , /0&

(3.22)

(3.23)

where are two unknown scalars. equations (3.22) and (3.23) yield a linear system of equations to be solved for the components of , , and constants and .

, / !21"

The cubic curve may degenerate to a lower order polynomial. In the case of a straight line segment, , the vector coefficients corresponding to the third order and the second order terms vanish,

(3.24)

As the second derivatives of the curve vanishes,

, / !21 6

13

(3.25)

System (3.23) can now be solved for

and

,

which yields the final form for the curve,

/ ! 1 /0& ! 1

In this degenerated case, the nodal data of the curve bility conditions,

(3.26)

! , ,

(3.27) and

must satisfy the compati-

Assume now that the cubic curve degenerates to a second order polynomial vector coefficient corresponding to the third order term vanishes,

, / !21 "

(3.28)

, the

(3.29)

The second derivative of the curve is a constant different from zero, which implies that,

, / 1 , /0& 1

We have

/

(3.30)

1 /

1

(3.31)

This implies condition (3.28) from which, in turn, follows that the curve must degenerate to a line segment. Consequently, the cubic can never degenerate to a second order polynomial. 3.1.2 Cubic Normals The third order normals along the curve can be written as:

(3.32) / ! 1 / ! 1 / ! 1

!21 !21 where are two unknown vector coefficients, / and / are usual linear shape functions !21 and / + % are integrated Legendre polynomials. Note that, (3.33)

78, / ! 1 / !21

14

They can be expressed as,

! & !! ! /! ! /

with derivatives,

&

! /! & & 1 & ! & !

& 11 ! 1 & / &

! &

(3.34)

! 1 7 , / ! 1 is a fifth order !21 & , / points. As / / ! 1 is equivalent to enforcing 1 / & 1 % / (3.35) / 1 / 1 / 1 / 1 1 / & 1 have already been satisfied in (3.22). Taking / 1 , we have, Conditions / / 1 7, / 1 1 , / 1 1 / / / , / 1 / , / 1 1 / 1 1 Similarly, taking / , we have, 1 7 , / 1 , / 1 7 , / 1 , / 1 / / / 1 1 1 / 1 1 , / 1 / 1 1 / 1 1 / 1 , / 1 / , / / / , / / 1 Then, taking / , we have, 1 7 , / 1 , 1 1/ 1 1 27 , / 1 , / 1 1 1 1 / / , / / / , / / 1 Finally, taking / , we have, 1 1 , / 1 / & 27 , / & / , / 1 1 / 1

An degree polynomial fits a curve to polynomial, enforcing condition (2.3)

15

The above equations give a linear system of four equations in terms of six unknown components of . In a matrix form,

where

:

are

(3.36)

matrices of coefficients,

, / 11 , / 1 , /

and

, / 11 2 , / 1 , /

6

(3.37)

is a known right-hand side vector,

1 , / 1 / 1 1 1 1 / , / / , / Two additional needed conditions are constructed by minimizing the norm of the derivative 7 , , !1 ! 7 , / * (3.38)

Using the orthogonal property of integrated Legendre polynomials, we get

(3.39) / 1 / 1 / 1 and / 1 ! ; is subject to the constraints where (3.36). To solve the minimization problem by the method of Lagrange multipliers, we form a

system of ten equations in terms of the components of the two unknown vectors , and four 1 Lagrange multipliers 4/ ' . The matrix equation is (3.40)

where is a & : & square matrix of coefficients; is the solution vector of ten unknowns; is a known right-hand side vector; are two : diagonal matrices,

(3.41) 6

16

& : &'

The matrix may degenerate to a singular matrix. This happens, in particular, when the general cubic curve degenerates to a straight line segment, as discussed in the previous section. The last two constraints in (3.36) then vanish which results in singular matrices , , and the global matrix .

It is for this reason, that we can not employ the standard Gauss elimination to solve system (3.40), and use the Singular Value Decomposition (SVD) technique instead. The implementation based on SVD is uniformly stable for all possible situations including the degenerated and nearly degenerated matrices. Details on the SVD implementation with illustrating numerical experiments are provided in Appendix A.

1 /0& '& 1 / & % 1 /0& '& '& 1 1 /0& %'& 1 / / 3 / 1 / 3 ! 1 1 , / ! 1 7 ,/

To illustrated the resulted curve net, we show a special example of a spherical rectangular patch in Fig.6. The four vertices has position vectors

and the corresponding normal vectors are

. Applying the two steps discussed above, we obtain the curve net with a Hermite curve illustrated in blue and cubic normals along the curve illustrated in red, see Fig.6. The curve net will be used later to construct a rectangular patch.

3.2

surface fitting

,.- & 5' / ! ! 16"

Once we have constructed the four edge functions , along with the corresponding cubic normal vectors , we proceed now with a construction of parameterizations can be written for the rectangular patches. A general biquartic rectangular patch as the sum of vertex nodes contributions , mid-edge nodes contributions , and middle node contribution ,

78,.- & ''

/ ! ! 1 / ! ! 1 5/ ! ! 1 /! ! 1

(3.42)

The biquartic parameterization is constructed by enforcing the Compatibility Conditions (2.4) and the Compatibility Conditions (2.5) on the boundaries. The in expression (3.42) can be written as,

5/ ! ! 1

)/ ! ! 1 )/ ! ! 1

where is a vector coefficient for the middle node contribution, and shape functions, see Fig.15,

/! 1 /! 1 17

(3.43)

is the corresponding face

(3.44)

1.4 1.2

n1 = (0.5, 0.5, −0.5) p3 = (1, 1, 1 )

n = (0.5, −0.5, 0.5) 1 p = (1, 0, 1 )

1

4

0.8 0.6

Z

0.4 0.2 0 −0.2

p1 = (1, 0, 0 ) n1 = (0.5, −0.5, −0.5)

−0.4 0.5

p = (1, 1, 0 ) n22 = (0.5, 0.5, −0.5)

x1 1.5

0.2

0

−0.2

−0.4

y

0.4

0.6

1.2

1

0.8

p = (1, 1, 1 ) n31 = (0.5, 0.5, −0.5)

1.4 1.2

n1 = (0.5, −0.5, 0.5) p4 = (1, 0, 1 )

1

Z0.8 0.6

p2 = (1, 1, 0 ) n2 = (0.5, 0.5, −0.5)

0.4 0.2 0

p = (1, 0, 0 ) n11 = (0.5, −0.5, −0.5)

−0.2

1.5

−0.4 0.8

1

0.9

x

1

0.5 1.1

1.2

0 1.3

y

−0.5

Figure 6: A spherical rectangular patch

! ! 1 ) / /! ! 1

7, -

Note that the term vanishes along all four edges. In other words, does not affect the behavior of on the boundary. On the other side, the first two contributions in (3.42) are uniquely determined by the boundary data - edge functions and normals . This results in a two step procedure:

,.-

18

−3

x 10 4 3.5 3

ψs

2.5 2 1.5 1 0.5 0 1 0.8

1 0.6

ξ

2

0.8 0.6

0.4 0.4

0.2

ξ

0.2 0

1

0

Figure 7: Face shape function

/ ! ! 1 / ! ! 1 5/ ! ! 1

Step1: Construct a surface parameterization lating the boundary data, Step2: Determine vector

interpo-

in (3.43) using a minimum energy principle.

3.2.1 Vertex Nodes and Mid-Edge Nodes Contributions

n3 n4

n t p t n p t p p t 1

2 4

1 4

4

2

3

2 3

ξ2

ξ1

1

t

2 1

t

2

1 1

1 2

t22

Figure 8: Vertex Node Interpolation

/ ! ! 1

t13

"

The in (3.42) is a standard bicubic Hermite surface interpolant involves only vertex data, see Fig.8. The total vertex contribution can be expressed as,

/ ! ! 1

! ! 1 ! ! 1 ! ! 1 / / /

19

.

It

Here

denote the position vectors for each of the four vertices . are the corresponding bicubic vertex shape functions, listed in Table 1 and illustrated in Fig.9.

Vertex number

Table 1: The shape functions for position vector at vertex

!

& ''

&

is the vertex number, and are eight tangent vectors, where is the index of . For any rectangular patch, the eight tangent vectors are obtained from the curve functions ,

, , / 1 , / 1

, 0/ & 1 , / 1

, / & 1 , / & 1

, / 1 , / & 1

!

(3.45)

are the corresponding shape functions for tangent vectors at vertex in terms of , see Figures 10, 11,12 and 13. The eight tangent vector shape functions are listed in Table 2. Vertex number

Table 2: The shape functions for tangent vectors at vertex

are the four unknown twist vectors (mixed derivatives) at each vertex .

are the bicubic shape corresponding to each twist vectors, see Fig.14. The four twist vector shape functions are listed in Table 3. The

5/ ! ! 1

! ! 1 ! ! 1 5/ /

in (3.42) denotes the mid-edge nodes contributions,

20

(3.46)

1 0.8

1

ψ v1 0.6

0.8

0.4

ψ v2

0.2

0.6

0 1

0.4 0.8

0.2 0.6 0 1

ξ

0.4

2

0.8 0.2

0

0.1

0

0.3

0.2

0.4

0.6

0.5

0.7

0.9

0.8

1

0.8 0.6

0.4 0.4

0.2

ξ1

0.2 0

1

ξ1

0

1

0.8

ψ v3

1 0.6

ξ2

0.8

ψ v4

0.6

0.6

0.4

0.4

0.2

0.2

0 1

0 1 0.8

0.8

1

ξ

2

0.6

0.8 0.4

0.6

0

0.4

0.2

0.2 0

ξ

0.2 0

0

1

Figure 9: Vertex shape functions Vertex number

,

0.8 0.6

0.4

2

0.4

0.2

1

ξ

0.6

ξ

1

Table 3: The shape functions for twist vector at vertex

where are eight vector coefficients to be determined, two per edge; shape functions, see Table 4.

are the corresponding

The Compatibility Conditions (2.4) and the fact that the curves have been reconstructed using cubic polynomials only, , imply that contributions corresponding to last four shape functions . Note that the condition must simply vanish, i.e.,

/! ! 1

, / ! 1 "

8% & 5'

21

0.16

0.16

0.14

0.14

0.12

ψt11

0.12

ψt12

0.1

0.1

0.08

0.08

0.06

0.06

0.04

0.04

0.02

0.02

0 1

0 1 0.8

0.8

1 0.6

ξ

1 0.6

0.8 0.4

0.2

2

0.8

0.6

0.4 0.2 0

0

ξ

ξ1

0.6

0.4 0.4

0.2

2

0.2 0

0

ξ1

Figure 10: Shape functions for tangent vectors at vertex 1

0

0.16

−0.02

0.14

−0.04

ψt21

0.12

ψt22

−0.06

0.1

−0.08

0.08

−0.1

0.06

−0.12

0.04

−0.14

0.02

−0.16 1

0 1 0.8

0.8

1 0.6

ξ2

1 0.6

0.8 0.6

0.4 0.4

0.2

0.2 0

0

ξ2

ξ1

0.8 0.6

0.4 0.4

0.2

0.2 0

0

ξ1

Figure 11: Shape functions for tangent vectors at vertex 2

/! ! 1

does not apply to the contributions of the first four shape functions which automatically vanish on the patch boundary and contribute only with non-zero normals, see Fig.15.

The Compatibility Conditions (2.4) require the knowledge of Cross Boundary Derivatives (CBDs) along the boundary. The CBDs have a crucial effect on the shape of the constructed patches; they allow for the patches to effectively reflect the variation of the normals . The CBD at any point on the patch boundary is perpendicular to both the corresponding normal and the tangent vectors. Given the reference coordinates of a point on the rectangular patch, we first identity four corresponding points on the patch edges with coordinates . We compute then the four CBD functions corresponding to each of

7 ,.-

/! ! 1

/ ! 1 /0& ! 1 / ! 5& 1 / % ! 1

22

0

0

ψt31

−0.05

ψt32

−0.05

−0.1

−0.1

−0.15

−0.15

−0.2 1

−0.2 1 1

0.8

1

0.8 0.8

0.6

0.8 0.6

0.6 0.4

ξ

0.2

2

0.6 0.4

0.4

ξ1

0.2 0

0.4

ξ

0.2

2

ξ1

0.2 0

0

0

Figure 12: Shape functions for tangent vectors at vertex 3

0.16 0 0.14 −0.02 0.12 −0.04

ψt41

0.1

ψt42

0.08

−0.06 −0.08

0.06 −0.1 0.04 −0.12 0.02 −0.14 0 1

−0.16 1

0.8

0.8 0.6

1 0.6

ξ2

1 0.6

0.8

0.4 0.2

0.4 0

0.2

ξ2

ξ1

0.8 0.6

0.4 0.4

0.2

0.2 0

0

0

ξ1

Figure 13: Shape functions for tangent vectors at vertex 4 the four edges of the patch:

/! 1 /! 1 / & ! 1 / ! 1 / ! '& 1 / ! 1 / % ! 1 / ! 1

/ !! 11 / ! 1 / ! 1 /

/! 1 /! 1 /! 1 /! 1 /! 1 /! 1 / ! 1 / ! 1

78, "

! 1 / !1

/! 1

/! 1

/

" -

(3.47)

Note that CBDs are fourth order polynomials along the edges and the normals along the curve are third order polynomials . In order to satisfy the compatibility conditions, the rectangular patch must satisfy the system of four equations in (2.5). With , the 23

! ! !

0.025

0

0.02

ψc2

ψc1 0.015

−0.005 −0.01 −0.015

0.01

−0.02

0.005

−0.025 1 0.9 0.8

0 1

1

0.7 0.6 0.8

0.5

1 0.6

ξ

0.8

0.4

0.2

2

0.2 0

0.3

ξ2

ξ

1

0

0.6

0.4

0.8 0.6

0.4

0.4 0.2

ξ

0.2

0.1

1

0

0

0.025

0.02

0

ψc4

ψc3 0.015

−0.005 −0.01 −0.015

0.01

−0.02 −0.025 1

0.005

1

0.8

0 1

0.8 0.6

0.8

1 0.6

ξ2

0.6

0.8

0.4

0.4

0.6

0.4 0.4

0.2

0.2 0

0

ξ2

ξ1

0.2 0

Figure 14: Shape functions for twist vectors

Edge number

0.2

ξ1

0

,

Table 4: The shape functions for edge

! 1 !21 !/ 1 ! / ! 1 7 , / ! 1 ! 1 ! / ! 1

7 ,/ ! /

7 ,/ , %

system of equations can be expressed as,

,

(3.48)

where and are functions in terms of four unknown twist vectors and four unknown vector edge coefficients , respectively. The parametric function for vertex contribution in equation

24

0.01 0 0.008

ψe21

ψe11

−0.002 −0.004

0.006

−0.006 0.004 −0.008 −0.01 1

0.002

0.8

0 1

1 0.6

0.8

0.8

1 0.6

0.8 0.6

0.4

ξ

0.6

0.4

0.4

0.2

2

0

ξ2

ξ

0.2

1

0

0.4

0.2

ξ

0.2 0

1

0

0.01

0

0.008

ψe31

ψe41 0.006

−0.005

0.004 −0.01 1

0.002 0.9 0.8 1

0.7 0.6

0 1

0.8 0.5

0.8 0.3

ξ2

1

0.6

0.4

0.6

0.8

0.4 0.2 0

ξ2

ξ1

0.2

0.1

0.6

0.4 0.4

0.2

0.2 0

0

Figure 15: Edge shape functions (3.45) can be reduced to a matrix form,

0

/ ! ! 1 $ / ! 1 / ! 1 / ! 1 / ! 1 ( $ / ! 1 / ! 1 / ! 1 / ! 1 (

We have,

/

7 , )/ !! 11 11 7 , /! 11 , //

, /! 11 / 7

7 , /

/

7 , )/ !! 11 11 / 7 , / ! 1 1 / 7 , / ! 1 1 / 7 , /

ξ1

(3.49)

$ / ! 1 / ! 1 ( +

(3.50) 25

and,

/ , / !2!211 11 //

,, '5// 2! 1 1 / ! 1 / , '/ !21 1

(3.51)

is a matrix prescribed in terms of tangent vectors

/

7 , / !! 11 11 7 , '/ ! 1 1 //

7 , /

/ 7 , / ! 1 1

/! 1

/ !21

,

/

7 , )/ !2!211 11 / 7 , / !21 1 / 7 , / !21 1 / 7 , /

$ / ! 1 / ! 1 (

(3.52)

/ !216"

As is a seventh order polynomial , vanishing at the endpoints of the edge, enforc ing condition (2.5) is equivalent to enforcing,

(3.53) / 1 % & & 1 equations, we get values of the components of with . Solving a system of / & eight vector coefficients and + & '' .

3.2.2 Middle node Contribution Mathematical formulation of any boundary value problem consists of a differential equation and boundary conditions. The connection between transfinite interpolation and the boundary values problems is explicit in the last section. In this case, we exactly interpolate the prescribed boundary conditions. However, the behavior of the interpolant away from the boundaries is quite arbitrary. Thus, the solution of a boundary value problem can be viewed as the construction of a function that extends the boundary conditions into the domain, with differential equations playing the role of a constraining or smoothing operator. The interpolation problem has no unique solution; there are infinitely many functions interpolating any given data. An unique construction can be established by putting additional constraints on the interpolant. Often such constraints appear as the minimization of some quantity. Many interpolation schemes use minimization of energy [27] as a means for controlling the shape of the interpolant,

/ / ! ! 1 1 *

26

(3.54)

Figure 16: A hexahedron with two adjacent faces on a spherical surface

where

,

is the rectangular reference domain, and

/ ! ! 1 '/ ! ! 1 5/ ! ! 1 (3.55) ! ! 1 interpolating the boundary conditions given in equation (3.42); with the first term / unknown for the middle node; are shape function described be : / 1 is the vector fore. The value of the coefficient is obtained by minimizing (3.54). Upon differentiating (3.55) with respect to & '' , we construct and solve a system of three linear equations for components of , (3.56) / 1

Now we conclude the discussion with some simple examples. Using the existing curve net for a spherical patch in Fig.6, and applying the surface fitting technique above, we can get a hexahedron with two adjacent faces on a spherical surface, illustrated in Fig.16. In the similar way, we can obtain a hexahedron with two adjacent faces on a cylindrical surface, see Fig.17.

3.3 Modeling of a Human Head The discussed continuous geometry reconstruction technique has been applied to model the geometry of a human head, necessary for high accuracy 3D hp FE simulation. The implementation 27

Figure 17: A hexahedron with two adjacent faces on a cylindrical surface has been done within our Geometrical Modeling Package (GMP) [23] interfaced with software LBIE - Mesh Level Set Boundary and Interior-Exterior Mesher, developed at Center for Computational Visualization (CCV) at ICES. Given an MRI scan of a human head, we proceed in the following steps [28, 29], Pre-Processing. We first use the anisotropic diffusion method coupled with bilateral prefiltering to remove noise from imaging data. Depending on the application, suitable isosurfaces are selected for the interval volume by using the contour spectrum and the contour tree. Hexahedral Meshing. We then begin to extract 3D meshes from the interval volume, and a feature sensitive error function is adopted to reduce the number of elements while preserving features. Quality Improvement. Finally, the edge contraction method is then used to improve the mesh quality. The resulting trilinear hex mesh provides a minimum coarse representation of the head topology and a geometrical information on the surface of the head. The geometry information is given by labeling the surface nodes and specifying the corresponding normals. The topological (hex to points connectivities), and geometrical (normals)data is then imported into GMP, where the actual continuous geometry reconstruction takes place.

28

η3

5

8

●

8

●

ξs 2 ξs 1

5

(2)

6

7

6

7

9 12 (6) ξs 2

(5)

ξs 2

ξs 2

(3)

11 10

ξs 1

(4) ●

ξs 1

1

ξs 1

●

ξs 1

ξs 2 ξs 1

ξs 2

2●

4

1 (1)

4

η2

3

●

2

3

η1

Figure 18: A hexahedron in the reference frame The recently rewritten and upgraded GMP [30] supports the construction of exact parameteri zations for a general class of 2D (BEM) and 3D (FEM) manifolds in [31, 30]. We have explored a number of novel geometric modeling techniques and implemented them in the GMP. The GMP not only supports the construction of exact parameterizations for a general class of 2D [32] and 3D [33] objects, but also provides the derivatives of the mappings with respect to reference coordinates for any given points in reference frame. In our geometric modeling, a 2D object is represented as a union of curvilinear triangles or rectangles, while a 3D object is represented with a FE-like mesh of curvilinear hexahedral blocks. The geometry of the object is described then by constructing parameterizations for each of the blocks. The GMP model later can be used to generate the actual FE meshes of arbitrary high order, and to support geometry updates during mesh refinements. Each of the local edges or faces in the GMP has its own global orientation. Adjusting edge and face parameterizations for orientations involves transforming local edge and face coordinates into the global coordinates. We must ensure the compatibility of parameterizations. For example, parameterizations for adjacent hexahedra must be compatible with each other. Roughly speaking, when a face is shared by two hexahedra, and two parameterizations are used, the same FE on both sides of the face mesh must be obtained.

29

The coordinate transformations are handled in a hierarchical manner. First the coordinates of points are set. Next, the curve parameterizations are constructed in such a way that they conform to the existing vertex coordinates. Next, the parameterizations for rectangles must conform to the parameterization of the curves. Finally, in the most complicated case, the parameterizations for hexahedra must conform to the parameterizations of the rectangles, see Fig.19. The construction of the human head model is based on two GMP parameterizations: Transfinite Interpolation Rectangles and Transfinite Interpolation Hexahedra, both based on the classical transfinite interpolation and linear blending functions technique [34, 35]. The interface between ζ2

ζ1

x3 x2 x1

xH

xR η3

ζ2

η2

η1 η3

ζ1

ξ2 η2

ξ (ζ)

ηf

ξ1

η1

Figure 19: Compatibility of parameterizations for a hexahedron GMP and geometric data extracted from an MRI scan [36] provides connectivities of for a trilinear hex-mesh and topology of the boundary surface. Using this information, we then reconstruct a curvilinear hex-mesh with a continuous representation of the surface. The obtained hex-mesh head model is then used to generate the actual meshes for hp-Adaptive FE simulations. We proceed in the following steps: Step1: Linear Hex-Mesh Model Reconstruction Before generating a full nonlinear model, it is convenient to check the topology of the GMP model by reproducing the piecewise trilinear mesh within the GMP using the simplest trilinear parametrization for the GMP hexahedra. Note that all twelve edges of each hexahedron are straight line segments connecting two adjacent nodes. The geometric data required by constructing the linear model includes the indexes and coordinates of the eight nodes for each hexahedron. Contrary to the minimal data provided in the CCV model, the GMP stores all topological connectivities in a hierarchical 30

manner: points to edges and edges to points; edges to faces and faces to edges; faces to hexas and hexas to faces. This extended connectivity information is generated through a global sequential search, processing one hexahedron at a time. First the local entities are connected to the corresponding global entities, i.e., vertices to points, edges to curves, and faces to rectangles, and then the overall global entities are updated[23]. The topological connectivity information is sufficient to generate the simplest geometry model based on linear curves, bilinear faces, and trilinear hexahedron parameterizations. Step2: Obtaining Topological Information on the Boundary Surface. The necessary topological information for the construction of curvilinear hexahedra is collected. The geometric data provided by CCV includes labeling the nodes on the head surface and providing the outward normals at the nodes. We loop through the GMP rectangles, and specify those with all four vertices on the boundary surface. The rectangles are redefined as a new rectangle type : biquartic rectangle (biqua), a rectangular patch on continuous surface. Rectangles with at least one vertex on the boundary surface are redefined as: transfinite interpolation rectangles (TraQua). The remaining rectangles are still stored as bilinear quadrilaterals( BilQua). In the same way, we designate those curves with both vertices on the surface as a new curve type: Points with Coordinates and Normals (CoorNrm). The remaining curves are marked with the label of a segment of a straight line (Seglin). All points are classified as regular points (Regular) and all hexahedra are of GMP hex type: transfinite interpolation hexahedra (TriLiHex). Step3: Curvilinear Hex-Mesh Model Reconstruction. The surface reconstruction scheme discussed in the last few sections is performed. A cubic curve-net with compatible cubic normals along the curves CoorNrm is first constructed, and then the surface fitting techniques is used to generate the biqua rectangles. Before generating the whole human head model, we first test the scheme on some feature parts of the head, e.g., the nose. Any hexahedron that constitutes a part of the nose is a special case because all its eight vertices are on the surface. Following the three steps for the linear nose model, see the first figure in Fig. 20(a), we obtain a curvilinear model, illustrated in in Fig. 20(b). Fig.21 displays a zoom on the bottom hexahedron on the nose for which the curvilinear geometry is the lest visible. The entire reconstructed human head model is presented in Fig.22, with a statistics of the model summarized in able 5. The motivation for this work is to eventually simulate EM waves in a human head. This involves enclosing the head within a truncating sphere, and meshing the entire volume within the 31

Figure 20: The nose model of the human head. NRSURFS NRPOINT 0 25744

NRCURVE 8448

NRTRIAN 0

NRRECTA 4424

NRPRISM 0

NRHEXAS 704

Table 5: Numbers of global entities in GMP for the human head model sphere, and the head. Special absorbing boundary conditions are imposed on the truncating sphere to model the interaction with the rest of the space. We will use the model later to generate the actual meshes for -Adaptive Finite Element (FE) simulation of Electromagnetic (EM) wave.

4 Other

&

Interpolation Methods

In this section, for a comparison, we report on two other (less successful) niques implemented earlier. 32

reconstruction tech-

Figure 21: The bottom hexahedron of the nose model

4.1 Inverse Distance Weighting (IDW) Interpolation In this scheme, we use one of the most common techniques for surface interpolation, Inverse Distance Weighting Interpolation (IDW) [37]. IDW assumes that the value of the interpolant should be influenced more by the value of the nearby points and less by more distant points. The main advantage of the inverse distance weighting method is that it does not place any restriction on the topology of the sets being interpolated. The starting point for IDW interpolation is a set of , representing a distance from specific so called normalized distance functions geometrical object e.g., points or lines. In our case, the interpolation takes place over the reference rectangle, and four normalized distance functions will simply represent the distance from the rectangle’s edges.

/ 1 &

1 & , we define the corre / 1 / + & 1 as: 1 - / (4.57) 1 - /

Given any set of the normalized distance functions sponding inverse distance weighting functions

/

/ 1

1

/

is the normalized distance to the point. The exponents where interpolating functions as follows: 1. 2.

&

&

: the interpolant is not differentiable points

: the interpolant is differentiable

33

&

where

/

times at the points.

control behavior of the

1

Figure 22: The head model as a union of curvilinear hexahedra

/ 1

The inverse distance weighting functions form a partition of unity in the sense that . This partition of unity also assures completeness of the system and guarantees reproduction of constant functions.

&

, / !21

7, / ! 1

We proceed now with the discussion of the actual interpolation technique. As in the previous and normals defined along the sections, we start with a reconstruction of curves curves. The curve reconstruction is once again based on Hermite representation (2.16). To obtain the two tangent vectors at both endpoints, we use the geometric information on the adjacent endpoints. We select for the tangent vector to be a vector perpendicular to lying in the plane spanned by and vector . A simple geometrical manipulation yields:

where data .

/ / / /

11

11

(4.58)

is directed toward point . Note that is directed toward point and that in case of , reduces to . In the same way, we obtain the tangent vector at point

34

Figure 23: The cross boundary derivatives (CBD’s).

!1 7 , / 7 , / 2! 1 8 7 , / ! 1 !1 !21 ! 2 1 , / , / (4.59) 7 , / , / 2! 1 , / ! 1 , / ! 1 , / ! 1 !21 and projecting it onto the plane perpendicular to the tangent vector , / , !, 21 7 !21 7 / , / (4.60) 7 , / ! 1 !21 7 , / !21 7 , / 2 ! 1 We emphasize that normals 78, / are not polynomials, in fact, they are not even rational functions.

As the IDW interpolation produces non-polynomial parametrizations anyway, we are not bounded to a polynomial variation of the normals , and we can use geometrically simpler techniques. We construct the normal vector along the curve by taking a linear combination of normal vectors and ,

Once we have selected the normals, we proceed now with the patch interpolation. The construction is based on the selection of appropriate Cross Boundary Derivative (CBD) vectors. For each of the four edge points, we choose the CBD vectors to be the cross products of the corresponding tangent and normal vectors at the point, as shown in Fig.23. For a point on the first curve, we have,

/! 1

, /! 1

7 , / ! 1 , / ! 1

(4.61)

In a similar way, we can get the CBDs for the other three curves as follows,

, '/ !! 11

7 , '/ !! 11 , , /! 1

7 , / ! 1 , , /

7 , / ,

35

/ !! 11 /! 1 /

(4.62)

, - : & ''

We have four edge functions , along with the corresponding CBD vectors . We use a IWD Interpolation to derive a parameterization for the patches. The . The four parametrization is constructed over the reference domain, with

,.- & ''

9 /! ! 1

! ! " $ % '&)(

Figure 24: Normalized distance functions in reference frame normalized distance functions represent the distance from the four edges, which in turn are in this case can be described on boundaries shown in Fig.24. Each portion of the boundary is represented implicitly by the functions as follows:

/ !! !! 11 ! 6 !

/ & 6

/ !! !! 11 8 !& !

/ In order to guarantee that the interpolants have first derivatives, we take in equation (4.57). 1 The corresponding weighting functions / & '' , can now be expressed in terms of 9, ! /0& ! 1 /0& ! 1 /! ! 1 ! / & ! 1 /0& ! 1 ! ! / & ! 1 ! ! /0& ! 1 ! /0& ! 1 / & ! 1 ! ! 1 The interpolating function / for the rectangular patch is constructed as a linear com bination of edge functions , - & 5' with the corresponding weight functions & 5' , illustrated in Fig. 25, ! ! 1 (4.63) / ,.-0/ ! ! 1 / ! ! 1

The selection of edge functions

,.- &

36

ensures

continuity and it is done in the

following way,

, / !! !! 11 , '/ ! ! 1 , / ! ! 1 , /

,.- : & ! ''

/ !! 11' ! , '/ ! 1 / & , / ! 1 '! / & , / ,

, / ! ! 1

, ! ! ,

, are the corresponding

, / ! 1 , / ! 1'! , / ! 1 , / ! 1 , / ! 1

(4.65)

, / ! 1 '! ! 1 !/ ! 1 ,! / ! 1 , / ! ! !/ ! 1 ! 1 ! 1 , / , /

Similarly for derivatives,

(4.64)

,.- & ''

where , are the Hermite curves, and CBDs. With point approaching the first edge, we have,

! 1 ! 1

, / !! 11 , 5/ ! 1 , / ! 1 , '/

(4.66)

The limitation behavior of the edge functions ensures continuity of the patch parametrization across edges shared by two patches. Consider the value of the interpolant of the rectangular patch as point approaches the first edge,

/! ! 1

9

! ! 1 / , / ! ! 1 , & , , / ! ! 1 ,

/! ! 1

/! !

/! 1

1

,

,

(4.67)

Consequently, two adjacent rectangular patches have the same function values at the common edge, the interpolant approaches the physical coordinates of the boundary edge when approaches the corresponding reference coordinates of that edge. The derivatives of with

9

37

Figure 25: Physical domain defined on a reference domain respect to

!

!

and

behave in a similar way,

/! ! 1 !

/! ! 1 !

,.-0/ ! ! 1 !

/! ! 1 ! ! , & ! , ! , ! , / ! ! 1 ,! ! , / ! ! 1 , / ! 1

! ,

(4.68)

Note that IDW interpolation does not place any restriction on the topology of the curves being interpolated. We exam this scheme with the special case depicted in Fig.21. Interpolating two rectangular patches over the implicitly defined curves, we get a curvilinear hexahedron with its three adjacent faces on a surface, see Fig.26. The IDW interpolation ensures continuity along the edges but, in general, guarantees only a -continuity at vertices. The main disadvantage of using the inverse weighting functions is their loss of regularity at vertices. At each vertex of the reference domain, the inverse weighting function is discontinuous and its derivative blow up to infinity, see Fig.27. It is not clear how this lack of regularity will affect the mesh generation and quality of meshes. 38

Figure 26: A

4.2 Bicubic

continuous hexahedron

Surface

/! ! 1 " 78, / ! 1

In this section, we report shortly on an earlier attempt of rectangular patch reconstruction using . With bicubic polynomial reconstruction, we have only bicubic polynomials, three free parameters per edge (the twist vectors defined at vertices). Consequently, enforcing to the curve tangents is possible only if the degree of and orthogonality of normals do not exceed four. In the presented attempt, we assume that we start with a simple linear interpolation of the normal,

7 , / !21

, / 2! 1

7 , / !21 /0& !21 ! (4.69) !21 and only then follow with a reconstruction of the curve , / done in such a way that the orthog!21 7 , / ! 1 , is satisfied. The order in which we reconstruct curves !21 , / onality condition / , / ! 1 and normals 7, / ! 1 is then opposite to the biquartic reconstruction scheme, first the normals

and then the curves.

Using the same Hermite curve representation (2.16), we have to determine two unknown vector coefficients and . In order to ensure compatibility condition (2.2), we need

39

(4.70)

The first weight function W1(ξ1,ξ2)

The second weight function W2(ξ1,ξ2)

1

1

0.8

0.8

0.6

0.6

W2

W1 0.4

0.4

0.2

0.2

0 1

0 1 0.8

0.8

1

ξ 0.6 2

0.4

0.4 ξ

0.2

2

0.8 0.4

1

0.2 0

1

ξ 0.6

0.8 0.6

0.4 ξ

0.2

The third weight function W3(ξ1,ξ2)

1

0.2 0

0

0.6

0

The fourth weight function W4(ξ1,ξ2)

1

1

0.8

0.8

0.6

0.6

W3

W4

0.4

0.4

0.2

0.2

0 1

0 1 0.8

0.8

1

ξ20.6

0.4 ξ

0.2

0.2 0

1

ξ20.6

0.8 0.4

0.6

1

0.8 0.4

0.4 ξ

0.2

0.2 0

0

0.6

1

0

Figure 27: Weight functions of IDW The scalar function

/! 1

/! 1

Note that to enforcing,

/ !21

can be written as a function of unknown vectors

and

as:

/ 11 / 11 ! $ / ! 1 : '/ !21 / ! 1 / !21 ( / 1 / 1 & ! (4.71) / 1 / 1 / / " as ./ !21 " and 7 / !21 " . Enforcing / ! 1 (2.3) is thus equivalent

"

/ & 1 / 1

(4.72)

Unknown vector coefficients can be written uniquely as a combination of the three linear independent basis vectors. The concept of the dual basis enables the decomposition of a 40

vector in a nonorthogonal basis. Looking at equation (4.70), we construct three basis vectors as:

(4.73)

and are not necessarily orthogonal, can not be given by the inner

Because basis vector product of and . In

order to calculate the corresponding components, we introduce another set of basis vectors , called the dual of , defined by the orthogonality relation where delta. In terms of the dual basis , . The components of the is the Kronecker

4/ 1 1 4/

vector can then be calculated as

. This implies that,

4/ 4/

11

(4.74)

The two unknown vector coefficients can now be written in terms of the dual basis,

(4.75)

(4.76) / ! 1 / !21 / !21 / ! 1 % with, / ! 1 / 11 ! / !2! 11 !21 ! 1 (4.77) / / 1 /0& ! 1 / !21 1 ! !21 ! 1 / / /0& / !21 / 1 ! / !21 1 ! 1 / /0& : '/ / / (4.78) 1 1 Requesting / and / , we get a linear system of equations, / 11 / 11 / 11 % (4.79) / / / and . Finally, the remaining degrees of freedom of the curve and solved for components the energy function, are determined by minimizing / 1 , ! !21 / ! 1 / 1 / ! 1 / 1 !21 !

/ / * (4.80) Consequently,

41

Taking derivatives of

gives: 1 ! 1 1 ! 1 1 / !21 / !21 / / / / / 1 / !21 / !21 / 1 / ! 1 / ! 1 !

1 !21 !21 / 1 / ! 1 1 / ! 1 / 1 / !21 / !21 / / 1 / ! / 1 / ! / 1 !

in terms of the unknowns

Solving the resulting system of linear equations yields,

/

1 / 1

/

1

(4.81)

(4.82)

curve under the condition that vectors The prescribed procedure results in a uniquely defined and , used to defined the original basis vectors , are linearly independent.

The bicubic rectangular patch can be written as a Hermite patch in equation (3.50), with the four CBDs computed in the standard way,

/ ! 1 / ! 1 / ! 1 / ! 1 (4.83) / ! 1 / ! 1 /! 1 /! 1 Note that vector CBDs are also third order polynomials . In order to interpolate a surface, !21

/ ! 1 / ! 1 / & ! 1 / ! 1 / ! 5& 1 / ! 1 / % ! 1 / ! 1

: / !! 11 : / ! 1 : / ! 1 /

the rectangular patch must satisfy the system of four equations in (2.5). Introducing product of CBDs and its corresponding normals,

/ !21

/

7 , / 7 , / 7 , / 7 , / 7 , / 7 , / 7 , / 7 ,

/

7 , / !! 11 11 /

7 , '/ ! 1 1 / 7 , / ! 1 1 / 7 , / !21 1 / 7 , / !21 1 / 7 , / !21 1 / 7 , / !21 1 / 7 , /

/ !2!211 11 / !21 1 $ / !21 1 / !21 1 / !21 1 / !21 1 $ / !21 1 /

/ ! 1 / ! 1 (

/

! 1 /! 1 (

"

/0& 1 / 1 / 1

is a dot

(4.84)

We reduce the compatibility condition for the rectangular patch to equation is a fourth order polynomial, , it is sufficient to enforcing,

/

42

/ ! 1 . Since / ! 1

(4.85)

8& ''

Solving the resulting system of three equations for four edges (twelve equations in total), we obtain . the components of the four twist vector The bicubic scheme produces results similar to the examples of the spherical and cylindrical rectangular patches shown in Fig.16 and Fig.17. The major drawback of the presented construction lies in a possibility of a degenerated case in which the procedure fails. This happens when normal vectors and are parallel to each other and, at the same time, are not orthogonal to . In the case of , the minimization principal can be used to yield a straight line segment. It was the degenerated case that demonstrated an insufficiency of linear variation of and prompted us to use the biquartic parameterizations.

7 , / !21

5 Conclusions and Future Works The paper presents results of a preliminary study on geometric reconstruction in context of geometries reproduced from MRI scans and mesh generation for high order FE discretization. The presented biquartic scheme seems to be the lowest order continuity construction for general unstructured meshes. The polynomial parametrizations are inexpensive to compute and guarantee computations. In the case of the presented high regularity of parametrization necessary in IWD technique, it may not not satisfied. It is not clear at this point, however, how the ( and not ) regular parametrization will affect the convergence rates of high order methods. The important property of the presented reconstruction scheme is that it remain uniformly stable in the case of degenerated geometrical data (see the Appendix). Among other tasks, we intend also to collaborate with CCV on multi-resolution techniques and hierarchical geometry reconstruction schemes. At this point, the information on geometry contained in the original fine mesh reconstruction, during the coarsing stage is reduced to normals only. Ideally, the geometry reconstruction on the coarse grid should conform to the fine grid representation in a more elaborate, multi-resolution model. We intend to address these topics in our future work.

43

A Singular Value Decomposition (SVD) n1 = (0, 0, 1)

p1 = (0, 0, 0)

p2 = (0, 1, 0)

0.04

n2 = (1, 0.01, 0) n2 = (1, 0.05, 0) n2 = (1, 0.1, 0) n2 = (1, 0. 05, 0) n2 = (1, 1, 0) 0

0.02

z 0

0.02

−0.02 0.04 −0.04 0.06

0

x

0.2 0.4

y

0.08

0.6 0.8 0.1

1

n1 = (0, 0, 1)

p1 = (0, 0, 0)

p2 = (0, 1, 0)

0.04 0.03 0.02 0.01

z

0

−0.01

n2 = n (1, =1,(1, 0) 0. 05, 0) 2 n = (1, 0.1, 0) 2 n = (1, 0.05, 0) 2 n2 = (1, 0.01, 0)

−0.02 −0.03 −0.04 0 0.2 0.4

y 0.6 0.8 1 0.1

0.08

0.09

Figure 28:

0.07

0.06

0.05

x

0.04

0.03

0.02

0.01

0

varies from (1, 0.01, 0) to (1, 1, 0)

A very powerful set of techniques dealing with sets of equations or matrices that are either 44

singular or numerically very close to singular is the so-called singular value decomposition (SVD) [38, 39, 40, 41]. SVD allows one to diagnose problems related to solution of linear systems provides numerical a answer as well.

matrix can be written as the product of an column-orthogonal Any matrix , an diagonal matrix with positive or zero elements, and the transpose of an orthogonal matrix :

where

The diagonal elements of matrix

(A.86)

(A.87)

are the singular values of matrix

;

& (A.88) is a square matrix, it is also row-orthogonal, & . Equation (3.40) defines Since " & , as a linear mapping from the vector space into itself. If is singular, , called the null space / 1 , mapped to zero, / / 1 1 . then there is some subspace of 1 is called The dimension of the null space is call the nullity of A. The dimension of the range / 1 / 1 . the rank of . The relevant theorem is “rank plus nullity equals N”, / 1 and / 1 . Specifically, the columns of SVD explicitly constructs orthonormal bases for / set of basis vectors that whose same-numbered elements are nonzero, are an orthonormal span the range; the columns of whose same-numbered elements are zero, are an orthonormal basis for the null space. " , one has to: To find the singular value decomposition [42] of the matrix

1. Find the eigenvalues of the matrix

and arrange them in descending order.

2. Find the number of nonzero eigenvalues of the matrix

.

3. Find the orthogonal eigenvectors of the matrix corresponding to the obtained eigenvalues, and arrange them in the same order to form the column-vectors of the matrix .

"

45

"

4. Form a diagonal matrix of first eigenvalues of the matrix

placing on the leading diagonal of it the square roots got in step 1 in the descending order.

5. Find the first column-vectors of the matrix 6. Add to the matrix process.

the rest of

" ,

& ''

vectors using the Gram-Schmidt orthogonalization

&'

In our case, matrices and are square , so they are orthogonal and their inverse matrices are equal to their transposes. The inverse matrix of becomes,

4/ 1 / 1

(A.89)

If one or more are zero or very close to zero, matrix is singular or numerically singular in this case. Equation (3.40) can then be used to obtain the solution:

/ 1

Where the diagonal elements of matrix

(A.90)

are given by

&

(A.91)

which being a singularity threshold.

/ 1

In the case when is singular, we need to check whether the vector on the right hand side lies in . If it does, the singular set of equations has multiple solutions. If it does not, the problem (3.40) has no solution. When ’s are very small but nonzero, the matrix is ill-conditioned. In that case, the direct solution methods of decomposition or Gaussian elimination may actually give a formal solution to the set of equations, i.e., a zero pivot may not be encountered, but the solution vector may have wildly large components whose algebraic cancellation, when multiplying by the matrix , may give a very poor approximation to the right-hand vector . In such cases, with 0 will the solution vector obtained by zeroing the small ’s is often better. Replacing provide the closest that minimizes the distance to in the least square sense,

*

&

(A.92)

whith denoting the residual of the solution.

, / ! 1 "

The simplest example is of an ill-conditioned situation deals with the parametrization of a cubic curve degenerating into a straight line segment . The geometric data for the 46

−3

x 10

n1 = (0, 0, 1)

1

p1 = (0, 0, 0) p2 = (0, 1, 0)

0.8 0.6 0.4 0.2

n = (1, 0, 0)

Z0

2

n = (1, 0.0001, 0) 2

−0.2

n = (1, 0.001, 0) 2 n = (1, 0.0005, 0)

−0.4

2

n = (1, 0.01, 0) 2

−0.6

n2 = (1, 0.005, 0)

−0.8 −1 0 0.2 0.4

y 0.6

0 0.5 0.8

1

x

−3

x 10

1.5 1

2

n1 = (0, 0, 1)

−3

x 10

p1 = (0, 0, 0) p2 = (0, 1, 0)

n2 = (1, 0.0001, 0) n2 = (1, 0, 0)

1 0.8

n2 =0) (1, 0.001, 0) n2 = (1, 0.0005,

Z0.6 0.4

n2 = (1, 0.005, 0) n2 = (1, 0.01, 0)

0.2 0

0

−0.2

0.5

−0.4

x

−0.6

1

−0.8 1.5 −1 0

0.2

Figure 29:

y0.4

0.6

0.8

2 1

varies from (1, 0.01, 0) to (1, 0, 0)

47

−3

x 10

/%%'& 1 6 /0& 1

/%% 1 /%5& 1

limity case are: and . The vanishing constraints in (3.36) result in a singular matrix . Using SVD, we study the behavior of the matrix as . We have use the curve reconstruction routine, with data . The code delivers uniformly stable results converging to the limity case. For illustration, Figures A A show the results of the curve reconstruction for values varies from to , and from to . Note we use the different scales in Figures to illustrated the convergence property.

1 * / % % 5 & ' & 1 '& '& / '& '& 1 /% %& '& 1

/ %& '& 1 /% '& 1

48

/ %'&

B Subroutines The following subroutines have been implemented into the GMP package. input head.f

the interface with data extracted from MRI that stores all the connectivity informations.

filter.f

A subroutine used to remove extra points that are redundant in the input file.

HermCur1.f

A subroutine which generates rational curve function with CBDs

ratio.f

A subroutine that reconstructs the

HermCur2.f

A subroutine which generates Hermite curve function with linear normals

HermRec.f

A subroutine that reconstructs the surface by using bicubic Hermite rectangular patch interpolations.

HermCur3.f

A subroutine which generates Hermite curve function with cubic normals

biqua.f

A subroutine that reconstructs the surface by using biquartic Hermite rectangular patch interpolations.

surface by using IDW interpolation.

matr.f

A subroutine that solves a system of linear vector equations using Guassian elimination.

svdcmp.f

A subroutine that solves a system of linear vector equations using SVD.

head

The txt file that is used as a input for the interface

49

References [1] G.Xu, C.Bajaj. “Adaptive Model Reconstruction by Triangular A-patches.” TICAM Report, 2003 [2] Bajaj C., Chen J., Xu G. “Modeling with Cubic A-patches.” ACM Transactions on Graphics, vol. 14, 103–133, April 1995 [3] Bajaj C., Chen J., Xu G. “Modeling with Cubic A-patches ACM Transactions on Graphic.” 14, vol. 2, 103–133, April,(1995) [4] Xu G., Bajaj C., Huang H. “C1 Modeling with A-patches from Rational Trivariate Functions Computer Aided Geometric Design.” 18, vol. 3, 221–243, 2001 [5] B.Guo. “Modeling Arbitrary Smooth Objects with Algebraic Surfaces.” PhD thesis, vol. Computer Sciences Cornell University, 1991 [6] T.W.Sederberg. “Piecewise Algebraic Surface Patches.” Computer Aided Geometric Design, vol. 2, 53–59, 1985 [7] H.P.Seidel. “Symmetric Recursive Algorithms for Surfaces: B-patches and the De Boor Algorithm for Polynomials Over Triangles.” Constr. Approx., vol. 7, 257–279, 1991 [8] C.T.Loop, T.DeRose. “Generalized B-spline Surfaces of Arbitray Topology.” Computer Graphics, vol. 24, 347–356, 1990 [9] W.Boehm, G.Farin J. “A Survey of Curve and Surface Method in CAGD.” Computer Aided Geometrc Design, vol. 1, 1–60, 1984 [10] G.Farin. “Curves and Surfaces for Computer Aided Geometric Design.” Academic Press, 1990 [11] Peters J. “C1 Surface Splines.” SIAM J. of Numer. Anal., vol. 32-2, 645–666, 1995 [12] Peters J. “Biquartic C1-surface Splines Over Irregular Meshes.” CAD, Jan 95 [13] Peters J. “Surfaces of Arbitrary Topology Constructed from Biquadratics and Bicubics.” Designing fair curves and surfaces, vol. SIAM Publications, N. Sapidis, 1994 [14] Peters J. “Smooth Free-form Surfaces Over Irregular Meshes Generalizing Quadratic Splines.” CAGD, pp. 347–361, 1993

50

[15] A.R.M.Piah. “Construction of Smooth Surfaces by Piecewise Tensor Product Polynomials.” CS Report, vol. 91/04, University of Dundee, UK, 1991 [16] O.Reif. “Biquadratic G-spline Surfaces.” Preprint 93-4, vol. Math Ins A, University Stuttgart, 1993 [17] M.Sabin. “Non-rectangular Surface Patches Suitable for Inclusion in a B-spline Surface.” Proceedings of Eurographics, vol. North Holland, 57–69, 1983 [18] T.N.T.Goodman. “Closed Surfaces Defined From Biquadratic Splines.” Constructive Approximation, vol. 7, 149–160, 1991 [19] E.Catmull, J.Clark. “Resursively Generated B-spline Surfaces on Arbitary Topological Meshes.” CAD, vol. 10, 350–355, 1978 [20] D.Doo. “A Subdivision Algorithm for Smoothing Down Irregularly Shaped Polyhedrons.” Proceedings on interactive techniques in computer aied design., vol. Bologna, 157–165, 1978 [21] L.Shirman, C.Sequin. “Local Surface Interploation with Bezier Patches.” Computer Aided Geometric Design, vol. 4, 279–296, 1990 [22] R.F.Sarraga. “ Interplation of Generally Unrestricted Cubic Bezier Curves.” Computer Aided Geometric Design, vol. 1-2, 23–40, 1987 [23] Xue D., Demkowicz L. “An Interface Between Geometrical Modeling Package(GMP) and Mesh-Based Geometry(MBG).” ICES Report, vol. 03, 2003 [24] M.A.Sabin. “Transfinite Surface Interpolation.” Proceedings of a sixth conference on Mathematics of Surfaces, 1996 [25] Peter J. “Smooth Interpolation of a Mesh of Curves.” Constructive Approximation, vol. 7, pp.221–246, 1991 [26] Mortenson M.E. Geometric Modeling. Wiley Computer Publishing, 1997 [27] Hoschek J., Lasser D. “Fundamentals of Computer Aided Geometric Design.” A K Peters, 1993 [28] Zhang Y., Bajaj C., Sohn B.S. “Adaptive and Quality 3D Meshing from Imaging Data.” Proceedings of 8th ACM Symposium on Solid Modeling and Applications, vol. June 16-20, 286–291, 2003

51

[29] Zhang Y., Bajaj C., Sohn B.S. “3D Finite Element Meshing from Image Data.” The special issure of Computational Methods in Applied Medhanics and Engineering on unstructral mesh generation, 2003 [30] Xue D., Demkowicz L. “Geometrical Modeling Package Version 2.0.” TICAM Report, vol. 02-30, 2002 [31] Demkowicz L., Bajer A., Banas K. “Geometrical Modeling Package.” TICAM Report, vol. 02-06, 1992 [32] Demkowicz L. “2D hp-Adaptive Finite Element Package(2Dhp90) Version 2.0.” TICAM Report, vol. 02-06, 2002 [33] Demkowicz L., Pardo D., Rachowicz W. “3D hp-Adaptive Finite Element Package(2Dhp90).” TICAM Report, vol. 02-24, 2002 [34] Gordon W.J., Hall C.A. “Transfinite Element Methods: Blending Function Interpolation over Arbitary Curved Element Domains.” Numer. Math., vol. 21, 109–129, 1973 [35] Gordon W.J. “Blending Function Methods of Bivariate and Multivariate Interpolation and Approximation.” SIAM J. Numer. Anal, vol. 8, 158–177, 1971 [36] Zhang Y., Bajaj C., Sohn B.S. “Adaptive Multiresolution and Quality 3D Meshing from Imaging Data.” TICAM Report, vol. 02-42, 2002 [37] Shepard D. “A Two-dimensional Interpolation Function for Irregularly spaced Data.” 23rd ACM National Conference, pp. 517–524, 1968 [38] Gentle J.E. “Singular Value Factorization.” Numerical Linear Algebra for Applications in Statistics, vol. Berlin: Springer-Verlag, 102–103, 1998 [39] Golub G.H., Van Loan C.F. “The Singular Value Decomposition.” 2.5.3 and 2.5.6 in Matrix Computations, vol. 3rd ed. Baltimore,Johns Hopkins University Press, 70–71, 1996 [40] Nash J.C. “The Singular-Value Decomposition and Its Use to Solve Least-Squares Problems.” Compact Numerical Methods for Computers: Linear Algebra and Function Minimisation, vol. Bristol, England: Adam Hilger, 30–48, 1990 [41] Press W.H., Flannery B.P., Teukolsky S. A.and Vetterling W.T. “Singular Value Decomposition.” Numerical Recipes in FORTRAN: The Art of Scientific Computing, vol. Cambridge University Press, 51–63, 1992

52

[42] Press.W.M, W.M.Flannery, B.P.Teukolsky. Numerical Recipes: The Art of Scientific Computing. New York, NY, Cambridge University Press, 1986

53

SURFACES RECONSTRUCTING OF WITH BIQUARTIC PATCHES AND OTHER TECHNIQUES Dong Xue, Leszek Demkowicz Texas Institute for Computational and Engineering Science The University of Texas at Austin Austin, TX 78712

Waldek Rachowicz Cracow University of Technology, Poland Abstract

We present in this paper an approach to reconstruct a surface that interpolates given surface is represented in a parametric form vertices with prescribed normal vectors. The implemented in two main steps: construction of a curve network using Hermite curves, followed by a construction of corresponding cubic normals, and surface fitting using biquartic rectangular patches interpolation. The proposed scheme is based on the idea of - and compatibility conditions for the curve net and the rectangular patches, respectively.

1

Key words: surface, rectangular patch, Hermite curve, interpolation Cross Boundary Derivatives(CDBs).

Acknowledgment The work has been done in collaboration with Zhang Yongjie from the Computational Visualization Center (CCV). All computations presented in the work were done through the National Science Foundation’s National Partnership for Advanced Computational Infrastructure (NPACI).

Contents 1 Introduction

4

1.1

Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.2

The Report Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2 Analysis

6

2.1

The

Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.2

Why cubic curve? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.3

Why cubic normal? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3 Biquartic 3.1

13

Curve Network Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.1.1

Hermite Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.1.2

Cubic Normals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3.2

3.3

Surface

surface fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.2.1

Vertex Nodes and Mid-Edge Nodes Contributions . . . . . . . . . . . . . . 19

3.2.2

Middle node Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Modeling of a Human Head . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Other

Interpolation Methods

32

4.1

Inverse Distance Weighting (IDW) Interpolation . . . . . . . . . . . . . . . . . . . 33

4.2

Bicubic

Surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2

5 Conclusions and Future Works

43

A Singular Value Decomposition (SVD)

44

B Subroutines

49

3

1 Introduction In a smooth surface reconstruction, two main forms are widely used to interpolate a parametric representations and implicit representations.

surface:

There are various successful approaches of using implicit representation in interpolating of a surface. Bajaj and Xu pioneered implicit algebraic surface patches, A-patches [1, 2, 3, 4] (see also [5],[6]). A-patches are defined by a manifold triangulation in . This technique is based on the idea of zero contour of trivariate function. Other implicit schemes include B-patches [7], and Spatches [8] which also provide elegant solutions to the smoothing problem at the cost of a currently non standard patch representation.

The proposed surface reconstruction scheme is based on an explicit parametric representation of rectangular patches focusing an arbitrary unstructured surface mesh. The scheme has been implemented within our Geometric Modeling Package (GMP), fitting into a general class of both explicit and implicit parameterizations. The GMP is used to support geometrical representation for Finite Element (FE) discretization. Compared with implicit parameterizations, parametric representations have the advantage of simplicity and better hold on the effect of the geometric data (in our case -an unstructured surface grid of bilinear quadrilaterals with normals prescribed at the grid points) on the ultimate shape of the reconstructed surface.

A number of methods for constructing smooth parametric surfaces has been proposed. The splines assembled from B-splines are widely used, because they have explicit formulas and built-in continuity. However the tensor product B-spline representation has a major shorting, the B-spline mesh must be a structured mesh. Using non uniform rational B-splines (NURBS) does not solve this problem since the trimming destroys one of the chief advatages of the B-spline representation, its built-in smoothness so that one ends up with the tricky task of smoothly joining the trimmed pieces [9, 10]. J. Peter pioneered the idea of surface splines [11, 12, 13, 14] (see also[15, 16]), which devises a representation that removes the regularity restrictions at the cost of creating a refined mesh of quardrilateral subcells. The mesh refinement can separate irregular vertices by implementing initial mesh refinement, edge cutting and quadratic meshing. The refined mesh guarantees that each orginal vertex is surrounded by vertices of degree four. The approach is then reduced to the B-spline paradigm by constructing surface parametrization using quadratic patches. The scheme is in the same conceptual frame work as tensor-product B-splines, and local evaluates by averaging and obeys the convex hull property [17, 18]. Other works include generalized subdivision scheme, [19, 20] and reparametrization scheme [21, 22].

4

1.1 Motivation The -adaptive FE method is one of the most powerful methodologies for simulating complex engineering problems, especially those involving singular solutions, complicated geometry, or multiple media with varying properties. The present work is motivated by an application of methods to simulate the absorption and diffraction of EM-waves in the human head. The EM waves are emitted by a dipole antenna which models a mobile phone. The issue of a precise geometric modeling is especially sensitive in high accuracy simulations. In our previous work, we interfaced GMP with the geometric data provided by J. Zhang at CCV, to obtain the connectivity information, and construct a piecewise trilinear model of the human head [23]. The data provided by Zhang is generated from MRI scans; the head model is represented in terms of the following nodal information: Nodal connectivities for each hexahedron. Coordinates of each node in a physical frame. Normal vectors

for each node on the surface of the head.

Obviously, a linear geometry model is insufficient for high accuracy simulations. The hp-adaptive method requires at least a continuous geometry reconstruction. If we only use continuous geometry reconstruction, the FE code adapts meshes to resolve the non-physical scattering of waves on edges resulting from the poor, low order geometry representation. The term interpolation is used here to describe ways of fitting a curve or surface to a set of data points or curves [24]. The objective of this project is to reconstruct a surface that interpolates given vertices with prescribed normal vectors. The proposed surface reconstruction scheme is based on a parametrized surface by implementing a local interpolation on rectangular patches, and it has the following properties: local control of shape, numerical stability, smoothness and continuity, ability to evaluate derivatives. The reconstruction of free-form geometry is the main modeling challenge. Generally speaking, there is no restriction on the number of cells meeting at a mesh point or the number of edges adjacent to a mesh cell. The surface Mesh cells need not be planar. 5

1.2 The Report Structure This paper is structured as follows: Section 1 explains the motivation for developing the surface reconstruction schemes. Section 2 lists the main steps for the surface reconstruction scheme, establishes the compatibility conditions, and explains the reason for using Hermite curve and cubic normals. Section 3 elaborates on the processes of constructing a curve network and biquartic rectangular patches. Several examples are given to demonstrate the effectiveness of the scheme. Section 4 presents details on implementing two other alternative surface interpolation schemes: Inverse Distance Weighting (IDW) 4.1 and Bicubic surface Interpolation 4.2. Section 5 summarizes and Appendix A illustrates the Singular Value Decomposition (SVD) used in section 3. Appendix B lists routines developed in the course of this project.

2 Analysis 2.1 The Constraints

Consider a rectangular patch with parameters and its neighbor patch with parameters as illustrated in Fig. 1. Concentrate on the common boundary where . If both rectangular patches are sufficiently smooth, and are compatible if and only if the surface normal to is well defined and agrees with the normal of at each point of the boundary:

(2.1) guarantees non-vanishing normals. For surfaces, continuity means oriented tangent plane continuity [25], smoothness for surfaces differs from smoothness for

bivariate functions. In the proposed methodology, the two steps,

!#"$ %'&)(+* ,.-0/ !213"

surface reconstruction will be done in

Step 1: Construction of a curve network. We construct admissible parametric curves , which are continuous at the vertices. Continuity is achieved by using an alternative sufficient constraint that forces the mesh curves to interpolate vertex specified at each data while having compatible normals point on the boundary of a rectangular patch, i.e., along the mesh curves, see Fig.2.

! "4$ %'&5(6* 78,.-0/.9 1:"

Step 2: surface fitting. We interpolate between curves of the network obtained in Step 1 using a smooth parametrized surface, by implementing an algorithm for local interpolation of rectangular patches, see Fig.3. A rectangular patch is the image of a bivariate polynomial 6

X1( u , v ) v

v

1

u

w

0

1

X 2 ( u, w )

u u 1

0

1

w

Figure 1: Parametrization of two abutting rectangular patches

N( ξ )

n1

n2

X’( ξ ) p2

X( ξ )

p1

, and its normal vector , ! ! 1 " $ % '&)( * parametrized with parameters 9 / restricted to a standard domain, 9 " Figure 2: A parametric curve

.

The compatibility conditions for curves should guarantee the mity to the boundaries. They consist of two parts,

smoothness as well as confor-

compatibility conditions for curves:

, / 11 , / 11 , / & , /0&

(2.2)

compatibility conditions for normals:

/ !21 , / ! 1 7, / ! 1 7

(2.3)

In a similar way, the compatibility conditions for rectangular patches also consist of two parts,

/0& 1 1 /

compatibility conditions for patches:

/ ! ! 1 1 , )/ ! ! 1 1

/ 5& , '/

(2.4)

compatibility conditions for patches:

! 1 ! 1 ! / 7 , )/

1/ ! / ! '& 1 7 , / ! 1

/0& 1

/ 1 / 1

7, / ! 1

/ 1 1 0/ & ! ! 1 1 , '/ ! ! 1 1

/ / % , /

where is the normal function along the curve; tives (CBDs).

! 1 ! 1 ! / & 7 , '/

! 1 ! 1 ! /% 7 , / % - & are the Cross Boundary Deriva

ξ2 1

X ( ξ 1, ξ 2 ) 1

0

ξ1

n3 n4

N c3 ( ξ 1 ) N c4 ( ξ 2 ) (ξ ) X n 2 p4 c3 1 n X c4 ( ξ 2 ) 1 p 3 N c2 ( ξ 2 ) Xc2 ( ξ 2 )

p1

N c1 ( ξ 1 )

X c1( ξ 1 )

p2

Figure 3: Interpolation of a rectangular patch

2.2 Why cubic curve? We usually want the reconstructed curve to be as smooth as possible, i.e., minimize the wiggles. avoid using high-degree polynomials.

8

Cubic curves are commonly used in graphics because curves of lower order commonly have too little flexibility, while curves of higher order are usually considered unnecessarily complex and can easily introduce undesired wiggles. The use of cubic curve can also be explained in an optimization sense. As the second derivative of curve approximates curvature, we determine an optimal curve reconstruction scheme by seeking the solution to the following variational problem:

, / !21

,

Given endpoints and normal vector , find a curve parametrization which satisfies compatibility conditions in equation (2.2) and minimizes the mean (linearized) curvature,

$ %'&)( *

, / ! 1 ! * (2.5) / %'& 1 , and it admits a unique solution which The problem can be stated formally in space

satisfies the following (equivalent to the minimization problem) variational statement,

(2.6) , , ! % for every parametrization, test function , satisfying homogeneous essential boundary condi-

tions,

, / 1 1 , / 11 , /0& , / & 6

Integration by parts leads to

,

,

!

(2.7)

, , ! , , , ! , , , , ,

(2.8)

Restricting ourselves first to test functions that vanish on the boundary along with their first order derivatives, we get

, , ! This implies that , and, therefore, ,

variational statement reduces to,

, "

/ %'& 1

must be a cubic polynomial

, / & 1 , / & 1 , / 1 , / 1 9

(2.9)

, " .

The

, / & 1 , , / 1 , we conclude that (2.10) , / & 1 , / & 1 !21 . The direction of unit binorThe unit tangent vector on the curve is defined as /

! 1 !21 / ! 1 / ! 1 . mal is determined by / . The unit principal normal / 1 1 Decomposing , / & and , /0& into their tangent and principal normal components, , /0& 11 / , /0& 1 1 1 1 / , /0& 1 1 1 1 , /0& / , /0& / , / & (2.11) Assuming

we reduce equation(2.10) to,

/ , /0& 1 1 / , /0& 1 1 / , /0& 1 1 / , / & 1 1 (2.12) This, together with essential conditions for test function in equation (2.8), implies the final natural ! boundary condition at & , , /0& 1 % 1 or equivalently, , / & 6 . In the same way, , / 1 , / 1 (2.13)

The cubic curve reconstruction has the following properties: it minimizes the mean (linearized) curvature, it provides the lowest order polynomials representation that interpolates two points and allows for the gradient at each point to be defined which makes the continuity possible it reduces automatically to linear polynomials ( a straight line segment), for appropriate vertex data. The algebraic form of a parametric cubic curve is given by the following vector equation:

! , / !21 ! !

(2.14)

where are vector algebraic coefficients. The coefficients are not the most convenient way of controlling the shape of the curve in typical modeling situations, nor do they contribute much to an intuitive understanding of the curve. A practical alternative is offered by the Hermite 10

n 2(ξ )

n 1(ξ )

t1 ( ξ )

t 2 (ξ ) p2

p1 Figure 4: A Hermite curve

interpolation, which allows us to define a curve segment in terms of its endpoints [26]. More specifically, we need to know the coordinates and tangent vectors at both endpoints, see , we obtain the algebraic coefficients Fig.4. Evaluating equation(2.14) and its derivative at in terms of the boundary conditions by solving a set of four equations in four unknowns,

! %'&

(2.15) Substituting the coefficients into equation (2.14) and rearranging the terms produces, , / !21 / ! ! ! ! & !21 1 /! ! ! 1 ! 1 / !21 ! 1 / !21 !21 (2.16) / / / /

where

!21 ! 1 ! 1 ! / ! 1 / ! / !21 / /

! / !21 !21 ! / / 2! 1

are the Hermite Basis Functions. Differentiating equation(2.16), we get,

, / 2! 1

(2.17)

Matrix algebra and its notation scheme offer a compact mathematical form for representing a curve. We can rewrite equation (2.16) as a product of two matrices. in matrix form, where

is a

, / ! 1

$ ! ! ! 5& ( and

& & &

Hermite basis transformation matrix,

&

(2.18)

&

$ (

(2.19)

is the Hermite control matrix,

Equation (2.18) is also known as the cubic Hermite spline equation. 11

(2.20)

2.3 Why cubic normal?

Figure 5: Five basis functions and their first derivatives To obtain the desired conditions in equation (2.5), we need to find the normal derivative along the mesh curves. The approach does away with position vectors and tangent vectors . In other words, the curve reconstruction and , and uses an independent normal function defined along involves not only defining the curve but also specifying the variation of normals the curve and matching the normals at the endpoints,

7 , /! 1

78, / 1

7 ,

(2.21) 7 ,/&1

! ! 1 The parametrization for a rectangular patch 4/ is then requested to satisfy compatibility conditions (2.4) and compatibility conditions (2.5). The use of bicubic polynomials " which in turn led to geometrically unacceptable situation. implied 7, ! ! Biquartic functions involve linear combinations of twenty five monomials in and . The scalar-valued interpolation involves tensor products of five one dimensional shape functions, the !1 !21 Hermite basis / + & '' and a fifth bubble function / , illustrated in Fig.5. The order of the surface interpolation introduces four additional edge shape functions which result in a total 1 scalar unknowns to satisfy condition. Thus, we have five coefficients per of / edge. As Condition (4.84) must be enforced along each edge, the maximum polynomial order of 7 , / !21 is three, 7, " .

12

3 Biquartic

&

Surface

3.1 Curve Network Construction

, / ! 1

78, / ! 1

The curve network construction includes Hermite curves interpolation along with its cubic normal interpolation . The curve network should satisfy compatibility conditions in compatibility conditions in (2.3) for normals. The curve net will later be (2.2) for curves and extended to a smooth surface. 3.1.1 Hermite Curves

of polynomial (2.16). A Hermite curve is constructed by determining the vector coefficients The two unknown vectors can be expressed in terms of their Euclidean components, and . The boundary conditions resulted from the variational principle are,

/

1

/

1

essential boundary conditions:

, / 11

, / & %

natural boundary conditions:

, / 11 , /0&

(3.22)

(3.23)

where are two unknown scalars. equations (3.22) and (3.23) yield a linear system of equations to be solved for the components of , , and constants and .

, / !21"

The cubic curve may degenerate to a lower order polynomial. In the case of a straight line segment, , the vector coefficients corresponding to the third order and the second order terms vanish,

(3.24)

As the second derivatives of the curve vanishes,

, / !21 6

13

(3.25)

System (3.23) can now be solved for

and

,

which yields the final form for the curve,

/ ! 1 /0& ! 1

In this degenerated case, the nodal data of the curve bility conditions,

(3.26)

! , ,

(3.27) and

must satisfy the compati-

Assume now that the cubic curve degenerates to a second order polynomial vector coefficient corresponding to the third order term vanishes,

, / !21 "

(3.28)

, the

(3.29)

The second derivative of the curve is a constant different from zero, which implies that,

, / 1 , /0& 1

We have

/

(3.30)

1 /

1

(3.31)

This implies condition (3.28) from which, in turn, follows that the curve must degenerate to a line segment. Consequently, the cubic can never degenerate to a second order polynomial. 3.1.2 Cubic Normals The third order normals along the curve can be written as:

(3.32) / ! 1 / ! 1 / ! 1

!21 !21 where are two unknown vector coefficients, / and / are usual linear shape functions !21 and / + % are integrated Legendre polynomials. Note that, (3.33)

78, / ! 1 / !21

14

They can be expressed as,

! & !! ! /! ! /

with derivatives,

&

! /! & & 1 & ! & !

& 11 ! 1 & / &

! &

(3.34)

! 1 7 , / ! 1 is a fifth order !21 & , / points. As / / ! 1 is equivalent to enforcing 1 / & 1 % / (3.35) / 1 / 1 / 1 / 1 1 / & 1 have already been satisfied in (3.22). Taking / 1 , we have, Conditions / / 1 7, / 1 1 , / 1 1 / / / , / 1 / , / 1 1 / 1 1 Similarly, taking / , we have, 1 7 , / 1 , / 1 7 , / 1 , / 1 / / / 1 1 1 / 1 1 , / 1 / 1 1 / 1 1 / 1 , / 1 / , / / / , / / 1 Then, taking / , we have, 1 7 , / 1 , 1 1/ 1 1 27 , / 1 , / 1 1 1 1 / / , / / / , / / 1 Finally, taking / , we have, 1 1 , / 1 / & 27 , / & / , / 1 1 / 1

An degree polynomial fits a curve to polynomial, enforcing condition (2.3)

15

The above equations give a linear system of four equations in terms of six unknown components of . In a matrix form,

where

:

are

(3.36)

matrices of coefficients,

, / 11 , / 1 , /

and

, / 11 2 , / 1 , /

6

(3.37)

is a known right-hand side vector,

1 , / 1 / 1 1 1 1 / , / / , / Two additional needed conditions are constructed by minimizing the norm of the derivative 7 , , !1 ! 7 , / * (3.38)

Using the orthogonal property of integrated Legendre polynomials, we get

(3.39) / 1 / 1 / 1 and / 1 ! ; is subject to the constraints where (3.36). To solve the minimization problem by the method of Lagrange multipliers, we form a

system of ten equations in terms of the components of the two unknown vectors , and four 1 Lagrange multipliers 4/ ' . The matrix equation is (3.40)

where is a & : & square matrix of coefficients; is the solution vector of ten unknowns; is a known right-hand side vector; are two : diagonal matrices,

(3.41) 6

16

& : &'

The matrix may degenerate to a singular matrix. This happens, in particular, when the general cubic curve degenerates to a straight line segment, as discussed in the previous section. The last two constraints in (3.36) then vanish which results in singular matrices , , and the global matrix .

It is for this reason, that we can not employ the standard Gauss elimination to solve system (3.40), and use the Singular Value Decomposition (SVD) technique instead. The implementation based on SVD is uniformly stable for all possible situations including the degenerated and nearly degenerated matrices. Details on the SVD implementation with illustrating numerical experiments are provided in Appendix A.

1 /0& '& 1 / & % 1 /0& '& '& 1 1 /0& %'& 1 / / 3 / 1 / 3 ! 1 1 , / ! 1 7 ,/

To illustrated the resulted curve net, we show a special example of a spherical rectangular patch in Fig.6. The four vertices has position vectors

and the corresponding normal vectors are

. Applying the two steps discussed above, we obtain the curve net with a Hermite curve illustrated in blue and cubic normals along the curve illustrated in red, see Fig.6. The curve net will be used later to construct a rectangular patch.

3.2

surface fitting

,.- & 5' / ! ! 16"

Once we have constructed the four edge functions , along with the corresponding cubic normal vectors , we proceed now with a construction of parameterizations can be written for the rectangular patches. A general biquartic rectangular patch as the sum of vertex nodes contributions , mid-edge nodes contributions , and middle node contribution ,

78,.- & ''

/ ! ! 1 / ! ! 1 5/ ! ! 1 /! ! 1

(3.42)

The biquartic parameterization is constructed by enforcing the Compatibility Conditions (2.4) and the Compatibility Conditions (2.5) on the boundaries. The in expression (3.42) can be written as,

5/ ! ! 1

)/ ! ! 1 )/ ! ! 1

where is a vector coefficient for the middle node contribution, and shape functions, see Fig.15,

/! 1 /! 1 17

(3.43)

is the corresponding face

(3.44)

1.4 1.2

n1 = (0.5, 0.5, −0.5) p3 = (1, 1, 1 )

n = (0.5, −0.5, 0.5) 1 p = (1, 0, 1 )

1

4

0.8 0.6

Z

0.4 0.2 0 −0.2

p1 = (1, 0, 0 ) n1 = (0.5, −0.5, −0.5)

−0.4 0.5

p = (1, 1, 0 ) n22 = (0.5, 0.5, −0.5)

x1 1.5

0.2

0

−0.2

−0.4

y

0.4

0.6

1.2

1

0.8

p = (1, 1, 1 ) n31 = (0.5, 0.5, −0.5)

1.4 1.2

n1 = (0.5, −0.5, 0.5) p4 = (1, 0, 1 )

1

Z0.8 0.6

p2 = (1, 1, 0 ) n2 = (0.5, 0.5, −0.5)

0.4 0.2 0

p = (1, 0, 0 ) n11 = (0.5, −0.5, −0.5)

−0.2

1.5

−0.4 0.8

1

0.9

x

1

0.5 1.1

1.2

0 1.3

y

−0.5

Figure 6: A spherical rectangular patch

! ! 1 ) / /! ! 1

7, -

Note that the term vanishes along all four edges. In other words, does not affect the behavior of on the boundary. On the other side, the first two contributions in (3.42) are uniquely determined by the boundary data - edge functions and normals . This results in a two step procedure:

,.-

18

−3

x 10 4 3.5 3

ψs

2.5 2 1.5 1 0.5 0 1 0.8

1 0.6

ξ

2

0.8 0.6

0.4 0.4

0.2

ξ

0.2 0

1

0

Figure 7: Face shape function

/ ! ! 1 / ! ! 1 5/ ! ! 1

Step1: Construct a surface parameterization lating the boundary data, Step2: Determine vector

interpo-

in (3.43) using a minimum energy principle.

3.2.1 Vertex Nodes and Mid-Edge Nodes Contributions

n3 n4

n t p t n p t p p t 1

2 4

1 4

4

2

3

2 3

ξ2

ξ1

1

t

2 1

t

2

1 1

1 2

t22

Figure 8: Vertex Node Interpolation

/ ! ! 1

t13

"

The in (3.42) is a standard bicubic Hermite surface interpolant involves only vertex data, see Fig.8. The total vertex contribution can be expressed as,

/ ! ! 1

! ! 1 ! ! 1 ! ! 1 / / /

19

.

It

Here

denote the position vectors for each of the four vertices . are the corresponding bicubic vertex shape functions, listed in Table 1 and illustrated in Fig.9.

Vertex number

Table 1: The shape functions for position vector at vertex

!

& ''

&

is the vertex number, and are eight tangent vectors, where is the index of . For any rectangular patch, the eight tangent vectors are obtained from the curve functions ,

, , / 1 , / 1

, 0/ & 1 , / 1

, / & 1 , / & 1

, / 1 , / & 1

!

(3.45)

are the corresponding shape functions for tangent vectors at vertex in terms of , see Figures 10, 11,12 and 13. The eight tangent vector shape functions are listed in Table 2. Vertex number

Table 2: The shape functions for tangent vectors at vertex

are the four unknown twist vectors (mixed derivatives) at each vertex .

are the bicubic shape corresponding to each twist vectors, see Fig.14. The four twist vector shape functions are listed in Table 3. The

5/ ! ! 1

! ! 1 ! ! 1 5/ /

in (3.42) denotes the mid-edge nodes contributions,

20

(3.46)

1 0.8

1

ψ v1 0.6

0.8

0.4

ψ v2

0.2

0.6

0 1

0.4 0.8

0.2 0.6 0 1

ξ

0.4

2

0.8 0.2

0

0.1

0

0.3

0.2

0.4

0.6

0.5

0.7

0.9

0.8

1

0.8 0.6

0.4 0.4

0.2

ξ1

0.2 0

1

ξ1

0

1

0.8

ψ v3

1 0.6

ξ2

0.8

ψ v4

0.6

0.6

0.4

0.4

0.2

0.2

0 1

0 1 0.8

0.8

1

ξ

2

0.6

0.8 0.4

0.6

0

0.4

0.2

0.2 0

ξ

0.2 0

0

1

Figure 9: Vertex shape functions Vertex number

,

0.8 0.6

0.4

2

0.4

0.2

1

ξ

0.6

ξ

1

Table 3: The shape functions for twist vector at vertex

where are eight vector coefficients to be determined, two per edge; shape functions, see Table 4.

are the corresponding

The Compatibility Conditions (2.4) and the fact that the curves have been reconstructed using cubic polynomials only, , imply that contributions corresponding to last four shape functions . Note that the condition must simply vanish, i.e.,

/! ! 1

, / ! 1 "

8% & 5'

21

0.16

0.16

0.14

0.14

0.12

ψt11

0.12

ψt12

0.1

0.1

0.08

0.08

0.06

0.06

0.04

0.04

0.02

0.02

0 1

0 1 0.8

0.8

1 0.6

ξ

1 0.6

0.8 0.4

0.2

2

0.8

0.6

0.4 0.2 0

0

ξ

ξ1

0.6

0.4 0.4

0.2

2

0.2 0

0

ξ1

Figure 10: Shape functions for tangent vectors at vertex 1

0

0.16

−0.02

0.14

−0.04

ψt21

0.12

ψt22

−0.06

0.1

−0.08

0.08

−0.1

0.06

−0.12

0.04

−0.14

0.02

−0.16 1

0 1 0.8

0.8

1 0.6

ξ2

1 0.6

0.8 0.6

0.4 0.4

0.2

0.2 0

0

ξ2

ξ1

0.8 0.6

0.4 0.4

0.2

0.2 0

0

ξ1

Figure 11: Shape functions for tangent vectors at vertex 2

/! ! 1

does not apply to the contributions of the first four shape functions which automatically vanish on the patch boundary and contribute only with non-zero normals, see Fig.15.

The Compatibility Conditions (2.4) require the knowledge of Cross Boundary Derivatives (CBDs) along the boundary. The CBDs have a crucial effect on the shape of the constructed patches; they allow for the patches to effectively reflect the variation of the normals . The CBD at any point on the patch boundary is perpendicular to both the corresponding normal and the tangent vectors. Given the reference coordinates of a point on the rectangular patch, we first identity four corresponding points on the patch edges with coordinates . We compute then the four CBD functions corresponding to each of

7 ,.-

/! ! 1

/ ! 1 /0& ! 1 / ! 5& 1 / % ! 1

22

0

0

ψt31

−0.05

ψt32

−0.05

−0.1

−0.1

−0.15

−0.15

−0.2 1

−0.2 1 1

0.8

1

0.8 0.8

0.6

0.8 0.6

0.6 0.4

ξ

0.2

2

0.6 0.4

0.4

ξ1

0.2 0

0.4

ξ

0.2

2

ξ1

0.2 0

0

0

Figure 12: Shape functions for tangent vectors at vertex 3

0.16 0 0.14 −0.02 0.12 −0.04

ψt41

0.1

ψt42

0.08

−0.06 −0.08

0.06 −0.1 0.04 −0.12 0.02 −0.14 0 1

−0.16 1

0.8

0.8 0.6

1 0.6

ξ2

1 0.6

0.8

0.4 0.2

0.4 0

0.2

ξ2

ξ1

0.8 0.6

0.4 0.4

0.2

0.2 0

0

0

ξ1

Figure 13: Shape functions for tangent vectors at vertex 4 the four edges of the patch:

/! 1 /! 1 / & ! 1 / ! 1 / ! '& 1 / ! 1 / % ! 1 / ! 1

/ !! 11 / ! 1 / ! 1 /

/! 1 /! 1 /! 1 /! 1 /! 1 /! 1 / ! 1 / ! 1

78, "

! 1 / !1

/! 1

/! 1

/

" -

(3.47)

Note that CBDs are fourth order polynomials along the edges and the normals along the curve are third order polynomials . In order to satisfy the compatibility conditions, the rectangular patch must satisfy the system of four equations in (2.5). With , the 23

! ! !

0.025

0

0.02

ψc2

ψc1 0.015

−0.005 −0.01 −0.015

0.01

−0.02

0.005

−0.025 1 0.9 0.8

0 1

1

0.7 0.6 0.8

0.5

1 0.6

ξ

0.8

0.4

0.2

2

0.2 0

0.3

ξ2

ξ

1

0

0.6

0.4

0.8 0.6

0.4

0.4 0.2

ξ

0.2

0.1

1

0

0

0.025

0.02

0

ψc4

ψc3 0.015

−0.005 −0.01 −0.015

0.01

−0.02 −0.025 1

0.005

1

0.8

0 1

0.8 0.6

0.8

1 0.6

ξ2

0.6

0.8

0.4

0.4

0.6

0.4 0.4

0.2

0.2 0

0

ξ2

ξ1

0.2 0

Figure 14: Shape functions for twist vectors

Edge number

0.2

ξ1

0

,

Table 4: The shape functions for edge

! 1 !21 !/ 1 ! / ! 1 7 , / ! 1 ! 1 ! / ! 1

7 ,/ ! /

7 ,/ , %

system of equations can be expressed as,

,

(3.48)

where and are functions in terms of four unknown twist vectors and four unknown vector edge coefficients , respectively. The parametric function for vertex contribution in equation

24

0.01 0 0.008

ψe21

ψe11

−0.002 −0.004

0.006

−0.006 0.004 −0.008 −0.01 1

0.002

0.8

0 1

1 0.6

0.8

0.8

1 0.6

0.8 0.6

0.4

ξ

0.6

0.4

0.4

0.2

2

0

ξ2

ξ

0.2

1

0

0.4

0.2

ξ

0.2 0

1

0

0.01

0

0.008

ψe31

ψe41 0.006

−0.005

0.004 −0.01 1

0.002 0.9 0.8 1

0.7 0.6

0 1

0.8 0.5

0.8 0.3

ξ2

1

0.6

0.4

0.6

0.8

0.4 0.2 0

ξ2

ξ1

0.2

0.1

0.6

0.4 0.4

0.2

0.2 0

0

Figure 15: Edge shape functions (3.45) can be reduced to a matrix form,

0

/ ! ! 1 $ / ! 1 / ! 1 / ! 1 / ! 1 ( $ / ! 1 / ! 1 / ! 1 / ! 1 (

We have,

/

7 , )/ !! 11 11 7 , /! 11 , //

, /! 11 / 7

7 , /

/

7 , )/ !! 11 11 / 7 , / ! 1 1 / 7 , / ! 1 1 / 7 , /

ξ1

(3.49)

$ / ! 1 / ! 1 ( +

(3.50) 25

and,

/ , / !2!211 11 //

,, '5// 2! 1 1 / ! 1 / , '/ !21 1

(3.51)

is a matrix prescribed in terms of tangent vectors

/

7 , / !! 11 11 7 , '/ ! 1 1 //

7 , /

/ 7 , / ! 1 1

/! 1

/ !21

,

/

7 , )/ !2!211 11 / 7 , / !21 1 / 7 , / !21 1 / 7 , /

$ / ! 1 / ! 1 (

(3.52)

/ !216"

As is a seventh order polynomial , vanishing at the endpoints of the edge, enforc ing condition (2.5) is equivalent to enforcing,

(3.53) / 1 % & & 1 equations, we get values of the components of with . Solving a system of / & eight vector coefficients and + & '' .

3.2.2 Middle node Contribution Mathematical formulation of any boundary value problem consists of a differential equation and boundary conditions. The connection between transfinite interpolation and the boundary values problems is explicit in the last section. In this case, we exactly interpolate the prescribed boundary conditions. However, the behavior of the interpolant away from the boundaries is quite arbitrary. Thus, the solution of a boundary value problem can be viewed as the construction of a function that extends the boundary conditions into the domain, with differential equations playing the role of a constraining or smoothing operator. The interpolation problem has no unique solution; there are infinitely many functions interpolating any given data. An unique construction can be established by putting additional constraints on the interpolant. Often such constraints appear as the minimization of some quantity. Many interpolation schemes use minimization of energy [27] as a means for controlling the shape of the interpolant,

/ / ! ! 1 1 *

26

(3.54)

Figure 16: A hexahedron with two adjacent faces on a spherical surface

where

,

is the rectangular reference domain, and

/ ! ! 1 '/ ! ! 1 5/ ! ! 1 (3.55) ! ! 1 interpolating the boundary conditions given in equation (3.42); with the first term / unknown for the middle node; are shape function described be : / 1 is the vector fore. The value of the coefficient is obtained by minimizing (3.54). Upon differentiating (3.55) with respect to & '' , we construct and solve a system of three linear equations for components of , (3.56) / 1

Now we conclude the discussion with some simple examples. Using the existing curve net for a spherical patch in Fig.6, and applying the surface fitting technique above, we can get a hexahedron with two adjacent faces on a spherical surface, illustrated in Fig.16. In the similar way, we can obtain a hexahedron with two adjacent faces on a cylindrical surface, see Fig.17.

3.3 Modeling of a Human Head The discussed continuous geometry reconstruction technique has been applied to model the geometry of a human head, necessary for high accuracy 3D hp FE simulation. The implementation 27

Figure 17: A hexahedron with two adjacent faces on a cylindrical surface has been done within our Geometrical Modeling Package (GMP) [23] interfaced with software LBIE - Mesh Level Set Boundary and Interior-Exterior Mesher, developed at Center for Computational Visualization (CCV) at ICES. Given an MRI scan of a human head, we proceed in the following steps [28, 29], Pre-Processing. We first use the anisotropic diffusion method coupled with bilateral prefiltering to remove noise from imaging data. Depending on the application, suitable isosurfaces are selected for the interval volume by using the contour spectrum and the contour tree. Hexahedral Meshing. We then begin to extract 3D meshes from the interval volume, and a feature sensitive error function is adopted to reduce the number of elements while preserving features. Quality Improvement. Finally, the edge contraction method is then used to improve the mesh quality. The resulting trilinear hex mesh provides a minimum coarse representation of the head topology and a geometrical information on the surface of the head. The geometry information is given by labeling the surface nodes and specifying the corresponding normals. The topological (hex to points connectivities), and geometrical (normals)data is then imported into GMP, where the actual continuous geometry reconstruction takes place.

28

η3

5

8

●

8

●

ξs 2 ξs 1

5

(2)

6

7

6

7

9 12 (6) ξs 2

(5)

ξs 2

ξs 2

(3)

11 10

ξs 1

(4) ●

ξs 1

1

ξs 1

●

ξs 1

ξs 2 ξs 1

ξs 2

2●

4

1 (1)

4

η2

3

●

2

3

η1

Figure 18: A hexahedron in the reference frame The recently rewritten and upgraded GMP [30] supports the construction of exact parameteri zations for a general class of 2D (BEM) and 3D (FEM) manifolds in [31, 30]. We have explored a number of novel geometric modeling techniques and implemented them in the GMP. The GMP not only supports the construction of exact parameterizations for a general class of 2D [32] and 3D [33] objects, but also provides the derivatives of the mappings with respect to reference coordinates for any given points in reference frame. In our geometric modeling, a 2D object is represented as a union of curvilinear triangles or rectangles, while a 3D object is represented with a FE-like mesh of curvilinear hexahedral blocks. The geometry of the object is described then by constructing parameterizations for each of the blocks. The GMP model later can be used to generate the actual FE meshes of arbitrary high order, and to support geometry updates during mesh refinements. Each of the local edges or faces in the GMP has its own global orientation. Adjusting edge and face parameterizations for orientations involves transforming local edge and face coordinates into the global coordinates. We must ensure the compatibility of parameterizations. For example, parameterizations for adjacent hexahedra must be compatible with each other. Roughly speaking, when a face is shared by two hexahedra, and two parameterizations are used, the same FE on both sides of the face mesh must be obtained.

29

The coordinate transformations are handled in a hierarchical manner. First the coordinates of points are set. Next, the curve parameterizations are constructed in such a way that they conform to the existing vertex coordinates. Next, the parameterizations for rectangles must conform to the parameterization of the curves. Finally, in the most complicated case, the parameterizations for hexahedra must conform to the parameterizations of the rectangles, see Fig.19. The construction of the human head model is based on two GMP parameterizations: Transfinite Interpolation Rectangles and Transfinite Interpolation Hexahedra, both based on the classical transfinite interpolation and linear blending functions technique [34, 35]. The interface between ζ2

ζ1

x3 x2 x1

xH

xR η3

ζ2

η2

η1 η3

ζ1

ξ2 η2

ξ (ζ)

ηf

ξ1

η1

Figure 19: Compatibility of parameterizations for a hexahedron GMP and geometric data extracted from an MRI scan [36] provides connectivities of for a trilinear hex-mesh and topology of the boundary surface. Using this information, we then reconstruct a curvilinear hex-mesh with a continuous representation of the surface. The obtained hex-mesh head model is then used to generate the actual meshes for hp-Adaptive FE simulations. We proceed in the following steps: Step1: Linear Hex-Mesh Model Reconstruction Before generating a full nonlinear model, it is convenient to check the topology of the GMP model by reproducing the piecewise trilinear mesh within the GMP using the simplest trilinear parametrization for the GMP hexahedra. Note that all twelve edges of each hexahedron are straight line segments connecting two adjacent nodes. The geometric data required by constructing the linear model includes the indexes and coordinates of the eight nodes for each hexahedron. Contrary to the minimal data provided in the CCV model, the GMP stores all topological connectivities in a hierarchical 30

manner: points to edges and edges to points; edges to faces and faces to edges; faces to hexas and hexas to faces. This extended connectivity information is generated through a global sequential search, processing one hexahedron at a time. First the local entities are connected to the corresponding global entities, i.e., vertices to points, edges to curves, and faces to rectangles, and then the overall global entities are updated[23]. The topological connectivity information is sufficient to generate the simplest geometry model based on linear curves, bilinear faces, and trilinear hexahedron parameterizations. Step2: Obtaining Topological Information on the Boundary Surface. The necessary topological information for the construction of curvilinear hexahedra is collected. The geometric data provided by CCV includes labeling the nodes on the head surface and providing the outward normals at the nodes. We loop through the GMP rectangles, and specify those with all four vertices on the boundary surface. The rectangles are redefined as a new rectangle type : biquartic rectangle (biqua), a rectangular patch on continuous surface. Rectangles with at least one vertex on the boundary surface are redefined as: transfinite interpolation rectangles (TraQua). The remaining rectangles are still stored as bilinear quadrilaterals( BilQua). In the same way, we designate those curves with both vertices on the surface as a new curve type: Points with Coordinates and Normals (CoorNrm). The remaining curves are marked with the label of a segment of a straight line (Seglin). All points are classified as regular points (Regular) and all hexahedra are of GMP hex type: transfinite interpolation hexahedra (TriLiHex). Step3: Curvilinear Hex-Mesh Model Reconstruction. The surface reconstruction scheme discussed in the last few sections is performed. A cubic curve-net with compatible cubic normals along the curves CoorNrm is first constructed, and then the surface fitting techniques is used to generate the biqua rectangles. Before generating the whole human head model, we first test the scheme on some feature parts of the head, e.g., the nose. Any hexahedron that constitutes a part of the nose is a special case because all its eight vertices are on the surface. Following the three steps for the linear nose model, see the first figure in Fig. 20(a), we obtain a curvilinear model, illustrated in in Fig. 20(b). Fig.21 displays a zoom on the bottom hexahedron on the nose for which the curvilinear geometry is the lest visible. The entire reconstructed human head model is presented in Fig.22, with a statistics of the model summarized in able 5. The motivation for this work is to eventually simulate EM waves in a human head. This involves enclosing the head within a truncating sphere, and meshing the entire volume within the 31

Figure 20: The nose model of the human head. NRSURFS NRPOINT 0 25744

NRCURVE 8448

NRTRIAN 0

NRRECTA 4424

NRPRISM 0

NRHEXAS 704

Table 5: Numbers of global entities in GMP for the human head model sphere, and the head. Special absorbing boundary conditions are imposed on the truncating sphere to model the interaction with the rest of the space. We will use the model later to generate the actual meshes for -Adaptive Finite Element (FE) simulation of Electromagnetic (EM) wave.

4 Other

&

Interpolation Methods

In this section, for a comparison, we report on two other (less successful) niques implemented earlier. 32

reconstruction tech-

Figure 21: The bottom hexahedron of the nose model

4.1 Inverse Distance Weighting (IDW) Interpolation In this scheme, we use one of the most common techniques for surface interpolation, Inverse Distance Weighting Interpolation (IDW) [37]. IDW assumes that the value of the interpolant should be influenced more by the value of the nearby points and less by more distant points. The main advantage of the inverse distance weighting method is that it does not place any restriction on the topology of the sets being interpolated. The starting point for IDW interpolation is a set of , representing a distance from specific so called normalized distance functions geometrical object e.g., points or lines. In our case, the interpolation takes place over the reference rectangle, and four normalized distance functions will simply represent the distance from the rectangle’s edges.

/ 1 &

1 & , we define the corre / 1 / + & 1 as: 1 - / (4.57) 1 - /

Given any set of the normalized distance functions sponding inverse distance weighting functions

/

/ 1

1

/

is the normalized distance to the point. The exponents where interpolating functions as follows: 1. 2.

&

&

: the interpolant is not differentiable points

: the interpolant is differentiable

33

&

where

/

times at the points.

control behavior of the

1

Figure 22: The head model as a union of curvilinear hexahedra

/ 1

The inverse distance weighting functions form a partition of unity in the sense that . This partition of unity also assures completeness of the system and guarantees reproduction of constant functions.

&

, / !21

7, / ! 1

We proceed now with the discussion of the actual interpolation technique. As in the previous and normals defined along the sections, we start with a reconstruction of curves curves. The curve reconstruction is once again based on Hermite representation (2.16). To obtain the two tangent vectors at both endpoints, we use the geometric information on the adjacent endpoints. We select for the tangent vector to be a vector perpendicular to lying in the plane spanned by and vector . A simple geometrical manipulation yields:

where data .

/ / / /

11

11

(4.58)

is directed toward point . Note that is directed toward point and that in case of , reduces to . In the same way, we obtain the tangent vector at point

34

Figure 23: The cross boundary derivatives (CBD’s).

!1 7 , / 7 , / 2! 1 8 7 , / ! 1 !1 !21 ! 2 1 , / , / (4.59) 7 , / , / 2! 1 , / ! 1 , / ! 1 , / ! 1 !21 and projecting it onto the plane perpendicular to the tangent vector , / , !, 21 7 !21 7 / , / (4.60) 7 , / ! 1 !21 7 , / !21 7 , / 2 ! 1 We emphasize that normals 78, / are not polynomials, in fact, they are not even rational functions.

As the IDW interpolation produces non-polynomial parametrizations anyway, we are not bounded to a polynomial variation of the normals , and we can use geometrically simpler techniques. We construct the normal vector along the curve by taking a linear combination of normal vectors and ,

Once we have selected the normals, we proceed now with the patch interpolation. The construction is based on the selection of appropriate Cross Boundary Derivative (CBD) vectors. For each of the four edge points, we choose the CBD vectors to be the cross products of the corresponding tangent and normal vectors at the point, as shown in Fig.23. For a point on the first curve, we have,

/! 1

, /! 1

7 , / ! 1 , / ! 1

(4.61)

In a similar way, we can get the CBDs for the other three curves as follows,

, '/ !! 11

7 , '/ !! 11 , , /! 1

7 , / ! 1 , , /

7 , / ,

35

/ !! 11 /! 1 /

(4.62)

, - : & ''

We have four edge functions , along with the corresponding CBD vectors . We use a IWD Interpolation to derive a parameterization for the patches. The . The four parametrization is constructed over the reference domain, with

,.- & ''

9 /! ! 1

! ! " $ % '&)(

Figure 24: Normalized distance functions in reference frame normalized distance functions represent the distance from the four edges, which in turn are in this case can be described on boundaries shown in Fig.24. Each portion of the boundary is represented implicitly by the functions as follows:

/ !! !! 11 ! 6 !

/ & 6

/ !! !! 11 8 !& !

/ In order to guarantee that the interpolants have first derivatives, we take in equation (4.57). 1 The corresponding weighting functions / & '' , can now be expressed in terms of 9, ! /0& ! 1 /0& ! 1 /! ! 1 ! / & ! 1 /0& ! 1 ! ! / & ! 1 ! ! /0& ! 1 ! /0& ! 1 / & ! 1 ! ! 1 The interpolating function / for the rectangular patch is constructed as a linear com bination of edge functions , - & 5' with the corresponding weight functions & 5' , illustrated in Fig. 25, ! ! 1 (4.63) / ,.-0/ ! ! 1 / ! ! 1

The selection of edge functions

,.- &

36

ensures

continuity and it is done in the

following way,

, / !! !! 11 , '/ ! ! 1 , / ! ! 1 , /

,.- : & ! ''

/ !! 11' ! , '/ ! 1 / & , / ! 1 '! / & , / ,

, / ! ! 1

, ! ! ,

, are the corresponding

, / ! 1 , / ! 1'! , / ! 1 , / ! 1 , / ! 1

(4.65)

, / ! 1 '! ! 1 !/ ! 1 ,! / ! 1 , / ! ! !/ ! 1 ! 1 ! 1 , / , /

Similarly for derivatives,

(4.64)

,.- & ''

where , are the Hermite curves, and CBDs. With point approaching the first edge, we have,

! 1 ! 1

, / !! 11 , 5/ ! 1 , / ! 1 , '/

(4.66)

The limitation behavior of the edge functions ensures continuity of the patch parametrization across edges shared by two patches. Consider the value of the interpolant of the rectangular patch as point approaches the first edge,

/! ! 1

9

! ! 1 / , / ! ! 1 , & , , / ! ! 1 ,

/! ! 1

/! !

/! 1

1

,

,

(4.67)

Consequently, two adjacent rectangular patches have the same function values at the common edge, the interpolant approaches the physical coordinates of the boundary edge when approaches the corresponding reference coordinates of that edge. The derivatives of with

9

37

Figure 25: Physical domain defined on a reference domain respect to

!

!

and

behave in a similar way,

/! ! 1 !

/! ! 1 !

,.-0/ ! ! 1 !

/! ! 1 ! ! , & ! , ! , ! , / ! ! 1 ,! ! , / ! ! 1 , / ! 1

! ,

(4.68)

Note that IDW interpolation does not place any restriction on the topology of the curves being interpolated. We exam this scheme with the special case depicted in Fig.21. Interpolating two rectangular patches over the implicitly defined curves, we get a curvilinear hexahedron with its three adjacent faces on a surface, see Fig.26. The IDW interpolation ensures continuity along the edges but, in general, guarantees only a -continuity at vertices. The main disadvantage of using the inverse weighting functions is their loss of regularity at vertices. At each vertex of the reference domain, the inverse weighting function is discontinuous and its derivative blow up to infinity, see Fig.27. It is not clear how this lack of regularity will affect the mesh generation and quality of meshes. 38

Figure 26: A

4.2 Bicubic

continuous hexahedron

Surface

/! ! 1 " 78, / ! 1

In this section, we report shortly on an earlier attempt of rectangular patch reconstruction using . With bicubic polynomial reconstruction, we have only bicubic polynomials, three free parameters per edge (the twist vectors defined at vertices). Consequently, enforcing to the curve tangents is possible only if the degree of and orthogonality of normals do not exceed four. In the presented attempt, we assume that we start with a simple linear interpolation of the normal,

7 , / !21

, / 2! 1

7 , / !21 /0& !21 ! (4.69) !21 and only then follow with a reconstruction of the curve , / done in such a way that the orthog!21 7 , / ! 1 , is satisfied. The order in which we reconstruct curves !21 , / onality condition / , / ! 1 and normals 7, / ! 1 is then opposite to the biquartic reconstruction scheme, first the normals

and then the curves.

Using the same Hermite curve representation (2.16), we have to determine two unknown vector coefficients and . In order to ensure compatibility condition (2.2), we need

39

(4.70)

The first weight function W1(ξ1,ξ2)

The second weight function W2(ξ1,ξ2)

1

1

0.8

0.8

0.6

0.6

W2

W1 0.4

0.4

0.2

0.2

0 1

0 1 0.8

0.8

1

ξ 0.6 2

0.4

0.4 ξ

0.2

2

0.8 0.4

1

0.2 0

1

ξ 0.6

0.8 0.6

0.4 ξ

0.2

The third weight function W3(ξ1,ξ2)

1

0.2 0

0

0.6

0

The fourth weight function W4(ξ1,ξ2)

1

1

0.8

0.8

0.6

0.6

W3

W4

0.4

0.4

0.2

0.2

0 1

0 1 0.8

0.8

1

ξ20.6

0.4 ξ

0.2

0.2 0

1

ξ20.6

0.8 0.4

0.6

1

0.8 0.4

0.4 ξ

0.2

0.2 0

0

0.6

1

0

Figure 27: Weight functions of IDW The scalar function

/! 1

/! 1

Note that to enforcing,

/ !21

can be written as a function of unknown vectors

and

as:

/ 11 / 11 ! $ / ! 1 : '/ !21 / ! 1 / !21 ( / 1 / 1 & ! (4.71) / 1 / 1 / / " as ./ !21 " and 7 / !21 " . Enforcing / ! 1 (2.3) is thus equivalent

"

/ & 1 / 1

(4.72)

Unknown vector coefficients can be written uniquely as a combination of the three linear independent basis vectors. The concept of the dual basis enables the decomposition of a 40

vector in a nonorthogonal basis. Looking at equation (4.70), we construct three basis vectors as:

(4.73)

and are not necessarily orthogonal, can not be given by the inner

Because basis vector product of and . In

order to calculate the corresponding components, we introduce another set of basis vectors , called the dual of , defined by the orthogonality relation where delta. In terms of the dual basis , . The components of the is the Kronecker

4/ 1 1 4/

vector can then be calculated as

. This implies that,

4/ 4/

11

(4.74)

The two unknown vector coefficients can now be written in terms of the dual basis,

(4.75)

(4.76) / ! 1 / !21 / !21 / ! 1 % with, / ! 1 / 11 ! / !2! 11 !21 ! 1 (4.77) / / 1 /0& ! 1 / !21 1 ! !21 ! 1 / / /0& / !21 / 1 ! / !21 1 ! 1 / /0& : '/ / / (4.78) 1 1 Requesting / and / , we get a linear system of equations, / 11 / 11 / 11 % (4.79) / / / and . Finally, the remaining degrees of freedom of the curve and solved for components the energy function, are determined by minimizing / 1 , ! !21 / ! 1 / 1 / ! 1 / 1 !21 !

/ / * (4.80) Consequently,

41

Taking derivatives of

gives: 1 ! 1 1 ! 1 1 / !21 / !21 / / / / / 1 / !21 / !21 / 1 / ! 1 / ! 1 !

1 !21 !21 / 1 / ! 1 1 / ! 1 / 1 / !21 / !21 / / 1 / ! / 1 / ! / 1 !

in terms of the unknowns

Solving the resulting system of linear equations yields,

/

1 / 1

/

1

(4.81)

(4.82)

curve under the condition that vectors The prescribed procedure results in a uniquely defined and , used to defined the original basis vectors , are linearly independent.

The bicubic rectangular patch can be written as a Hermite patch in equation (3.50), with the four CBDs computed in the standard way,

/ ! 1 / ! 1 / ! 1 / ! 1 (4.83) / ! 1 / ! 1 /! 1 /! 1 Note that vector CBDs are also third order polynomials . In order to interpolate a surface, !21

/ ! 1 / ! 1 / & ! 1 / ! 1 / ! 5& 1 / ! 1 / % ! 1 / ! 1

: / !! 11 : / ! 1 : / ! 1 /

the rectangular patch must satisfy the system of four equations in (2.5). Introducing product of CBDs and its corresponding normals,

/ !21

/

7 , / 7 , / 7 , / 7 , / 7 , / 7 , / 7 , / 7 ,

/

7 , / !! 11 11 /

7 , '/ ! 1 1 / 7 , / ! 1 1 / 7 , / !21 1 / 7 , / !21 1 / 7 , / !21 1 / 7 , / !21 1 / 7 , /

/ !2!211 11 / !21 1 $ / !21 1 / !21 1 / !21 1 / !21 1 $ / !21 1 /

/ ! 1 / ! 1 (

/

! 1 /! 1 (

"

/0& 1 / 1 / 1

is a dot

(4.84)

We reduce the compatibility condition for the rectangular patch to equation is a fourth order polynomial, , it is sufficient to enforcing,

/

42

/ ! 1 . Since / ! 1

(4.85)

8& ''

Solving the resulting system of three equations for four edges (twelve equations in total), we obtain . the components of the four twist vector The bicubic scheme produces results similar to the examples of the spherical and cylindrical rectangular patches shown in Fig.16 and Fig.17. The major drawback of the presented construction lies in a possibility of a degenerated case in which the procedure fails. This happens when normal vectors and are parallel to each other and, at the same time, are not orthogonal to . In the case of , the minimization principal can be used to yield a straight line segment. It was the degenerated case that demonstrated an insufficiency of linear variation of and prompted us to use the biquartic parameterizations.

7 , / !21

5 Conclusions and Future Works The paper presents results of a preliminary study on geometric reconstruction in context of geometries reproduced from MRI scans and mesh generation for high order FE discretization. The presented biquartic scheme seems to be the lowest order continuity construction for general unstructured meshes. The polynomial parametrizations are inexpensive to compute and guarantee computations. In the case of the presented high regularity of parametrization necessary in IWD technique, it may not not satisfied. It is not clear at this point, however, how the ( and not ) regular parametrization will affect the convergence rates of high order methods. The important property of the presented reconstruction scheme is that it remain uniformly stable in the case of degenerated geometrical data (see the Appendix). Among other tasks, we intend also to collaborate with CCV on multi-resolution techniques and hierarchical geometry reconstruction schemes. At this point, the information on geometry contained in the original fine mesh reconstruction, during the coarsing stage is reduced to normals only. Ideally, the geometry reconstruction on the coarse grid should conform to the fine grid representation in a more elaborate, multi-resolution model. We intend to address these topics in our future work.

43

A Singular Value Decomposition (SVD) n1 = (0, 0, 1)

p1 = (0, 0, 0)

p2 = (0, 1, 0)

0.04

n2 = (1, 0.01, 0) n2 = (1, 0.05, 0) n2 = (1, 0.1, 0) n2 = (1, 0. 05, 0) n2 = (1, 1, 0) 0

0.02

z 0

0.02

−0.02 0.04 −0.04 0.06

0

x

0.2 0.4

y

0.08

0.6 0.8 0.1

1

n1 = (0, 0, 1)

p1 = (0, 0, 0)

p2 = (0, 1, 0)

0.04 0.03 0.02 0.01

z

0

−0.01

n2 = n (1, =1,(1, 0) 0. 05, 0) 2 n = (1, 0.1, 0) 2 n = (1, 0.05, 0) 2 n2 = (1, 0.01, 0)

−0.02 −0.03 −0.04 0 0.2 0.4

y 0.6 0.8 1 0.1

0.08

0.09

Figure 28:

0.07

0.06

0.05

x

0.04

0.03

0.02

0.01

0

varies from (1, 0.01, 0) to (1, 1, 0)

A very powerful set of techniques dealing with sets of equations or matrices that are either 44

singular or numerically very close to singular is the so-called singular value decomposition (SVD) [38, 39, 40, 41]. SVD allows one to diagnose problems related to solution of linear systems provides numerical a answer as well.

matrix can be written as the product of an column-orthogonal Any matrix , an diagonal matrix with positive or zero elements, and the transpose of an orthogonal matrix :

where

The diagonal elements of matrix

(A.86)

(A.87)

are the singular values of matrix

;

& (A.88) is a square matrix, it is also row-orthogonal, & . Equation (3.40) defines Since " & , as a linear mapping from the vector space into itself. If is singular, , called the null space / 1 , mapped to zero, / / 1 1 . then there is some subspace of 1 is called The dimension of the null space is call the nullity of A. The dimension of the range / 1 / 1 . the rank of . The relevant theorem is “rank plus nullity equals N”, / 1 and / 1 . Specifically, the columns of SVD explicitly constructs orthonormal bases for / set of basis vectors that whose same-numbered elements are nonzero, are an orthonormal span the range; the columns of whose same-numbered elements are zero, are an orthonormal basis for the null space. " , one has to: To find the singular value decomposition [42] of the matrix

1. Find the eigenvalues of the matrix

and arrange them in descending order.

2. Find the number of nonzero eigenvalues of the matrix

.

3. Find the orthogonal eigenvectors of the matrix corresponding to the obtained eigenvalues, and arrange them in the same order to form the column-vectors of the matrix .

"

45

"

4. Form a diagonal matrix of first eigenvalues of the matrix

placing on the leading diagonal of it the square roots got in step 1 in the descending order.

5. Find the first column-vectors of the matrix 6. Add to the matrix process.

the rest of

" ,

& ''

vectors using the Gram-Schmidt orthogonalization

&'

In our case, matrices and are square , so they are orthogonal and their inverse matrices are equal to their transposes. The inverse matrix of becomes,

4/ 1 / 1

(A.89)

If one or more are zero or very close to zero, matrix is singular or numerically singular in this case. Equation (3.40) can then be used to obtain the solution:

/ 1

Where the diagonal elements of matrix

(A.90)

are given by

&

(A.91)

which being a singularity threshold.

/ 1

In the case when is singular, we need to check whether the vector on the right hand side lies in . If it does, the singular set of equations has multiple solutions. If it does not, the problem (3.40) has no solution. When ’s are very small but nonzero, the matrix is ill-conditioned. In that case, the direct solution methods of decomposition or Gaussian elimination may actually give a formal solution to the set of equations, i.e., a zero pivot may not be encountered, but the solution vector may have wildly large components whose algebraic cancellation, when multiplying by the matrix , may give a very poor approximation to the right-hand vector . In such cases, with 0 will the solution vector obtained by zeroing the small ’s is often better. Replacing provide the closest that minimizes the distance to in the least square sense,

*

&

(A.92)

whith denoting the residual of the solution.

, / ! 1 "

The simplest example is of an ill-conditioned situation deals with the parametrization of a cubic curve degenerating into a straight line segment . The geometric data for the 46

−3

x 10

n1 = (0, 0, 1)

1

p1 = (0, 0, 0) p2 = (0, 1, 0)

0.8 0.6 0.4 0.2

n = (1, 0, 0)

Z0

2

n = (1, 0.0001, 0) 2

−0.2

n = (1, 0.001, 0) 2 n = (1, 0.0005, 0)

−0.4

2

n = (1, 0.01, 0) 2

−0.6

n2 = (1, 0.005, 0)

−0.8 −1 0 0.2 0.4

y 0.6

0 0.5 0.8

1

x

−3

x 10

1.5 1

2

n1 = (0, 0, 1)

−3

x 10

p1 = (0, 0, 0) p2 = (0, 1, 0)

n2 = (1, 0.0001, 0) n2 = (1, 0, 0)

1 0.8

n2 =0) (1, 0.001, 0) n2 = (1, 0.0005,

Z0.6 0.4

n2 = (1, 0.005, 0) n2 = (1, 0.01, 0)

0.2 0

0

−0.2

0.5

−0.4

x

−0.6

1

−0.8 1.5 −1 0

0.2

Figure 29:

y0.4

0.6

0.8

2 1

varies from (1, 0.01, 0) to (1, 0, 0)

47

−3

x 10

/%%'& 1 6 /0& 1

/%% 1 /%5& 1

limity case are: and . The vanishing constraints in (3.36) result in a singular matrix . Using SVD, we study the behavior of the matrix as . We have use the curve reconstruction routine, with data . The code delivers uniformly stable results converging to the limity case. For illustration, Figures A A show the results of the curve reconstruction for values varies from to , and from to . Note we use the different scales in Figures to illustrated the convergence property.

1 * / % % 5 & ' & 1 '& '& / '& '& 1 /% %& '& 1

/ %& '& 1 /% '& 1

48

/ %'&

B Subroutines The following subroutines have been implemented into the GMP package. input head.f

the interface with data extracted from MRI that stores all the connectivity informations.

filter.f

A subroutine used to remove extra points that are redundant in the input file.

HermCur1.f

A subroutine which generates rational curve function with CBDs

ratio.f

A subroutine that reconstructs the

HermCur2.f

A subroutine which generates Hermite curve function with linear normals

HermRec.f

A subroutine that reconstructs the surface by using bicubic Hermite rectangular patch interpolations.

HermCur3.f

A subroutine which generates Hermite curve function with cubic normals

biqua.f

A subroutine that reconstructs the surface by using biquartic Hermite rectangular patch interpolations.

surface by using IDW interpolation.

matr.f

A subroutine that solves a system of linear vector equations using Guassian elimination.

svdcmp.f

A subroutine that solves a system of linear vector equations using SVD.

head

The txt file that is used as a input for the interface

49

References [1] G.Xu, C.Bajaj. “Adaptive Model Reconstruction by Triangular A-patches.” TICAM Report, 2003 [2] Bajaj C., Chen J., Xu G. “Modeling with Cubic A-patches.” ACM Transactions on Graphics, vol. 14, 103–133, April 1995 [3] Bajaj C., Chen J., Xu G. “Modeling with Cubic A-patches ACM Transactions on Graphic.” 14, vol. 2, 103–133, April,(1995) [4] Xu G., Bajaj C., Huang H. “C1 Modeling with A-patches from Rational Trivariate Functions Computer Aided Geometric Design.” 18, vol. 3, 221–243, 2001 [5] B.Guo. “Modeling Arbitrary Smooth Objects with Algebraic Surfaces.” PhD thesis, vol. Computer Sciences Cornell University, 1991 [6] T.W.Sederberg. “Piecewise Algebraic Surface Patches.” Computer Aided Geometric Design, vol. 2, 53–59, 1985 [7] H.P.Seidel. “Symmetric Recursive Algorithms for Surfaces: B-patches and the De Boor Algorithm for Polynomials Over Triangles.” Constr. Approx., vol. 7, 257–279, 1991 [8] C.T.Loop, T.DeRose. “Generalized B-spline Surfaces of Arbitray Topology.” Computer Graphics, vol. 24, 347–356, 1990 [9] W.Boehm, G.Farin J. “A Survey of Curve and Surface Method in CAGD.” Computer Aided Geometrc Design, vol. 1, 1–60, 1984 [10] G.Farin. “Curves and Surfaces for Computer Aided Geometric Design.” Academic Press, 1990 [11] Peters J. “C1 Surface Splines.” SIAM J. of Numer. Anal., vol. 32-2, 645–666, 1995 [12] Peters J. “Biquartic C1-surface Splines Over Irregular Meshes.” CAD, Jan 95 [13] Peters J. “Surfaces of Arbitrary Topology Constructed from Biquadratics and Bicubics.” Designing fair curves and surfaces, vol. SIAM Publications, N. Sapidis, 1994 [14] Peters J. “Smooth Free-form Surfaces Over Irregular Meshes Generalizing Quadratic Splines.” CAGD, pp. 347–361, 1993

50

[15] A.R.M.Piah. “Construction of Smooth Surfaces by Piecewise Tensor Product Polynomials.” CS Report, vol. 91/04, University of Dundee, UK, 1991 [16] O.Reif. “Biquadratic G-spline Surfaces.” Preprint 93-4, vol. Math Ins A, University Stuttgart, 1993 [17] M.Sabin. “Non-rectangular Surface Patches Suitable for Inclusion in a B-spline Surface.” Proceedings of Eurographics, vol. North Holland, 57–69, 1983 [18] T.N.T.Goodman. “Closed Surfaces Defined From Biquadratic Splines.” Constructive Approximation, vol. 7, 149–160, 1991 [19] E.Catmull, J.Clark. “Resursively Generated B-spline Surfaces on Arbitary Topological Meshes.” CAD, vol. 10, 350–355, 1978 [20] D.Doo. “A Subdivision Algorithm for Smoothing Down Irregularly Shaped Polyhedrons.” Proceedings on interactive techniques in computer aied design., vol. Bologna, 157–165, 1978 [21] L.Shirman, C.Sequin. “Local Surface Interploation with Bezier Patches.” Computer Aided Geometric Design, vol. 4, 279–296, 1990 [22] R.F.Sarraga. “ Interplation of Generally Unrestricted Cubic Bezier Curves.” Computer Aided Geometric Design, vol. 1-2, 23–40, 1987 [23] Xue D., Demkowicz L. “An Interface Between Geometrical Modeling Package(GMP) and Mesh-Based Geometry(MBG).” ICES Report, vol. 03, 2003 [24] M.A.Sabin. “Transfinite Surface Interpolation.” Proceedings of a sixth conference on Mathematics of Surfaces, 1996 [25] Peter J. “Smooth Interpolation of a Mesh of Curves.” Constructive Approximation, vol. 7, pp.221–246, 1991 [26] Mortenson M.E. Geometric Modeling. Wiley Computer Publishing, 1997 [27] Hoschek J., Lasser D. “Fundamentals of Computer Aided Geometric Design.” A K Peters, 1993 [28] Zhang Y., Bajaj C., Sohn B.S. “Adaptive and Quality 3D Meshing from Imaging Data.” Proceedings of 8th ACM Symposium on Solid Modeling and Applications, vol. June 16-20, 286–291, 2003

51

[29] Zhang Y., Bajaj C., Sohn B.S. “3D Finite Element Meshing from Image Data.” The special issure of Computational Methods in Applied Medhanics and Engineering on unstructral mesh generation, 2003 [30] Xue D., Demkowicz L. “Geometrical Modeling Package Version 2.0.” TICAM Report, vol. 02-30, 2002 [31] Demkowicz L., Bajer A., Banas K. “Geometrical Modeling Package.” TICAM Report, vol. 02-06, 1992 [32] Demkowicz L. “2D hp-Adaptive Finite Element Package(2Dhp90) Version 2.0.” TICAM Report, vol. 02-06, 2002 [33] Demkowicz L., Pardo D., Rachowicz W. “3D hp-Adaptive Finite Element Package(2Dhp90).” TICAM Report, vol. 02-24, 2002 [34] Gordon W.J., Hall C.A. “Transfinite Element Methods: Blending Function Interpolation over Arbitary Curved Element Domains.” Numer. Math., vol. 21, 109–129, 1973 [35] Gordon W.J. “Blending Function Methods of Bivariate and Multivariate Interpolation and Approximation.” SIAM J. Numer. Anal, vol. 8, 158–177, 1971 [36] Zhang Y., Bajaj C., Sohn B.S. “Adaptive Multiresolution and Quality 3D Meshing from Imaging Data.” TICAM Report, vol. 02-42, 2002 [37] Shepard D. “A Two-dimensional Interpolation Function for Irregularly spaced Data.” 23rd ACM National Conference, pp. 517–524, 1968 [38] Gentle J.E. “Singular Value Factorization.” Numerical Linear Algebra for Applications in Statistics, vol. Berlin: Springer-Verlag, 102–103, 1998 [39] Golub G.H., Van Loan C.F. “The Singular Value Decomposition.” 2.5.3 and 2.5.6 in Matrix Computations, vol. 3rd ed. Baltimore,Johns Hopkins University Press, 70–71, 1996 [40] Nash J.C. “The Singular-Value Decomposition and Its Use to Solve Least-Squares Problems.” Compact Numerical Methods for Computers: Linear Algebra and Function Minimisation, vol. Bristol, England: Adam Hilger, 30–48, 1990 [41] Press W.H., Flannery B.P., Teukolsky S. A.and Vetterling W.T. “Singular Value Decomposition.” Numerical Recipes in FORTRAN: The Art of Scientific Computing, vol. Cambridge University Press, 51–63, 1992

52

[42] Press.W.M, W.M.Flannery, B.P.Teukolsky. Numerical Recipes: The Art of Scientific Computing. New York, NY, Cambridge University Press, 1986

53