A Parameterization Transporting Surface Offset Construction Method ...
A Parameterization Transporting Surface Offset Construction Method Based on the Eikonal Equation Yuanli Wang∗ Fran¸ccois Guibault† Ricardo Camarero‡ and Ko-Foa Tchon§ ´ Ecole Polytechnique de Montr´eal,Montr´eal, C.P. 6079, Succ. Centre-ville, Qu´ebec, H3C 3A7, Canada
EXTENDED ABSTRACT This paper presents a new method to construct offset surfaces and transport the parameterization of the original surface to the propagated surface. This method, based on the solution of the equation ∇φ · ∇φ = 1, derived from the Eikonal equation, comprises three essential steps : (1) construct an uniform Cartesian grid fully covering the domain of interest, consisting of discretized curves and surfaces; (2) solve the Eikonal equation using a modified fast sweeping numerical scheme on the Cartesian grid, and, (3) march in space from the original surface, starting from the discretized geometric front points, along the local normal direction, by iteratively solving a propagation equation, ~x t = ∇φ/|∇φ|, using a fourth-order Runge-Kutta scheme. The most attractive feature of this technique is that it can exactly map the parameterization information from the original front to the offset front. This technique can be used to offset curves and surfaces, and to generate meshes in boundary layer areas. Numerical experiments are presented to demonstrate the accuracy and efficiency of the method. Technical Topic or Category : CFD Student Paper Competition in Grid generation and geometry definition for complex configurations
Nomenclature φ Minimum Euclidean distance between each Cartesian grid node to geometric front ~x Geometric front position vector t Time step ~xt Time derivative d Given propagated distance Γ Original geometric front Γp Propagated geometric front Ω Interested domain Subscript (i, j, k) Index for one Cartesian grid node in 3D
I.
Introduction
Boundary Offset consists in propagating or advancing the boundaries of a physical domain either inward or outward along their normal direction. It is sometimes called parallel curve or surface construction, 1 and solution to this problem has applications in many fields, including CAD modeling, mold forming and NC machining, just to name a few. In the present paper, we focus on the automatic partitioning of a computational domain for analysis purposes. Extending earlier work by Guibault and al., 2 our goal is to ∗ PhD
Student, Department of Mechanical Engineering, [email protected] Department of Computer Engineering, [email protected], AIAA Member. ‡ Director, Department of Mechanical Engineering, [email protected], AIAA Member. § Research Associate, Department of Mechanical Engineering, [email protected], AIAA Member.
† Professor,
1 of 15 American Institute of Aeronautics and Astronautics
enhance offset surface construction methods that are used in support of automatic blocking procedures for hybrid mesh generation. Dealing with offset surface construction, two important problems arise, for which no perfectly satisfactory method has been proposed yet. The first problem relates to the potential self-intersection of the front as it grows from the original surface in concave regions, and the second, to the establishment of a common parameterization for the original and offset surfaces. Local warping may occur when the offset distance is greater than the local curvature radius in the concave regions. Global self-intersection, on the other hand, may occur, when the distance between two distinct points on the curve or surface reach a local minimum. This second problem, related to parameterization, is not cited as often as the first as a potential issue in offset construction. Numerous methods have however been developed to deal with incidental problems that could be cured if parameterization could be transported from the original surface to its offset. These problems include transport of trim curves3 and tool path generation using offset curves and surfaces,? and, in the present case, semi-structured mesh generation in near wall regions. The purpose of this work is to develop a propagation method which not only prevents local and global self-intersection during offsetting but also can entirely convey parameterization from the original surfaces to the propagated surfaces. A variety of contributions deal with the computation of offset curves and surfaces. They can be classified in two types: • direct offset methods (DOM), which propagate curves or surfaces directly based on original geometric data; • indirect offset methods (IOM), which cast curves or surfaces offset problem into partial differential equations (PDE), with the numerical solution to PDE being the offset curves or surfaces. A.
Direct Offset Methods (DOM)
The advantage of DOM is that the entire or partial original parameterization information can be well preserved, but the self-intersections problem can not be avoided, and extra care for removing self-intersections is necessary. The capacity to effectively eliminate these self-intersections is an important criterion to evaluate each DOM method. Many efforts have been made to eliminate self-intersections. One early work in DOM is the AdvancingLayers Method developed by Pirzadeh4 based on a grid-marching strategy. The layers stop advancing before self-intersection occurs on the marching front. Hoschek5 developed a curves offset method aimed at deleting undesirable portions of offset curves with cusps or with self-intersection points. Farouki 6 described algorithms which first decompose the original surfaces into parametric patches, and then use Hermite interpolation to construct offset surfaces. Blomgren,7 Tiller and Hanson,8 and Coquillar9 described various ways to offset control polygon for NURBS curves, where offset curves were generated based on offset of the control polygon. Nachman10 extended this idea to propagate surfaces. Guibault2 developed an algorithm to offset NURBS BRep represented curves or surfaces, where self-intersections can be detected and removed from the offsetting curves or surfaces. Piegl and Tiller11 sampled offset curves and surfaces based on bounds on the second derivatives to avoid self-intersections. Kumar et al.3 presented an approach for offsetting a trimmed NURBS surface based on the assumption that self-intersections shall not occur. In the method developed by Sun and al.,12 in areas where local self-intersections may occur, control points are repositioned to reduce local curvature, and the rest of the control points remain unchanged. B.
Indirect Offset Methods (IOM)
Most IOM methods can effectively avoid self-intersections problems, but at the cost of partially or completely losing the parameterization information stored in the original geometric front after propagation. In general, to restore a similar connectivity between orignal and offset parameterization is not a trivial task using IOM methods. In IOM, there are mainly three types of PDEs used to propagate curves or surfaces. The earliest one is the Marker Particle Method .13–15 The major contribution of this work is that it is the first effort to explicitly model front propagation problems using PDEs.
2 of 15 American Institute of Aeronautics and Astronautics
~xt = F ~n(~xs )
(1)
In this equation, ~x is the geometric front position vector, ~xt is its time derivative, ~xs its space derivative along parameter s, F is the offset speed, which is a given function, and ~n is the normal vector. Replacing ~xt and ~xs with their finite difference approximations, the original equation is then transferred into a set of algebraic equations. The offset curves or surfaces are generated by successively solving this transformed equation. Since this method captures a strong solution to Eqn. (1), self-intersections are unavoidable, and an extra scheme to detect and remove the warped areas from the result is required. Another effort in IOM is the Level set method developed by Sethian and Osher,15 which models front propagation problems as a hyperbolic differential equation. The significant improvement of this work is that it intrinsically prevents self-intersections by constructing a weak solution to the offset problem. Numerical scheme and optimization method for solving Eqn. (2) are discussed in works by Sethian. 16, 17 φt + ∇φ · ~xt = 0
(2)
∇φ · ∇φ = 1
(3)
Kimmel18 solved the level set equation to offset NURBS curves and surfaces in 1993. The weaknesses of this method include: (1) The fact that result of offsetting curves or surfaces are dependent on many other factors, such as offset speed F , curvature k, re-initialization of the level sets during calculation to maintain a nice level set representation; and (2) The fact that two-dimensional parameterization of the original surface is completely lost. Since the offsetting surfaces are represented by triangulated surfaces, a costly task is to transform these offsetting surfaces to parameterized surfaces. Recently, we have proposed a third type of PDE for surface offset construction. ? The Offset Distance Equation, derived from the Eikonal equation is stated as :
Solution of this equation represents the shortest distance from any space position to the boundary. It is thus a simplification of the level set mathematical model. The Fast Sweeping Scheme (FSS) is used to solve this equation (please see? for details), as described in the works of Tsai19 and Zhao.20 FSS improves the computational complexity from Ω(log(N )) in the level set method to Ω(N ), where N is the total number of Cartesian grid points. In our previous work, as in the level set method, offset surfaces where extracted from the computed φ solution as isovalue surfaces, resulting in triangulated surface representation, which had to be reparameterized. Even though many valuable works have been done in parameterizing triangulated surfaces, such as the works by Sheffer,? and Hormann,? these methods supply a way to re-construct two-dimensional parameterization based on triangulated surfaces but cannot be used to establish a parametric relationship between the original and offset surface. Overall, among the reviewed DOM and IOM methods, DOM can entirely or partially transfer parameterization describing the original surface, but can not supply an elegant way to eliminate self-intersections from the resulting offset. In IOM, on the other hand, the level set equation and the Offset Distance Equation can efficiently prevent self-intersections, although entirely loosing the original parameterization information. The present work proposes to combine the advantages of the two types of methods to build a new offset construction method that maintains parametric connectivity between the original and offset surfaces, and still avoid self-intersections through the use of a weak solution to the shortest distance problem. The remainder of the paper is structured as follows: in the next section, the computational method is presented. Preliminary results are then discussed, both in terms of accuracy and speed, and conclusions are sketched out. Each of the following sections will be thoroughly reviewed and significantly enhanced in the final version of the paper.
II.
Method
The propagation process proposed here is directly inspired from DOM. To enforce parametric consistency between Γ, the original surface, and Γp , its offset, we start the propagation from Γ. All the points located on Γ are propagated along their local normal directions. The main difference between the current method and traditionnal DOM methods lies in the way local normal directions are computed. In DOM, the local normal 3 of 15 American Institute of Aeronautics and Astronautics
directions are computed using geometric information, such as the positions of points located on Γ. In the method proposed, we instead use the φ field to compute local normal directions. The function φ represents the weak solution to Eqn. (3). This property guarantees that local normal directions will not intercross each other during propagation when a proper time step is used. The central idea of the method is thus to propagate each discretized point of the original geometric front Γ along it’s normal direction in φ space rather than geometric space. The propagation equation is given as follows : ~xn+1 = ~xn + F · dt (4) In Eqn. (4), F is the unit normal propagation velocity, F = ∇φ/|∇φ|. For each new propagated point ~xn+1 , there are two stopping criteria; when one of them is satisfied, the propagation is stopped : (1) If φ is equal to or greater than the specified propagation distance; (2) If φn+1 − φn < 0, or in other words, if φ starts to decrease. After propagation, we sequentially connect all the final ~x n+1 to construct Γp . All the points on Γ and Γp share a one to one parametric mapping. There are three essential steps in this method : 1. find a domain Ω, and discretize Ω into an uniform Cartesian grid; 2. calculate φ for each Cartesian grid node using fast sweeping algorithm;? 3. offset each node locating on Γ along its local normal direction according to the given distance using fourth-order Runge-Kutta method.? A.
Discretize Ω
The purpose of computing a Cartesian frame is to assign a scalar value φ at each Cartesian grid node. φ is the minimum Euclidean distance from each Cartesian grid node to Γ. This φ field can be used : • to compute the local unit normal propagation velocity ∇φ/|∇φ| for each node on Γ; • to compute φ values for each propagated point using the φ values at Cartesian grid points using interpolation. The φ function can thus directly be used as a stopping criterion. If φ ≥ d, stop the propagation; otherwise, continue propagation, until the criterion is reached; • as an indirect stopping criterion in the propagation process. If the surfaces to be propagated are too close, the user may give a large propagation distance, and collisions between two offset surfaces may occur. In this case, we can use φ to prevent collisions. When φn+1 − φn < 0, we stop propagation Fig. 2. Suppose we are given a geometry Γ (Fig. 1(a)). We first compute an interested domain Ω which should be large enough to cover Γ, and then discretize Ω into an uniform Cartesian grid (Fig. 1(b)). B. 1.
Solve φ Equation for representing φ
We are now faced with the problem of representing the relationship between the offset distance φ and the space coordinates. ~0 be a point on Γ; X ~ be a Cartesian grid point in Ω, Γ ⊂ Ω. Suppose that the distance between Let X ~ ~ ~ this distance is represented by φ. Then the X0 and X is the smallest Euclidean distance between Γ and X, distance map function is: ~ −X ~0 | φ = |X (5) rewrite eqn. (5), we get φ=
p
(x − x0 )2 + (y − y0 )2 + (z − z0 )2
(6)
so the offset distance equation can be written as follows (
dφ dφ dφ 2 ) + ( )2 + ( )2 = 1 dx dy dz
4 of 15 American Institute of Aeronautics and Astronautics
(7)
(a) Original Geometry
(b) Uniform Cartesian Grid
(c) Initializing exact φ for all near boundary nodes
Figure 1. Discretizing Ω and initializing φ
in which, φ is 0 when (x, y, z) lies on Γ. By now, we have deducted a general mathematical form to model the offset problem: φ2x + φ2y + φ2z = 1 φ=0 if (x, y, z) ∈ Γ
(8)
which can be restated more generally as equation 3. 2.
Numerical Scheme for solving φ
The unique solution to the Eqn. (8) is √ 2 ( φt+1 i,j,k =
(a+b+c)+
3dh +2(ab+bc+ac−a2 −b2 −c2 ) 3
min(a, b, c) + dh
if ab + bc + ac − a2 − b2 − c2 ≥ 0 if ab + bc + ac − a2 − b2 − c2 < 0
t+1 t φnew i,j,k = min(φi,j,k , φi,j,k )
(9)
t The solution process consists in always using φnew i,j,k to replace the old φi,j,k values in the sweeping calculation process. A more detailed description of the fast sweeping algorithm can be found in. ?, 20 Before solving Eqn. 8, we first assign exact values at Cartesian grid points in or near Γ (Fig. 1(c)) to enforce the boundary condition. We then apply the fast sweeping algorithm to calculate the φ field for the remainder of the domain.
C.
Propagate parameters
The propagation equation is ~xt =
∇φ |∇φ|
(10)
replacing ~xt by a finite difference scheme, Eqn. (10) can be rewritten as ~xn+1 = ~xn ± The term ∇φ |∇φ|
∇φ |∇φ|
∇φ ∗ dt |∇φ|
(11)
can be viewed as a unit propagation speed in the normal direction. For an outward propaga-
tion, is positive; for an inward propagation, into three steps:
∇φ |∇φ|
is negative. The propagation process can be divided
5 of 15 American Institute of Aeronautics and Astronautics
cbr1
cbr2
Figure 2. Stop propagation when φn+1 − φn < 0
1. calculate φx , φy , and φz We first calculate φx , φy , and φz at each Cartesian grid node using a two-sided finite difference scheme for internal nodes, and one-sided finite difference scheme for boundary nodes. Then, we calculate φ x , φy , and φz for all the points to be propagated at Γ using interpolation (Fig. 3(a)). 2. propagate points on Γ We use a forth-order Runge-Kutta method? to obtain the solution to Eqn. (11). We stop propagation once we detect that the value of φ at position ~xn+1 is greater or equal to the given propagation distance. If the value of φ at the final position ~xn+1 is greater than the given propagation distance, we use linear interpolation to get the exact position ~xf inal (Fig. 3(b)). 3. construct Γp Finally, we construct Γp by sequentially connecting all points on Γp using the parameterization information conveyed from Γ.
III.
Results
We verify the effectiveness of this scheme from two aspects: • Accuracy, in which we verify that the error for the propagated front has first order accuracy. • Efficiency, in which we verity that the complexity of this scheme is Ω(N ). We use two test cases to verity this algorithm, a circle, which is centered at (0.0), radius is 0.5; and a cube, which is centered at (0.0.0), side length is 0.5. We propagate these test cases to their inward direction. The offset distance is 0.1. The time increment is 0.01. We fix the number of discretizing points at front, and vary the Cartesian grid size from 202 to 3202 for the circle, from 203 to 1603 for the cube. For the circle, the front is discretized into 40 points; for the cube, the front is discretized into 82308 triangles. These tests were carried out on a standard desktop PC running Linux, and equiped with an AMD Athlon processor running at 1400 MHz, and a memory size of 512Mb. A.
Accuracy Validation
To measure accuracy and order of convergence of the method, we use the L 2 and L∞ norms21 to measure the error on the position of the propagated front. As shown in the results presented, our method converges 6 of 15 American Institute of Aeronautics and Astronautics
(φ x1, φ y1)
(φ x2, φ y2)
(x0, y 0)
(φ x , φ y )
Γ (φ x4, φ y4)
(φ x3, φ y3)
(x1, y1) (xfinal, y final) (x2 , y 2)
(a) We first calculate (φxi , φyi ) at each Cartesian grid node using a finite difference scheme, and then interpolate (φx , φy ), which are the points located on Γ, using (φx1 , φy1 ), (φx2 , φy2 ), (φx3 , φy3 ) and (φx4 , φy4 )
(b) If φ ≥ d at the last propagation point (x2 , y2 ), then calculate the final position of (xf inal , yf inal ) using a linearly interpoled φ value between (x1 , y1 ) and (x2 , y2 ). So φ = d at (xf inal , yf inal ).
Figure 3. Propagation process
for increasing mesh size with first order accuracy. Table 1 illustrates the comparisons of L 2 and L∞ for inwardly propagating a circle; Table 2 illustrates the comparisons of L2 and L∞ for inwardly propagating a cube. From Fig. 4 and Fig. 6, we can see that both L2 and L∞ converge when we increase the mesh size. Fig. 5(a) illustrates the result of propagating a circle in its inward direction using a coarse Cartesian mesh (mesh size=202); Fig. 5(b) illustrates the same propagation of a circle using a refined Cartesian mesh (mesh size=3202). Fig. 7(a) shows the result of propagating a cube in its inward direction using a coarse Cartesian mesh (mesh size=203); and Fig. 7(b) is the result for propagating the same cube in its inward direction using a refined Cartesian mesh (mesh size=1603). From these figures, we can see that more precise offset surface definitions are obtained when we increase the Cartesian mesh density. From Fig. 5 and Fig. 7, we can intuitively observe that mesh size severely influences the propagated front. The output circle in Fig. 5(b) looks smoother than the circle in Fig. 5(a), and the behavior for sharp corners in the output cube in Fig. 7(b) looks much better than the cube in Fig. 7(a). As shown by these experiments, the propagated front position directly depends on mesh size, since the ∇φ unit propagation velocity in the normal direction, |∇φ| , is entirely dependent on the precision of the φ field, which is the numerical solution to Eqn.(8). The fast sweeping scheme used to solve Eqn.(8) is first order accurate with respect to mesh size, this dominates the propagation accuracy. The convergence of the φ field has been discussed in more details in the previous work.? 1.
Circle Mesh Size L2 L∞
202 0.00497535 0.00548393
402 0.00315108 0.00367538
802 0.00244957 0.00319467
1602 0.00182176 0.00246195
Table 1. Norms comparison (inwardly propagate a circle)
7 of 15 American Institute of Aeronautics and Astronautics
3202 0.00136155 0.00172537
Nomrs comparision for circle(inward propagation) 0.01
log(norm)
norm2 of propagated circle norminfty of propagated circle
0.001 0.001
0.01
0.1
1
log(dh)
Figure 4. Norms comparison (inwardly propagate a circle)
2.
Cube Mesh Size
203
403
803
1603
L2 L∞
0.0476718 0.17472
0.0393768 0.128057
0.0368557 0.100638
0.036115 0.0862622
Table 2. Norm Comparison for the propagated surface (Cube Inward Propagation)
B.
Efficiency Validation
In order to evaluate the efficiency of the scheme, we compare the timing results obtained for the calculations presented in our previous work? with those performed specifically using the new propagation method. The total calculation time includes three parts: 1. initializing time, which is the time used to initialize the exact distance for the points near or on front; 2. fast sweeping time, which is the time used to solve numerical solution φ to equation Eqn. 8; 3. iso-front capturing time or propagation time iso-front capturing time is the time used to capture iso-front, which is called propagated front in this paper. This algorithm was used in our previous work.? propagation time is the time used to propagate all the points located on the original front to their destination positions. This algorithm is used in the present work. Thus the total time consumed is calculated either by the formula of total time = initializing time + f ast sweeping time + iso − f ront capturing time 8 of 15 American Institute of Aeronautics and Astronautics
(12)
cbr
cbr
(a) mesh size = 202
(b) mesh size = 3202
Figure 5. Propagate a circle (inward direction)
in our old algorithm,? or calculated by the formula of total time = initializing time + f ast sweeping time + propagation time
(13)
in this work. To compare the efficiency, we will measure these times consumed for each part separately. Table 3 illustrates the time consumed for each part to propagate a circle in its inward direction; Table 4 illustrates the time consumed for each part to propagate a cube to its inward direction. Fig. 8 displays the time slope consumed in each part for inwardly propagating a circle, and Fig. 9 displays the time slope consumed in each part for inwardly propagating a cube. From Fig. 8 and Fig. 9, we can see that initializing time, fast sweeping time and iso-surface capturing time are linearly dependent on the mesh size as discussed in our previous paper;? propagation time is dependent on the number of points for front initialization and time increment. Once these two factors are fixed, propagation time is a constant value. The reason for this is that in the propagation process, we start our propagation from the original front. Once the number of initial points in the original front is fixed, the next factor to influence the propagation time is time increment dt. If dt is too small, the points on Γ have to be propagated many times to reach their destinations. Thus the propagation time is not be influenced by mesh size. It is only a function to the number of original front points and time increment. 1.
Circle Mesh Size
202
402
802
1602
3202
Initializing Time Fast Sweeping Time iso-surface capturing time Propagation Time
0.00457 0.003768 0.000212 0.016103
0.007592 0.023966 0.000712 0.016017
0.015335 0.160162 0.001863 0.015958
0.040478 1.2056 0.005838 0.016986
0.131439 9.11685 0.022492 0.017038
Table 3. Time Comparison for Circle (Inward Propagation)
2.
Cube
9 of 15 American Institute of Aeronautics and Astronautics
Nomrs comparision for cube(inward propagation)
log(norm)
1
norm2 of propagated surface norminfty of propagated surface
0.1
0.01 0.01
0.1
1
log(dh)
Figure 6. Norms comparison (inwardly propagate a cube)
Mesh Size
203
403
803
1603
initializing time Fast Sweeping time Propagation Time 1 Propagation Time 2
1.15137 0.044765 0.01799 0.088228
2.87791 0.54249 0.418908 0.082101
13.9404 9.05241 8.01082 0.082883
97.7702 162.96 147.378 0.088191
Table 4. Time Comparison for a Cube (Inward Propagation)
10 of 15 American Institute of Aeronautics and Astronautics
(a) mesh size 203
(b) mesh size 1603
Figure 7. propagate a cube (inward direction)
Times comparision for circle (inward propagation) 1000
initializing time (s) fast sweeping time (s) iso-surface capturing time (s) propagation time (s)
100
log(time)
10
1
0.1
0.01
0.001
0.0001 0.001
0.01
0.1 log(dh)
Figure 8. Propagate a circle (inward direction)
11 of 15 American Institute of Aeronautics and Astronautics
1
Times comparision for cube(inward propagation) 1000
initializing time (s) fast sweeping time (s) iso-surface capturing time (s) propagation time (s)
100
log(time)
10
1
0.1
0.01 0.01
0.1 log(dh)
Figure 9. Propagate a cube (inward direction)
12 of 15 American Institute of Aeronautics and Astronautics
1
IV.
Method characteristics and Applications
The major advantages of this method is that the propagation lines will not intercross each other when proper time step t is used (as shown in Fig. 10); two-dimensional parameterization stored in Γ are entirely carried out to the propagated front.
Figure 10. Propagation lines at sharp conner
This method can be applied to many types of surfaces, i.e. parameterized surfaces, tessellated surfaces, trimmed surfaces. It also supplies a better flexibility to offset partial surfaces for a geometry. It has a strong ability to control the offset distance to obtain uneven distanced surfaces. For exemple, we can establish a certain connection between local parameters (ui , vi ) of the original surface and propagation distance to control the offset process. This method can be used in many applications, for example, to shrink or enlarge a geometry with any complexity (Fig. 11(a)), to generate a structured grid in the area of boundary layers for viscous flow (Fig. 11(b)). We present results from some applications at hand to illustrate the use of the method on concrete application. Fig. 12(a) is a mesh generation example in a thin layer around the body of a diffuser; Fig. 12(b) is an example to offset the body of a diffuser in the inward direction to generate a thin layer around the body of diffuser.
(a) Offsetting front
(b) Generating mesh for boundary layers
Figure 11. Examples of application
13 of 15 American Institute of Aeronautics and Astronautics
(a) Mesh generation result
(b) Surface propagation result
14 of 15 American Institute of Aeronautics and Astronautics
V.
Conclusion
The main features of the proposed method are : • The parameterized information stored in the original geometry front can be entirely transformed to the offset curves and surfaces. • The local normal directions at front points will not inter-cross each other. • This technique can be applied to many types of front, i.e. a tesselated front, parameterized front, trimmed surfaces, etc. • The accuracy of this method is first order to mesh size. • The complexity of this method is Ω(N ), N is the total number of Cartesian grid nodes. The method can be extended to: (1) generate structured grids near boundary area for studying viscous flow. It guarantees that the mesh is orthogonal to its geometric front, since the propagating paths are developed along local geometric normal directions; (2) generate geometric information, for example, to shrink or enlarge a geometry; and (3) generate an uneven distanced skin by properly setting the offset distance function.
References 1 Maekawa, 2 Guibault,
T., “An Overview of Offset Curves and Surfaces,” Computer-Aided Design, Vol. 31, No. 3, 1999, pp. 165–173. F., Labb, P., Garon, A., and Camarero, R., “Automatic Topologically Based Blocking for Hybrid Grid Gener-
ation,” . 3 Kumar, R., Shastry, K., and Prakash, B., “Computing offsets of trimmed NURBS surfaces,” Computer-Aided Design, Vol. 35, No. 5, 2003, pp. 411–420. 4 Pirzadeh, S., “Unstructured Viscous Grid Generation by Advancing-Layers Method,” No. AIAA-93-3453, Monterey, CA, August. 5 Hoschek, J., “Offset Curves in the Plane,” Computer-Aided Design, Vol. 17, No. 2, 1989, pp. 77–82. 6 Farouki, R., “The approximation of non-degenerate offset surfaces,” Computer Aided Geometric Design, Vol. 3, No. 1, 1986, pp. 15–43. 7 Blomgren, R., “B-Spline Curves,” 1981, Boeing Document, Class Notes, B-7150-BB-WP-281/D-441.2. 8 Tiller, W. and Hanson, E., “Offsets of Two-Dimensional Profiles,” IEEE Computer Graphics and Applications, Vol. 4, No. 9, 1984, pp. 61–69. 9 Coquillart, S., “Computing offsets of B-spline curves,” Computer-Aided Design, Vol. 19, No. 6, 1987, pp. 305–309. 10 Nachman, L., “Approximating Offset Surfaces in a NURBS Environment,” Second SIAM Conference on Geometric Design, Tempe, Arizona, November 1991. 11 Piegl, L. A. and Tiller, W., “Computing offsets of NURBS curves and surfaces,” Computer-Aided Design, Vol. 31, No. 2, February 1999, pp. 147–156. 12 Sun, Y., Nee, A., and Lee, K., “Modifying free-formed NURBS curves and surfaces for offsetting without local selfintersection,” Computer-Aided Design,In Press, 2004. 13 Sethian, J., An Analysis of Flame Propagation, Ph.d. dissertation, University of California, Berkeley, CA, 1982. 14 Sethian, J., “Curvature and the Evolution of Fronts,” Comm. in Math. Phys., Vol. 101, 1985, pp. 487–499. 15 Osher, S. and Sethian, J., “Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations,” Journal of Computational Physics, Vol. 79, 1988, pp. 12–49. 16 Sethian, J., “A Fast Marching Level Set Method for Monotonically Advancing Fronts,” Proc. Nat. Acad. Sci, Vol. 93, No. 4, 1996, pp. 1591–1595. 17 Sethian, J., Level set methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, Cambridge University Press, 1999. 18 Kimmel, R. and Bruckstein, A., “Shape offsets via level sets,” Computer-Aided Design, Vol. 25, No. 3, 1993, pp. 154–162. 19 Tsai, Y., Cheng, L., Osher, S., and Zhao, H., “Fast Sweeping Algorithms for a Class of Hamilton-Jacobi Equations,” SIAM Journal on Numerical Analysis, Vol. 41, No. 2, 2003, pp. 673–694. 20 Zhao, H., “A fast sweeping method for Eikonal equations,” to appear in Math. Comp. (2004), May 2004. 21 Fortin, A. and Garon, A., Les ´ ´ el´ ements finis: de la th´ eorie ` a la pratique, Ecole Polytechnique de Montr´eal, 2000.
15 of 15 American Institute of Aeronautics and Astronautics
EXTENDED ABSTRACT This paper presents a new method to construct offset surfaces and transport the parameterization of the original surface to the propagated surface. This method, based on the solution of the equation ∇φ · ∇φ = 1, derived from the Eikonal equation, comprises three essential steps : (1) construct an uniform Cartesian grid fully covering the domain of interest, consisting of discretized curves and surfaces; (2) solve the Eikonal equation using a modified fast sweeping numerical scheme on the Cartesian grid, and, (3) march in space from the original surface, starting from the discretized geometric front points, along the local normal direction, by iteratively solving a propagation equation, ~x t = ∇φ/|∇φ|, using a fourth-order Runge-Kutta scheme. The most attractive feature of this technique is that it can exactly map the parameterization information from the original front to the offset front. This technique can be used to offset curves and surfaces, and to generate meshes in boundary layer areas. Numerical experiments are presented to demonstrate the accuracy and efficiency of the method. Technical Topic or Category : CFD Student Paper Competition in Grid generation and geometry definition for complex configurations
Nomenclature φ Minimum Euclidean distance between each Cartesian grid node to geometric front ~x Geometric front position vector t Time step ~xt Time derivative d Given propagated distance Γ Original geometric front Γp Propagated geometric front Ω Interested domain Subscript (i, j, k) Index for one Cartesian grid node in 3D
I.
Introduction
Boundary Offset consists in propagating or advancing the boundaries of a physical domain either inward or outward along their normal direction. It is sometimes called parallel curve or surface construction, 1 and solution to this problem has applications in many fields, including CAD modeling, mold forming and NC machining, just to name a few. In the present paper, we focus on the automatic partitioning of a computational domain for analysis purposes. Extending earlier work by Guibault and al., 2 our goal is to ∗ PhD
Student, Department of Mechanical Engineering, [email protected] Department of Computer Engineering, [email protected], AIAA Member. ‡ Director, Department of Mechanical Engineering, [email protected], AIAA Member. § Research Associate, Department of Mechanical Engineering, [email protected], AIAA Member.
† Professor,
1 of 15 American Institute of Aeronautics and Astronautics
enhance offset surface construction methods that are used in support of automatic blocking procedures for hybrid mesh generation. Dealing with offset surface construction, two important problems arise, for which no perfectly satisfactory method has been proposed yet. The first problem relates to the potential self-intersection of the front as it grows from the original surface in concave regions, and the second, to the establishment of a common parameterization for the original and offset surfaces. Local warping may occur when the offset distance is greater than the local curvature radius in the concave regions. Global self-intersection, on the other hand, may occur, when the distance between two distinct points on the curve or surface reach a local minimum. This second problem, related to parameterization, is not cited as often as the first as a potential issue in offset construction. Numerous methods have however been developed to deal with incidental problems that could be cured if parameterization could be transported from the original surface to its offset. These problems include transport of trim curves3 and tool path generation using offset curves and surfaces,? and, in the present case, semi-structured mesh generation in near wall regions. The purpose of this work is to develop a propagation method which not only prevents local and global self-intersection during offsetting but also can entirely convey parameterization from the original surfaces to the propagated surfaces. A variety of contributions deal with the computation of offset curves and surfaces. They can be classified in two types: • direct offset methods (DOM), which propagate curves or surfaces directly based on original geometric data; • indirect offset methods (IOM), which cast curves or surfaces offset problem into partial differential equations (PDE), with the numerical solution to PDE being the offset curves or surfaces. A.
Direct Offset Methods (DOM)
The advantage of DOM is that the entire or partial original parameterization information can be well preserved, but the self-intersections problem can not be avoided, and extra care for removing self-intersections is necessary. The capacity to effectively eliminate these self-intersections is an important criterion to evaluate each DOM method. Many efforts have been made to eliminate self-intersections. One early work in DOM is the AdvancingLayers Method developed by Pirzadeh4 based on a grid-marching strategy. The layers stop advancing before self-intersection occurs on the marching front. Hoschek5 developed a curves offset method aimed at deleting undesirable portions of offset curves with cusps or with self-intersection points. Farouki 6 described algorithms which first decompose the original surfaces into parametric patches, and then use Hermite interpolation to construct offset surfaces. Blomgren,7 Tiller and Hanson,8 and Coquillar9 described various ways to offset control polygon for NURBS curves, where offset curves were generated based on offset of the control polygon. Nachman10 extended this idea to propagate surfaces. Guibault2 developed an algorithm to offset NURBS BRep represented curves or surfaces, where self-intersections can be detected and removed from the offsetting curves or surfaces. Piegl and Tiller11 sampled offset curves and surfaces based on bounds on the second derivatives to avoid self-intersections. Kumar et al.3 presented an approach for offsetting a trimmed NURBS surface based on the assumption that self-intersections shall not occur. In the method developed by Sun and al.,12 in areas where local self-intersections may occur, control points are repositioned to reduce local curvature, and the rest of the control points remain unchanged. B.
Indirect Offset Methods (IOM)
Most IOM methods can effectively avoid self-intersections problems, but at the cost of partially or completely losing the parameterization information stored in the original geometric front after propagation. In general, to restore a similar connectivity between orignal and offset parameterization is not a trivial task using IOM methods. In IOM, there are mainly three types of PDEs used to propagate curves or surfaces. The earliest one is the Marker Particle Method .13–15 The major contribution of this work is that it is the first effort to explicitly model front propagation problems using PDEs.
2 of 15 American Institute of Aeronautics and Astronautics
~xt = F ~n(~xs )
(1)
In this equation, ~x is the geometric front position vector, ~xt is its time derivative, ~xs its space derivative along parameter s, F is the offset speed, which is a given function, and ~n is the normal vector. Replacing ~xt and ~xs with their finite difference approximations, the original equation is then transferred into a set of algebraic equations. The offset curves or surfaces are generated by successively solving this transformed equation. Since this method captures a strong solution to Eqn. (1), self-intersections are unavoidable, and an extra scheme to detect and remove the warped areas from the result is required. Another effort in IOM is the Level set method developed by Sethian and Osher,15 which models front propagation problems as a hyperbolic differential equation. The significant improvement of this work is that it intrinsically prevents self-intersections by constructing a weak solution to the offset problem. Numerical scheme and optimization method for solving Eqn. (2) are discussed in works by Sethian. 16, 17 φt + ∇φ · ~xt = 0
(2)
∇φ · ∇φ = 1
(3)
Kimmel18 solved the level set equation to offset NURBS curves and surfaces in 1993. The weaknesses of this method include: (1) The fact that result of offsetting curves or surfaces are dependent on many other factors, such as offset speed F , curvature k, re-initialization of the level sets during calculation to maintain a nice level set representation; and (2) The fact that two-dimensional parameterization of the original surface is completely lost. Since the offsetting surfaces are represented by triangulated surfaces, a costly task is to transform these offsetting surfaces to parameterized surfaces. Recently, we have proposed a third type of PDE for surface offset construction. ? The Offset Distance Equation, derived from the Eikonal equation is stated as :
Solution of this equation represents the shortest distance from any space position to the boundary. It is thus a simplification of the level set mathematical model. The Fast Sweeping Scheme (FSS) is used to solve this equation (please see? for details), as described in the works of Tsai19 and Zhao.20 FSS improves the computational complexity from Ω(log(N )) in the level set method to Ω(N ), where N is the total number of Cartesian grid points. In our previous work, as in the level set method, offset surfaces where extracted from the computed φ solution as isovalue surfaces, resulting in triangulated surface representation, which had to be reparameterized. Even though many valuable works have been done in parameterizing triangulated surfaces, such as the works by Sheffer,? and Hormann,? these methods supply a way to re-construct two-dimensional parameterization based on triangulated surfaces but cannot be used to establish a parametric relationship between the original and offset surface. Overall, among the reviewed DOM and IOM methods, DOM can entirely or partially transfer parameterization describing the original surface, but can not supply an elegant way to eliminate self-intersections from the resulting offset. In IOM, on the other hand, the level set equation and the Offset Distance Equation can efficiently prevent self-intersections, although entirely loosing the original parameterization information. The present work proposes to combine the advantages of the two types of methods to build a new offset construction method that maintains parametric connectivity between the original and offset surfaces, and still avoid self-intersections through the use of a weak solution to the shortest distance problem. The remainder of the paper is structured as follows: in the next section, the computational method is presented. Preliminary results are then discussed, both in terms of accuracy and speed, and conclusions are sketched out. Each of the following sections will be thoroughly reviewed and significantly enhanced in the final version of the paper.
II.
Method
The propagation process proposed here is directly inspired from DOM. To enforce parametric consistency between Γ, the original surface, and Γp , its offset, we start the propagation from Γ. All the points located on Γ are propagated along their local normal directions. The main difference between the current method and traditionnal DOM methods lies in the way local normal directions are computed. In DOM, the local normal 3 of 15 American Institute of Aeronautics and Astronautics
directions are computed using geometric information, such as the positions of points located on Γ. In the method proposed, we instead use the φ field to compute local normal directions. The function φ represents the weak solution to Eqn. (3). This property guarantees that local normal directions will not intercross each other during propagation when a proper time step is used. The central idea of the method is thus to propagate each discretized point of the original geometric front Γ along it’s normal direction in φ space rather than geometric space. The propagation equation is given as follows : ~xn+1 = ~xn + F · dt (4) In Eqn. (4), F is the unit normal propagation velocity, F = ∇φ/|∇φ|. For each new propagated point ~xn+1 , there are two stopping criteria; when one of them is satisfied, the propagation is stopped : (1) If φ is equal to or greater than the specified propagation distance; (2) If φn+1 − φn < 0, or in other words, if φ starts to decrease. After propagation, we sequentially connect all the final ~x n+1 to construct Γp . All the points on Γ and Γp share a one to one parametric mapping. There are three essential steps in this method : 1. find a domain Ω, and discretize Ω into an uniform Cartesian grid; 2. calculate φ for each Cartesian grid node using fast sweeping algorithm;? 3. offset each node locating on Γ along its local normal direction according to the given distance using fourth-order Runge-Kutta method.? A.
Discretize Ω
The purpose of computing a Cartesian frame is to assign a scalar value φ at each Cartesian grid node. φ is the minimum Euclidean distance from each Cartesian grid node to Γ. This φ field can be used : • to compute the local unit normal propagation velocity ∇φ/|∇φ| for each node on Γ; • to compute φ values for each propagated point using the φ values at Cartesian grid points using interpolation. The φ function can thus directly be used as a stopping criterion. If φ ≥ d, stop the propagation; otherwise, continue propagation, until the criterion is reached; • as an indirect stopping criterion in the propagation process. If the surfaces to be propagated are too close, the user may give a large propagation distance, and collisions between two offset surfaces may occur. In this case, we can use φ to prevent collisions. When φn+1 − φn < 0, we stop propagation Fig. 2. Suppose we are given a geometry Γ (Fig. 1(a)). We first compute an interested domain Ω which should be large enough to cover Γ, and then discretize Ω into an uniform Cartesian grid (Fig. 1(b)). B. 1.
Solve φ Equation for representing φ
We are now faced with the problem of representing the relationship between the offset distance φ and the space coordinates. ~0 be a point on Γ; X ~ be a Cartesian grid point in Ω, Γ ⊂ Ω. Suppose that the distance between Let X ~ ~ ~ this distance is represented by φ. Then the X0 and X is the smallest Euclidean distance between Γ and X, distance map function is: ~ −X ~0 | φ = |X (5) rewrite eqn. (5), we get φ=
p
(x − x0 )2 + (y − y0 )2 + (z − z0 )2
(6)
so the offset distance equation can be written as follows (
dφ dφ dφ 2 ) + ( )2 + ( )2 = 1 dx dy dz
4 of 15 American Institute of Aeronautics and Astronautics
(7)
(a) Original Geometry
(b) Uniform Cartesian Grid
(c) Initializing exact φ for all near boundary nodes
Figure 1. Discretizing Ω and initializing φ
in which, φ is 0 when (x, y, z) lies on Γ. By now, we have deducted a general mathematical form to model the offset problem: φ2x + φ2y + φ2z = 1 φ=0 if (x, y, z) ∈ Γ
(8)
which can be restated more generally as equation 3. 2.
Numerical Scheme for solving φ
The unique solution to the Eqn. (8) is √ 2 ( φt+1 i,j,k =
(a+b+c)+
3dh +2(ab+bc+ac−a2 −b2 −c2 ) 3
min(a, b, c) + dh
if ab + bc + ac − a2 − b2 − c2 ≥ 0 if ab + bc + ac − a2 − b2 − c2 < 0
t+1 t φnew i,j,k = min(φi,j,k , φi,j,k )
(9)
t The solution process consists in always using φnew i,j,k to replace the old φi,j,k values in the sweeping calculation process. A more detailed description of the fast sweeping algorithm can be found in. ?, 20 Before solving Eqn. 8, we first assign exact values at Cartesian grid points in or near Γ (Fig. 1(c)) to enforce the boundary condition. We then apply the fast sweeping algorithm to calculate the φ field for the remainder of the domain.
C.
Propagate parameters
The propagation equation is ~xt =
∇φ |∇φ|
(10)
replacing ~xt by a finite difference scheme, Eqn. (10) can be rewritten as ~xn+1 = ~xn ± The term ∇φ |∇φ|
∇φ |∇φ|
∇φ ∗ dt |∇φ|
(11)
can be viewed as a unit propagation speed in the normal direction. For an outward propaga-
tion, is positive; for an inward propagation, into three steps:
∇φ |∇φ|
is negative. The propagation process can be divided
5 of 15 American Institute of Aeronautics and Astronautics
cbr1
cbr2
Figure 2. Stop propagation when φn+1 − φn < 0
1. calculate φx , φy , and φz We first calculate φx , φy , and φz at each Cartesian grid node using a two-sided finite difference scheme for internal nodes, and one-sided finite difference scheme for boundary nodes. Then, we calculate φ x , φy , and φz for all the points to be propagated at Γ using interpolation (Fig. 3(a)). 2. propagate points on Γ We use a forth-order Runge-Kutta method? to obtain the solution to Eqn. (11). We stop propagation once we detect that the value of φ at position ~xn+1 is greater or equal to the given propagation distance. If the value of φ at the final position ~xn+1 is greater than the given propagation distance, we use linear interpolation to get the exact position ~xf inal (Fig. 3(b)). 3. construct Γp Finally, we construct Γp by sequentially connecting all points on Γp using the parameterization information conveyed from Γ.
III.
Results
We verify the effectiveness of this scheme from two aspects: • Accuracy, in which we verify that the error for the propagated front has first order accuracy. • Efficiency, in which we verity that the complexity of this scheme is Ω(N ). We use two test cases to verity this algorithm, a circle, which is centered at (0.0), radius is 0.5; and a cube, which is centered at (0.0.0), side length is 0.5. We propagate these test cases to their inward direction. The offset distance is 0.1. The time increment is 0.01. We fix the number of discretizing points at front, and vary the Cartesian grid size from 202 to 3202 for the circle, from 203 to 1603 for the cube. For the circle, the front is discretized into 40 points; for the cube, the front is discretized into 82308 triangles. These tests were carried out on a standard desktop PC running Linux, and equiped with an AMD Athlon processor running at 1400 MHz, and a memory size of 512Mb. A.
Accuracy Validation
To measure accuracy and order of convergence of the method, we use the L 2 and L∞ norms21 to measure the error on the position of the propagated front. As shown in the results presented, our method converges 6 of 15 American Institute of Aeronautics and Astronautics
(φ x1, φ y1)
(φ x2, φ y2)
(x0, y 0)
(φ x , φ y )
Γ (φ x4, φ y4)
(φ x3, φ y3)
(x1, y1) (xfinal, y final) (x2 , y 2)
(a) We first calculate (φxi , φyi ) at each Cartesian grid node using a finite difference scheme, and then interpolate (φx , φy ), which are the points located on Γ, using (φx1 , φy1 ), (φx2 , φy2 ), (φx3 , φy3 ) and (φx4 , φy4 )
(b) If φ ≥ d at the last propagation point (x2 , y2 ), then calculate the final position of (xf inal , yf inal ) using a linearly interpoled φ value between (x1 , y1 ) and (x2 , y2 ). So φ = d at (xf inal , yf inal ).
Figure 3. Propagation process
for increasing mesh size with first order accuracy. Table 1 illustrates the comparisons of L 2 and L∞ for inwardly propagating a circle; Table 2 illustrates the comparisons of L2 and L∞ for inwardly propagating a cube. From Fig. 4 and Fig. 6, we can see that both L2 and L∞ converge when we increase the mesh size. Fig. 5(a) illustrates the result of propagating a circle in its inward direction using a coarse Cartesian mesh (mesh size=202); Fig. 5(b) illustrates the same propagation of a circle using a refined Cartesian mesh (mesh size=3202). Fig. 7(a) shows the result of propagating a cube in its inward direction using a coarse Cartesian mesh (mesh size=203); and Fig. 7(b) is the result for propagating the same cube in its inward direction using a refined Cartesian mesh (mesh size=1603). From these figures, we can see that more precise offset surface definitions are obtained when we increase the Cartesian mesh density. From Fig. 5 and Fig. 7, we can intuitively observe that mesh size severely influences the propagated front. The output circle in Fig. 5(b) looks smoother than the circle in Fig. 5(a), and the behavior for sharp corners in the output cube in Fig. 7(b) looks much better than the cube in Fig. 7(a). As shown by these experiments, the propagated front position directly depends on mesh size, since the ∇φ unit propagation velocity in the normal direction, |∇φ| , is entirely dependent on the precision of the φ field, which is the numerical solution to Eqn.(8). The fast sweeping scheme used to solve Eqn.(8) is first order accurate with respect to mesh size, this dominates the propagation accuracy. The convergence of the φ field has been discussed in more details in the previous work.? 1.
Circle Mesh Size L2 L∞
202 0.00497535 0.00548393
402 0.00315108 0.00367538
802 0.00244957 0.00319467
1602 0.00182176 0.00246195
Table 1. Norms comparison (inwardly propagate a circle)
7 of 15 American Institute of Aeronautics and Astronautics
3202 0.00136155 0.00172537
Nomrs comparision for circle(inward propagation) 0.01
log(norm)
norm2 of propagated circle norminfty of propagated circle
0.001 0.001
0.01
0.1
1
log(dh)
Figure 4. Norms comparison (inwardly propagate a circle)
2.
Cube Mesh Size
203
403
803
1603
L2 L∞
0.0476718 0.17472
0.0393768 0.128057
0.0368557 0.100638
0.036115 0.0862622
Table 2. Norm Comparison for the propagated surface (Cube Inward Propagation)
B.
Efficiency Validation
In order to evaluate the efficiency of the scheme, we compare the timing results obtained for the calculations presented in our previous work? with those performed specifically using the new propagation method. The total calculation time includes three parts: 1. initializing time, which is the time used to initialize the exact distance for the points near or on front; 2. fast sweeping time, which is the time used to solve numerical solution φ to equation Eqn. 8; 3. iso-front capturing time or propagation time iso-front capturing time is the time used to capture iso-front, which is called propagated front in this paper. This algorithm was used in our previous work.? propagation time is the time used to propagate all the points located on the original front to their destination positions. This algorithm is used in the present work. Thus the total time consumed is calculated either by the formula of total time = initializing time + f ast sweeping time + iso − f ront capturing time 8 of 15 American Institute of Aeronautics and Astronautics
(12)
cbr
cbr
(a) mesh size = 202
(b) mesh size = 3202
Figure 5. Propagate a circle (inward direction)
in our old algorithm,? or calculated by the formula of total time = initializing time + f ast sweeping time + propagation time
(13)
in this work. To compare the efficiency, we will measure these times consumed for each part separately. Table 3 illustrates the time consumed for each part to propagate a circle in its inward direction; Table 4 illustrates the time consumed for each part to propagate a cube to its inward direction. Fig. 8 displays the time slope consumed in each part for inwardly propagating a circle, and Fig. 9 displays the time slope consumed in each part for inwardly propagating a cube. From Fig. 8 and Fig. 9, we can see that initializing time, fast sweeping time and iso-surface capturing time are linearly dependent on the mesh size as discussed in our previous paper;? propagation time is dependent on the number of points for front initialization and time increment. Once these two factors are fixed, propagation time is a constant value. The reason for this is that in the propagation process, we start our propagation from the original front. Once the number of initial points in the original front is fixed, the next factor to influence the propagation time is time increment dt. If dt is too small, the points on Γ have to be propagated many times to reach their destinations. Thus the propagation time is not be influenced by mesh size. It is only a function to the number of original front points and time increment. 1.
Circle Mesh Size
202
402
802
1602
3202
Initializing Time Fast Sweeping Time iso-surface capturing time Propagation Time
0.00457 0.003768 0.000212 0.016103
0.007592 0.023966 0.000712 0.016017
0.015335 0.160162 0.001863 0.015958
0.040478 1.2056 0.005838 0.016986
0.131439 9.11685 0.022492 0.017038
Table 3. Time Comparison for Circle (Inward Propagation)
2.
Cube
9 of 15 American Institute of Aeronautics and Astronautics
Nomrs comparision for cube(inward propagation)
log(norm)
1
norm2 of propagated surface norminfty of propagated surface
0.1
0.01 0.01
0.1
1
log(dh)
Figure 6. Norms comparison (inwardly propagate a cube)
Mesh Size
203
403
803
1603
initializing time Fast Sweeping time Propagation Time 1 Propagation Time 2
1.15137 0.044765 0.01799 0.088228
2.87791 0.54249 0.418908 0.082101
13.9404 9.05241 8.01082 0.082883
97.7702 162.96 147.378 0.088191
Table 4. Time Comparison for a Cube (Inward Propagation)
10 of 15 American Institute of Aeronautics and Astronautics
(a) mesh size 203
(b) mesh size 1603
Figure 7. propagate a cube (inward direction)
Times comparision for circle (inward propagation) 1000
initializing time (s) fast sweeping time (s) iso-surface capturing time (s) propagation time (s)
100
log(time)
10
1
0.1
0.01
0.001
0.0001 0.001
0.01
0.1 log(dh)
Figure 8. Propagate a circle (inward direction)
11 of 15 American Institute of Aeronautics and Astronautics
1
Times comparision for cube(inward propagation) 1000
initializing time (s) fast sweeping time (s) iso-surface capturing time (s) propagation time (s)
100
log(time)
10
1
0.1
0.01 0.01
0.1 log(dh)
Figure 9. Propagate a cube (inward direction)
12 of 15 American Institute of Aeronautics and Astronautics
1
IV.
Method characteristics and Applications
The major advantages of this method is that the propagation lines will not intercross each other when proper time step t is used (as shown in Fig. 10); two-dimensional parameterization stored in Γ are entirely carried out to the propagated front.
Figure 10. Propagation lines at sharp conner
This method can be applied to many types of surfaces, i.e. parameterized surfaces, tessellated surfaces, trimmed surfaces. It also supplies a better flexibility to offset partial surfaces for a geometry. It has a strong ability to control the offset distance to obtain uneven distanced surfaces. For exemple, we can establish a certain connection between local parameters (ui , vi ) of the original surface and propagation distance to control the offset process. This method can be used in many applications, for example, to shrink or enlarge a geometry with any complexity (Fig. 11(a)), to generate a structured grid in the area of boundary layers for viscous flow (Fig. 11(b)). We present results from some applications at hand to illustrate the use of the method on concrete application. Fig. 12(a) is a mesh generation example in a thin layer around the body of a diffuser; Fig. 12(b) is an example to offset the body of a diffuser in the inward direction to generate a thin layer around the body of diffuser.
(a) Offsetting front
(b) Generating mesh for boundary layers
Figure 11. Examples of application
13 of 15 American Institute of Aeronautics and Astronautics
(a) Mesh generation result
(b) Surface propagation result
14 of 15 American Institute of Aeronautics and Astronautics
V.
Conclusion
The main features of the proposed method are : • The parameterized information stored in the original geometry front can be entirely transformed to the offset curves and surfaces. • The local normal directions at front points will not inter-cross each other. • This technique can be applied to many types of front, i.e. a tesselated front, parameterized front, trimmed surfaces, etc. • The accuracy of this method is first order to mesh size. • The complexity of this method is Ω(N ), N is the total number of Cartesian grid nodes. The method can be extended to: (1) generate structured grids near boundary area for studying viscous flow. It guarantees that the mesh is orthogonal to its geometric front, since the propagating paths are developed along local geometric normal directions; (2) generate geometric information, for example, to shrink or enlarge a geometry; and (3) generate an uneven distanced skin by properly setting the offset distance function.
References 1 Maekawa, 2 Guibault,
T., “An Overview of Offset Curves and Surfaces,” Computer-Aided Design, Vol. 31, No. 3, 1999, pp. 165–173. F., Labb, P., Garon, A., and Camarero, R., “Automatic Topologically Based Blocking for Hybrid Grid Gener-
ation,” . 3 Kumar, R., Shastry, K., and Prakash, B., “Computing offsets of trimmed NURBS surfaces,” Computer-Aided Design, Vol. 35, No. 5, 2003, pp. 411–420. 4 Pirzadeh, S., “Unstructured Viscous Grid Generation by Advancing-Layers Method,” No. AIAA-93-3453, Monterey, CA, August. 5 Hoschek, J., “Offset Curves in the Plane,” Computer-Aided Design, Vol. 17, No. 2, 1989, pp. 77–82. 6 Farouki, R., “The approximation of non-degenerate offset surfaces,” Computer Aided Geometric Design, Vol. 3, No. 1, 1986, pp. 15–43. 7 Blomgren, R., “B-Spline Curves,” 1981, Boeing Document, Class Notes, B-7150-BB-WP-281/D-441.2. 8 Tiller, W. and Hanson, E., “Offsets of Two-Dimensional Profiles,” IEEE Computer Graphics and Applications, Vol. 4, No. 9, 1984, pp. 61–69. 9 Coquillart, S., “Computing offsets of B-spline curves,” Computer-Aided Design, Vol. 19, No. 6, 1987, pp. 305–309. 10 Nachman, L., “Approximating Offset Surfaces in a NURBS Environment,” Second SIAM Conference on Geometric Design, Tempe, Arizona, November 1991. 11 Piegl, L. A. and Tiller, W., “Computing offsets of NURBS curves and surfaces,” Computer-Aided Design, Vol. 31, No. 2, February 1999, pp. 147–156. 12 Sun, Y., Nee, A., and Lee, K., “Modifying free-formed NURBS curves and surfaces for offsetting without local selfintersection,” Computer-Aided Design,In Press, 2004. 13 Sethian, J., An Analysis of Flame Propagation, Ph.d. dissertation, University of California, Berkeley, CA, 1982. 14 Sethian, J., “Curvature and the Evolution of Fronts,” Comm. in Math. Phys., Vol. 101, 1985, pp. 487–499. 15 Osher, S. and Sethian, J., “Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations,” Journal of Computational Physics, Vol. 79, 1988, pp. 12–49. 16 Sethian, J., “A Fast Marching Level Set Method for Monotonically Advancing Fronts,” Proc. Nat. Acad. Sci, Vol. 93, No. 4, 1996, pp. 1591–1595. 17 Sethian, J., Level set methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, Cambridge University Press, 1999. 18 Kimmel, R. and Bruckstein, A., “Shape offsets via level sets,” Computer-Aided Design, Vol. 25, No. 3, 1993, pp. 154–162. 19 Tsai, Y., Cheng, L., Osher, S., and Zhao, H., “Fast Sweeping Algorithms for a Class of Hamilton-Jacobi Equations,” SIAM Journal on Numerical Analysis, Vol. 41, No. 2, 2003, pp. 673–694. 20 Zhao, H., “A fast sweeping method for Eikonal equations,” to appear in Math. Comp. (2004), May 2004. 21 Fortin, A. and Garon, A., Les ´ ´ el´ ements finis: de la th´ eorie ` a la pratique, Ecole Polytechnique de Montr´eal, 2000.
15 of 15 American Institute of Aeronautics and Astronautics
