- Email: [email protected]
Orthogonal and conformal surface grid generation
Applied Numerical Mathematics 46 (2003) 249–262 www.elsevier.com/locate/apnum
Orthogonal and conformal surface grid generation Renzo Arina Dip. di Ingegneria Aeronautica e Spaziale, Politecnico di Torino, C.so Duca degli Abruzzi 24, 10129 Torino, Italy
Abstract Surface grid-generation algorithms are based on a projection strategy, consisting in generating a grid on the surface parametric plane, and then to project it on the surface. However during the projection step the grid characteristics, such as point distribution and orthogonality, are often not correctly transported. To overcome these difficulties, a new approach is proposed which can be divided into two main steps: (i) the parameterization of the surface by isothermal coordinates, (ii) the generation of a two-dimensional grid on the conformal parametric plane and its projection on the surface. In this way the projection is based on a conformal map, preserving angles and scale-length ratios. 2003 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Numerical grid generation; Orthogonal mappings; Quasiconformal mappings
1. Introduction Many surface grid-generation algorithms are based on a projection strategy. This approach consists in generating a grid on the surface parametric plane, and then to project it on the surface. In this way, it is possible to obtain surface grids with a reasonable computational effort. However in most cases during the projection step the grid characteristics, such as point distribution and orthogonality, are often not correctly transported. This fact makes the resulting surface grid strongly influenced by the surface local topology. In this paper we present a possible way to overcome this problem. It is known that if the surface is parameterized by isothermal coordinates, obtained by a conformal transformation of the surface domain onto a rectangle, the metric of any grid generated on the conformal parametric plane will be correctly projected on the surface. This statement is based on the fact that a conformal map preserves angles and scale-length ratios. In other words the resulting surface grid will display the appropriate clustering and orthogonality requirements, independently from the surface local topology. E-mail address: [email protected] (R. Arina). 0168-9274/03/$ – see front matter 2003 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/S0168-9274(03)00045-X
250
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
The present approach can be divided into two main steps: (i) the parameterization of the surface by isothermal coordinates, (ii) the generation of a two-dimensional grid on the conformal parametric plane and its projection on the surface. The first step is the most expensive part, but it is important to remark that it must be accomplished only once for a given surface geometry. Once the parameterization is obtained, it can be used as surface data set on which the user will work. The cost of the second step will depend upon the specific grid generation technique employed on the parametric plane. The projection consists in a simple algebraic interpolation and consequently its cost is negligible. In the next two sections, we review the problem of introducing isothermal coordinates on a surface. In Section 4 we show how a grid, generated by an algebraic method can be efficiently projected on the surface. Because of the angle preserving property of the conformal projection, it is worthwhile to develop a technique capable of generating orthogonal grids on the parametric plane, and to project them conformally on the surface. This issue is addressed in the last section of the paper, where it is described as a particular case of the more general problem of two-dimensional orthogonal grid generation. It is also proved that it is possible to generate orthogonal grids specifying the boundary point distribution along the entire surface boundary.
2. Conformal and orthogonal surface coordinate transformations Conformal and orthogonal mappings are particular cases of a wider family of mappings termed quasiconformal [1,6]. The theory of quasiconformal mappings on surfaces is well developed and examples on its application in grid generation can be found in [7]. In this section we will follow and extend an other application of quasiconformal mappings reported in [4]. The most general set of equations which can be used for grid generation in three-dimensional Euclidean spaces, for r¯ = x i e¯i = ξ i g¯i , where {x i } denotes Cartesian coordinates, {ξ i } a general set of curvilinear coordinates, {e¯i } and {g¯i } the tangent vectors along the coordinate lines, is given in tensor notation by g ij
r ∂ 2xr ij k ∂x − g Γ = 0, ij ∂ξ i ∂ξ j ∂ξ k
(1)
with i, j, k, r = 1, 2, 3, where gij are the components of the metric tensor and Γijk the Christoffel symbols corresponding to the curvilinear coordinate system, and g ik gkj = δji . Details on the derivation of Eq. (1) are reported in [3]. A surface in the Euclidean space is individuated by two coordinates, say ξ 1 and ξ 2 , holding the third fixed. If in addition the ξ 3 -coordinate lines are orthogonal to this surface (g13 = g23 = 0), system (1) reads, with r = 1, 2, 3 and α, β, γ = 1, 2, g αβ
r ∂ 2xr αβ γ ∂x − g Γ = (K1 + K2 )nr , αβ ∂ξ α ∂ξ β ∂ξ γ
(2)
where K1 and K2 are the principal curvatures of the surface, and ni e¯i the local unit normal vector forming with the tangent vectors of the surface coordinate lines ξ 1 and ξ 2 , a right-handed frame. The same system could also be obtained starting from the Gauss and Weingarten equations for a surface, as shown by Warsi [9].
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
251
Eq. (2) is valid for any general surface curvilinear coordinate system {ξ α }. The specification of a particular mapping is obtained by imposing an appropriate form of the metric tensor components gαβ relative to the curvilinear system under study. In three dimensions the symmetric tensor gαβ has six independent terms, being gij = gj i if i = j . If the three-dimensional space has to be Euclidean, as in our case, only three of them are independent, and consequently we have only three degrees of freedom in fixing the metric tensor components. One condition has already been imposed by stating that g13 = 0. Hence we still have two conditions for specifying g11 , g12 and g22 . For an orthogonal surface grid {ξ α } we have g12 = 0, then Eq. (2) take the form, with F 2 = g22 /g11 , ξ 1 ≡ ξ and ξ 2 ≡ η, ∂x i ∂ 1 ∂x i ∂ F + = (K1 + K2 )ni . (3) ∂ξ ∂ξ ∂η F ∂η The last condition is employed in specifying the function F (ξ, η), that is the scale-factor ratio. For each specific application the control of the grid is then obtained by acting on F (ξ, η). In the case of conformal mappings we have F = 1. System (3) is formed by three equations, but it is possible to reduce it to a system of two equations by introducing an arbitrary surface coordinate system {uα } such that x i (uα ), with surface metric tensor α . By the chain rule for differentiating composite functions components aαβ , and Christoffel symbols Γβγ 1 2 and considering u and u as dependent variables, from (3) we obtain, with α, β, γ = 1, 2, ∂uα ∂ 1 ∂uα 1 ∂uγ ∂uδ ∂uγ ∂uδ ∂ α F + = −Γ γ δ F + . (4) ∂ξ ∂ξ ∂η F ∂η ∂ξ ∂ξ F ∂η ∂η An alternative way for deriving Eq. (4), starts directly from the constraints on the surface metric tensor components gαβ with respect to the orthogonal coordinates {ξ α }, which are rewritten here g22 = F 2 (ξ, η)g11
and
g12 = g21 = 0,
(5)
where the function F is non-zero and continuously differentiable. If aγ δ are the metric tensor components with respect to the uα -coordinates, by the transformation law for tensor components, ∂uγ ∂uδ , ∂ξ α ∂ξ β Eqs. (5) yield the relations aγ β ∂uγ ∂uα = ε αβ √ , F ∂ξ a ∂η gαβ = aγ δ
(6)
(7)
2 and ε αβ is the permutation symbol (ε 11 = ε 22 = 0, ε 12 = −ε 21 = 1). A similar where a = a11 a22 − a12 expression can be obtained for ∂uα /∂η, and precisely
√ ∂uβ ∂uα = εβγ F aαγ a . (8) ∂η ∂ξ Eqs. (7) and (8) are the generalized Cauchy–Riemann equations. The reason for this terminology is that they reduce to the well known Cauchy–Riemann equations of complex analysis when the surface is flat ({uα } Cartesian coordinates and aαβ = δαβ ) and the coordinates {ξ α } represent conformal coordinates. Differentiating by ξ Eq. (7), by η Eq. (8), and adding them, we obtain γ ∂uγ ∂ aβγ ∂uα ∂ 1 ∂uα aβγ ∂ αβ ∂u ∂ − . (9) F + =ε √ √ ∂ξ ∂ξ ∂η F ∂η ∂η ∂ξ ∂ξ ∂η a a
252
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
Eqs. (4) and (9) are equivalent, but Eq. (9) are more appropriate as mapping system, because for their γ solution it is not necessary to compute explicitly the Christoffel symbols Γαβ . It is worth to note that Eq. (9) represent the extension to a surface domain of the equations proposed by Ryskin and Leal [8], for the numerical generation of two-dimensional orthogonal grids in Euclidean (that is flat) domains. On flat surfaces it is possible to introduce Cartesian coordinates {uα }, in this case the r.h.s of Eq. (9) vanishes. The same is true for the r.h.s of Eq. (4). Hence the system takes the same form as the mapping system proposed in [8]. Following Ryskin and Leal we call function F (ξ, η), the distortion function.
3. Surface conformal mapping Surface grid generation consists in the mapping of a domain D, defined on the surface x i (uα ), onto the rectangular domain D defined on the mathematical plane ξ α . The coordinate transformation uα (ξ β ) is the solution of the boundary-value problem formed by a mapping system obtained from the general differential system (2), by specifying appropriate constraints on the curvilinear metric tensor components gij , and by a set of boundary conditions. In the present case the constraints of orthogonality yield the mapping system (9). We are interested in parameterizing the surface with isothermal coordinates. That is, we want to map conformally D over the plane D; hence, requiring g22 = g11 , F must be constant and equal to unity. The mapping uα (ξ β ) will then be an analytic function. The boundary conditions must ensure that the image of the boundary ∂D coincides with the boundary ∂D on the surface. This correspondence can be obtained by specifying Dirichlet boundary conditions. However the conditions of analyticity (7), (8) must hold through the boundary ∂D, consequently it is not possible to fix an arbitrary distribution of the points along the boundary at the same time, being the problem over determined. It is possible to overcome this problem by representing the boundary ∂D in parametric form, and to enforce a ‘shape correspondence’ of ∂D with the image of ∂D, by moving the boundary grid points along the boundary ∂D, between the corresponding corner points, in order to satisfy Eqs. (7), (8). If the mapping uα = f (ξ β ) represents an allowable mapping, the Jacobian determinant of the function f : D → D must not vanish on D. In the previous section it has been shown that conformal mappings satisfy in each point of D the Cauchy–Riemann equations (7), then it follows that the Jacobian determinant reads aβγ ∂uγ ∂uβ ∂uα ∂uβ =√ . (10) J = εαβ ∂ξ ∂η a ∂η ∂η γ
are zero. This quadratic form is positive definite, and it vanishes only if both the partial derivatives ∂u ∂η But this cannot be true for a regular solution, then the grid will always be unfolded. The existence of surface conformal mappings, solutions of the boundary value problem formed by Eq. (9) and the boundary-point relocation procedure ensuring analytic continuation at the boundaries, can be proved as follows. It can be seen that the problem of finding isothermal parameters is reduced to finding a solution of the Beltrami equation, and that for F (ξ, η) = 1, this solution exists [1]. Isothermal coordinates define a natural conformal structure for an orientable surface, which thus becomes a Riemann surface [5]. By the Koebe uniformization theorem for Riemann surfaces, it follows that an arbitrary region D of the surface bounded by a simple closed curve (Jordan domain), can be conformally mapped
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
253
onto a simply-connected plane domain with one-to-one correspondence. The Koebe theorem states that the surface Jordan domain is conformally equivalent to a plane domain. Moreover the theorem on correspondence of boundaries states that the conformal mapping of a Jordan domain onto another Jordan domain has a homeomorphic extension to the boundary. For such a mapping, the images of three boundary points can be prescribed arbitrarily on the boundary of the image domain. Then, given the Jordan domains D and D, whose positive orientations are given by the points p1 , p2 , p3 on ∂D, and q1 , q2 , q3 on ∂D, then there is a uniquely determined conformal mapping of D onto D, which carries the boundary points pi into the boundary points qi . This result is not only a proof of the existence of a regular solution, globally one-to-one, but it also provides a geometric insight. Being interested in non-degenerate quadrilaterals domains D, it follows that it is not possible to map D on an a priori specified rectangular region D, with the vertices coinciding pairwise. All quadrilaterals are therefore divided into conformal equivalence classes. Every quadrilateral possesses a canonical mapping, transforming four distinguished points to the vertices of a rectangle, and every conformal equivalence class of quadrilaterals contains all similar rectangles characterized by the same ratio of sides M, termed conformal module. For a rectangular region, the conformal module coincides with the side ratio, then it is possible to normalize the region D into the unit square D ! and to treat M(D) as an unknown stretching parameter. The conformal module of the region D is usually unknown. Moreover it is not an easy task to compute it exactly. As pointed out in [7], in many cases a reasonable good estimate is sufficient. A possibility is to update the estimate during the iterative procedure by the relation ! ! g /g22 dξ ! dη! ! , (11) M 2 = D 11 ! ! ! ! D ! g22 /g11 dξ dη obtained by Eqs. (7), (8), where gij! are the metric tensor components with respect to the coordinates {ξ !α } defined on D ! . It is possible to arrive at an expression for M through an extremal problem, by use of a length-area method [6]. Considering the canonical mapping f of the quadrilateral D, with distinguished points (z1 , z2 , z3 , z4 ), z being a complex variable, onto the rectangle D = {u + iv | 0 < u < a, 0 < v < b}. Then 2 f (z) dx dy = ab. (12) D
Let Γ be the family of all locally rectifiable Jordan arcs in D which join the sides (z1 , z2 ) and (z3 , z4 ). Then f (z) |dz| b (13) γ
for every γ ∈ Γ , with equality if γ is the inverse image of a vertical line segment of D joining its horizontal sides. Hence (z)|2 dx dy D |f . (14) M= (infγ ∈Γ γ |f (z)||dz|)2 Eqs. (9), (11) and the boundary point relocation procedure form the boundary-value problem whose solution represents a conformal surface mapping. For its numerical solution, we have adopted an iterative
254
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
Fig. 1. Conformal coordinates on a paraboloid.
procedure, where the basic iteration consists of two steps. In the first one, for a fixed distribution of the boundary points and given conformal module, the finite-difference discretizations of Eq. (9) are solved by an approximate factorization technique. Then the positions of the boundary points are adjusted in order to satisfy Eq. (7). The treatment of the boundary points, explained in [2] for the case of plane domains, can be directly extended to the present case by applying the planar procedure on the local tangent plane of the surface at the points of the boundary line ∂D and individuated by the two tangent vectors of the coordinate lines {ξ α }. Before computing again a new iterate for the interior points, the estimate of the conformal module is updated by using Eq. (11). In Figs. 1 and 2, examples of conformal coordinates on an paraboloid and on a canopy-like configuration are shown.
Fig. 2. Conformal coordinates on a canopy-like surface.
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
255
4. Algebraic surface grids We analyze now the second step of the approach. That is the generation of a grid on the surface by applying a two-dimensional grid algorithm for planar domains, on the domain D, the isothermal parametric rectangle of aspect ratio M(D), and then its projection on the surface. This step can be repeated several times, generating different grids, corresponding to different control functions or different boundary-point distributions, but the parametric rectangle will remain unchanged. The specification of the boundary-point distribution can be made on the surface, where we have a direct control. Being the parametric plane related to the surface domain by the conformal mapping, the corresponding point distribution on D is automatically obtained. Once the grid has been generated, the projection consists in finding the functions {x i }, representing the surface geometry, at the new grid point positions. This operation can be made very efficiently, because the background conformal grid {ξ α } represents a Cartesian net on D, and as a consequence, the procedure for finding the location of the new grid points with respect to the conformal grid points, is trivial. Knowing the new point location on D, the value of the functions x i can be found by interpolation. In Figs. 3 and 4, we display the grids obtained on the paraboloid and canopy-like geometries by a transfinite interpolation, with control of orthogonality at the boundaries by means of Hermite interpolants. It is interesting to note how the orthogonality along the boundaries, obtained on D, has not been destroyed by the conformal projection. Similarly the specified grid clustering is not affected by the surface topology. The projection maps are based on the isothermal coordinates shown in Figs. 1 and 2.
5. Orthogonal surface grid generation In this section we describe a method for generating orthogonal grids, with fixed boundary points distribution. We will apply it to the mapping between the rectangle D, of aspect ratio M(D),
Fig. 3. Algebraic grid on a paraboloid.
256
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
Fig. 4. Algebraic grid on a canopy-like surface.
and a unit square D ! . Being now D ! our computational domain and D the isothermal parametric domain. It must be pointed out that the approach presented in this section, does not apply only to the particular case of rectangular domains but it is valid for general simply-connected non-degenerate quadrilaterals. Then, more generally, we will describe a grid-generator mapping between the twodimensional simply connected planar domain D, representing the physical non-degenerate quadrilateral, and the computational plane D, usually the unit square. 5.1. Quasiconformal mappings From Eq. (9), we have that an orthogonal coordinate transformation is obtained as a solution of the system ∂x i ∂ 1 ∂x i ∂ F + =0 (15) ∂ξ ∂ξ ∂η F ∂η with {x i } Cartesian coordinates, and i = 1, 2. As pointed out in [8], there is a strong connection between orthogonal mappings and quasiconformal mappings. Introducing the complex planes z = x + iy and ζ = ξ + iη, the mapping functions x(ξ, η) and y(ξ, η), can be interpreted as the real and imaginary part of a complex function z = f (ζ ) = x(ζ ) + iy(ζ ), defined on the computational domain D ⊂ ζ , and with image on the physical domain D ⊂ z. Let dso2 = gij dξ i dξ j be a given metric defined on D, it can be expressed in terms of the complex variables ζ and ζ¯ = ξ − iη, as follows 2 2 (16) dso2 = ρ(ξ, η) dζ + µ(ξ, η) dζ¯ , where
√ g11 + g22 + 2 g , |ρ| = 4 2
(17)
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
and g11 − g22 + 2ig12 µ= √ , g11 + g22 + 2 g
√ g11 + g22 − 2 g |µ| = √ 1. g11 + g22 + 2 g 2
Similarly, on D, for the Cartesian metric ds 2 = dx 2 + dy 2 we have, recalling that z = f (ζ ), fζ¯ 2 2 2 2 ¯ ds = |dz| = |fζ | dζ + dζ , fζ
257
(18)
(19)
with fζ = 12 (fξ − ifη ) and fζ¯ = 12 (fξ + ifη ) being the formal complex derivatives. So it appears that if we require, for fζ = 0, fζ¯ = µfζ ,
(20)
then ds 2 = λ(ζ ) dso2 ; that is the mapping z = f (ζ ), is conformal with respect to a given metric gij . The connection between Eq. (15) and the Beltrami equation (20) follows by showing that Eq. (20) corresponds to the Cauchy–Riemann equations. In fact, let µ = µr + iµi , from (20) we have (1 − µr )xξ = (1 + µr )yη − µi (yξ − xη ), (1 + µr )xη = −(1 − µr )yξ + µi (xξ + yη ), let 1−F , µi = 0, (21) 1+F we have the generalized Cauchy–Riemann equations (from now on referred as GCR), which are the plane analog of Eqs. (7) and (8), µr =
F xξ = yη , xη = −Fyξ .
(22) (23)
Hence, solving system (15) amounts to solve the Beltrami equation for µ real. From the theory of quasiconformal mappings [6,1], every homeomorphic function which is a solution of Eq. (20), with K −1 1, (24) K +1 holding at every point of D, is a K-quasiconformal mapping, and µ is termed the complex dilatation. The connection between orthogonal mapping and K-quasiconformal mappings is more evident considering the geometric definition of K-quasiconformal mappings. Geometrically a quadrilateral D is the conformal image of a rectangle D, if the ratio of the sides of D is equal to the modulus of D. If z = f (ζ ) maps D onto D, with modulus m(D) = M(D), then f is said to be K-quasiconformal if the maximal dilatation K is finite and M(D) 1 K. (25) K m(D) |µ| =
In the case of grid generation, D is usually the unit square, then m(D) = 1, and for non-degenerate quadrilaterals M(D) is always finite. Moreover if µ is real and continuous, as in our case, the mapping z = f (ζ ) is a diffeomorphism, and then it represents a coordinate transformation mapping.
258
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
For quasiconformal mappings, it is possible to state a generalized Riemann mapping theorem [1]. In fact for µ regular, the metric gij defines a structure of Riemann surface on the domain D, and the general uniformization theorem asserts the existence of a conformal mapping. Moreover under suitable boundary conditions, the quasiconformal mapping can be extended to the mapping between the closures D and D. 5.2. Distortion function The previous characterization of quasiconformal mappings as conformal mapping between a
In the frame of grid generation,
with a metric given by gij , and D, implies M(D) = M(D). rectangle D, we are interested in finding a transformation between the unit square D and the quadrilateral D. Then we
→ D is can see that the quasiconformal mapping x i = f (ξ j ) will be of the form f = fˆ ◦ h, where fˆ : D conformal. It follows that the choice of a distortion function F (ξ, η) will amount to specify the mapping
And the relation is given by the expressions (18), (21) between the distortion function F , h : D → D.
the complex dilatation µ and the curvilinear metric gij induced by the mapping h on the rectangle D. The following cases can be analyzed: (a) If the function h is the identity map, gij = δij , then F = 1. We are in the case of a conformal transformation between D and the rectangle D. Of course the two domains must belong to the same conformal equivalence class.
(b) If the function h is represented by the transformation, for (ξˆ , η) ˆ defined on D, ξˆ = M(D)ξ, ηˆ = η; 1 . The transformation f is quasiconformal then we have g11 = M(D)2 , g12 = 0, g22 = 1, and F = M(D) with K = M(D). This case corresponds to the calculation of conformal mappings with unknown conformal module, as described in the previous sections (see also [2]). (c) Another possibility is to represent the function h as composed by the one-dimensional stretching
functions, always with (ξˆ , η) ˆ defined on D,
ξˆ = h1 (ξ ), ηˆ = h2 (η); so g11 = (h1ξ )2 , g12 = 0, and g22 = (h2η )2 . From relations (18) and (21) we have that µ= F =
h1ξ − h2η h1ξ + h2η h2η h1ξ
,
.
This case, combined with the conformal module stretching (case (b)), has been adopted for orthogonal mappings with one-dimensional control [3,4]. (d) It is not possible to specify arbitrary functions hi (ξ, η), because it must be g12 = 0, otherwise µ will not be real. However in [3] it has been shown that it is possible to obtain quasiconformal mappings
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
259
with a more effective control than the one-dimensional stretchings of case (c). This control is obtained specifying along each coordinate line a function governing the arc length distribution between the grid points. This function must be normalized with the condition that the grid points should cover the entire grid line. In this way the functions hi (ξ j ) are more general than simple one-dimensional stretchings, because the distribution functions can vary passing from one to another curve of the same family. Once the function F is specified, system (15) can be solved by adding suitable boundary conditions, in order to respect the requirement of quasiconformal continuation, represented by the GCR equations, and computing the conformal module M(D) by Eq. (11). The numerical algorithm is similar to that described for the calculation of the isothermal coordinates. However in most applications, we are interested in fixing the boundary points. In this case, in order to satisfy the requirement of quasiconformal continuation throughout the boundary, it is not possible to specify the distortion function a priori. Both the distortion function along the boundary and the conformal module must be computed from the boundary distribution of the grid points, in such a way that the GCR are still valid along the boundary. In order to devise a consistent procedure, let introduce the new variables {ξ¯ , η}, ¯ defined by the relations 1 dξ, F dη¯ = dη;
dξ¯ =
then the GCR relations become xξ¯ = yη¯ , xη¯ = −yξ¯ . ¯ let now ω be the differential Defining on D a continuous function g(ζ¯ ) = −y(ζ¯ ) + ix(ζ¯ ) with ζ¯ = ξ¯ + i η, form ω = −y dξ¯ − x dη, ¯ and ω∗ its conjugate, given by ¯ ω∗ = x dξ¯ − y dη. From the GCR relations, it follows that the two differential forms are closed, that is dω = −(xξ¯ − yη¯ ) dξ¯ dη¯ = 0, dω∗ = −(yξ¯ + xη¯ ) dξ¯ dη¯ = 0, then ω is a harmonic form [5]. Moreover for the Stokes’ theorem, dω = dω∗ = 0 on D, implies ω = 0, ω∗ = 0 γ
γ
for any closed curve γ in D. From the previous definitions, it follows that
ω + iω∗ = i(x + iy) dξ¯ + i dη¯ = g(ζ¯ ) dζ¯ .
(26)
260
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
Hence
g ζ¯ dζ¯ =
γ
ω+i
γ
ω∗ = 0,
γ
which implies that g(ζ¯ ) is analytic in the region D. This is a classical result (Morera’s theorem). The complex differential dω + i dω∗ , being zero, is closed, then this form may represent the differential of a function f (ζ¯ ), and from (26) we have df (ζ¯ ) = g(ζ¯ ) dζ¯ .
(27)
Comparing Eqs. (26) and (27), we infer that f (ζ¯ ) = x(ζ¯ ) + iy(ζ¯ ). Being its derivative g(ζ¯ ) analytic, f will be analytic in D. Indeed f (ζ¯ ) represents the conformal mapping between the computational plane D and the domain D. The problem in which we are interested, consists into determining f (ζ¯ ), or equivalently ω(ζ¯ ), from boundary data. From the theorems of existence of Abelian differential forms [5], we have that ω(ζ¯ ) on D is uniquely determined by the condition ω∗ = 0. (28) ∂D
This condition, together with the fact that ω is a closed form, and then on ∂D it must be dω = 0, enable us to compute the conformal module M(D) and the values of F along the boundary, knowing the boundary-point distribution. By Eq. (28), and introducing again the coordinates {ξ, η}, we have that
1 x dξ − y dη = 0, (29) x dξ¯ − y dη¯ = F ∂D
∂D
where ∂D = ∂D1 + ∂D2 (∂D1 with ξ = 0 or 1, ∂D2 with η = 0 or 1). The values of the distortion function F along the boundaries are given by the condition dω = 0 along ∂D, which is satisfied by the GCR relation, which in {ξ, η} coordinates is given by the expression F xξ = yη . Along each side of the computational rectangle D, only yη (on ∂D1 ) or xξ (on ∂D2 ) is given by the boundary distribution. The other comes from the computation of the interior points. Then the boundary values of F must be updated iteratively during the computation, and they must satisfy the integral condition (29). Examples of such a method are given in Figs. 5 and 6, where orthogonal grids have been generated, keeping the same distribution of the boundary points employed in the examples of Figs. 3 and 4. Comparing Figs. 3 and 5, it is evident that the orthogonality in the interior points is more pronounced in the last case (Fig. 5), this because the orthogonality constraint has been imposed in all the domain, and not only along the boundaries.
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
261
Fig. 5. Orthogonal grid on a paraboloid.
Fig. 6. Orthogonal grid on a canopy-like surface.
6. Conclusions In this paper we have presented a novel technique for generating surface grids with specified metric characteristics. The technique consists into two main steps: (i) the parameterization of the surface by isothermal coordinates, (ii) the generation of a two-dimensional grid on the conformal parametric plane and its projection on the surface. In this way, the projection is based on a conformal map, preserving angles and scale-length ratios. It has been proved, by numerical examples, that orthogonal surface grids may be generated on quadrilaterals surface domains, employing well established, and computationally efficient, numerical
262
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
grid generation methods developed for two-dimensional flat domains. In this context we have described the more general problem of two-dimensional orthogonal grid generation, and it has been shown the possibility of generating orthogonal grids specifying the boundary point distribution along the entire domain boundary.
References [1] L. Ahlfors, Quasiconformal mappings and their applications, in: T.L. Saaty (Ed.), Lecture Notes in Modern Mathematics, Vol. II, Wiley, New York, 1964, pp. 151–164. [2] R. Arina, Orthogonal grids with adaptive control, in: J. Hauser, C. Taylor (Eds.), Numerical Grid Generation in CFD, Pineridge Press Ltd., Swansea, 1986. [3] R. Arina, Adaptive orthogonal curvilinear coordinates, in: K. Morton, M.J. Baines (Eds.), Numerical Methods for Fluid Dynamics III, Oxford Univ. Press, 1988, pp. 353–359. [4] R. Arina, Adaptive orthogonal surface coordinates, in: S. Sengupta, J. Hauser, P.R. Eiseman, J.F. Thompson (Eds.), Numerical Grid Generation in CFD ’88, Pineridge Press Ltd., Swansea, 1988, pp. 351–359. [5] H. Cohn, Conformal Mappings on Riemann Surfaces, Dover, New York, 1980. [6] O. Letho, K.I. Virtanen, Quasiconformal Mappings in the Plane, Springer-Verlag, New York, 1973. [7] C.W. Mastin, J.F. Thompson, Quasiconformal mappings and grid generation, SIAM J. Sci. Statist. Comput. 5 (1984) 52–62. [8] G. Ryskin, L.G. Leal, Orthogonal mapping, J. Comput. Phys. 50 (1983) 71–100. [9] Z.U.A. Warsi, Surface grid generation and differential geometry, in: J.E. Castillo (Ed.), Mathematical Aspects of Numerical Grid Generation, SIAM, Philadelphia, PA, 1991, pp. 99–104.
Orthogonal and conformal surface grid generation Renzo Arina Dip. di Ingegneria Aeronautica e Spaziale, Politecnico di Torino, C.so Duca degli Abruzzi 24, 10129 Torino, Italy
Abstract Surface grid-generation algorithms are based on a projection strategy, consisting in generating a grid on the surface parametric plane, and then to project it on the surface. However during the projection step the grid characteristics, such as point distribution and orthogonality, are often not correctly transported. To overcome these difficulties, a new approach is proposed which can be divided into two main steps: (i) the parameterization of the surface by isothermal coordinates, (ii) the generation of a two-dimensional grid on the conformal parametric plane and its projection on the surface. In this way the projection is based on a conformal map, preserving angles and scale-length ratios. 2003 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Numerical grid generation; Orthogonal mappings; Quasiconformal mappings
1. Introduction Many surface grid-generation algorithms are based on a projection strategy. This approach consists in generating a grid on the surface parametric plane, and then to project it on the surface. In this way, it is possible to obtain surface grids with a reasonable computational effort. However in most cases during the projection step the grid characteristics, such as point distribution and orthogonality, are often not correctly transported. This fact makes the resulting surface grid strongly influenced by the surface local topology. In this paper we present a possible way to overcome this problem. It is known that if the surface is parameterized by isothermal coordinates, obtained by a conformal transformation of the surface domain onto a rectangle, the metric of any grid generated on the conformal parametric plane will be correctly projected on the surface. This statement is based on the fact that a conformal map preserves angles and scale-length ratios. In other words the resulting surface grid will display the appropriate clustering and orthogonality requirements, independently from the surface local topology. E-mail address: [email protected] (R. Arina). 0168-9274/03/$ – see front matter 2003 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/S0168-9274(03)00045-X
250
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
The present approach can be divided into two main steps: (i) the parameterization of the surface by isothermal coordinates, (ii) the generation of a two-dimensional grid on the conformal parametric plane and its projection on the surface. The first step is the most expensive part, but it is important to remark that it must be accomplished only once for a given surface geometry. Once the parameterization is obtained, it can be used as surface data set on which the user will work. The cost of the second step will depend upon the specific grid generation technique employed on the parametric plane. The projection consists in a simple algebraic interpolation and consequently its cost is negligible. In the next two sections, we review the problem of introducing isothermal coordinates on a surface. In Section 4 we show how a grid, generated by an algebraic method can be efficiently projected on the surface. Because of the angle preserving property of the conformal projection, it is worthwhile to develop a technique capable of generating orthogonal grids on the parametric plane, and to project them conformally on the surface. This issue is addressed in the last section of the paper, where it is described as a particular case of the more general problem of two-dimensional orthogonal grid generation. It is also proved that it is possible to generate orthogonal grids specifying the boundary point distribution along the entire surface boundary.
2. Conformal and orthogonal surface coordinate transformations Conformal and orthogonal mappings are particular cases of a wider family of mappings termed quasiconformal [1,6]. The theory of quasiconformal mappings on surfaces is well developed and examples on its application in grid generation can be found in [7]. In this section we will follow and extend an other application of quasiconformal mappings reported in [4]. The most general set of equations which can be used for grid generation in three-dimensional Euclidean spaces, for r¯ = x i e¯i = ξ i g¯i , where {x i } denotes Cartesian coordinates, {ξ i } a general set of curvilinear coordinates, {e¯i } and {g¯i } the tangent vectors along the coordinate lines, is given in tensor notation by g ij
r ∂ 2xr ij k ∂x − g Γ = 0, ij ∂ξ i ∂ξ j ∂ξ k
(1)
with i, j, k, r = 1, 2, 3, where gij are the components of the metric tensor and Γijk the Christoffel symbols corresponding to the curvilinear coordinate system, and g ik gkj = δji . Details on the derivation of Eq. (1) are reported in [3]. A surface in the Euclidean space is individuated by two coordinates, say ξ 1 and ξ 2 , holding the third fixed. If in addition the ξ 3 -coordinate lines are orthogonal to this surface (g13 = g23 = 0), system (1) reads, with r = 1, 2, 3 and α, β, γ = 1, 2, g αβ
r ∂ 2xr αβ γ ∂x − g Γ = (K1 + K2 )nr , αβ ∂ξ α ∂ξ β ∂ξ γ
(2)
where K1 and K2 are the principal curvatures of the surface, and ni e¯i the local unit normal vector forming with the tangent vectors of the surface coordinate lines ξ 1 and ξ 2 , a right-handed frame. The same system could also be obtained starting from the Gauss and Weingarten equations for a surface, as shown by Warsi [9].
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
251
Eq. (2) is valid for any general surface curvilinear coordinate system {ξ α }. The specification of a particular mapping is obtained by imposing an appropriate form of the metric tensor components gαβ relative to the curvilinear system under study. In three dimensions the symmetric tensor gαβ has six independent terms, being gij = gj i if i = j . If the three-dimensional space has to be Euclidean, as in our case, only three of them are independent, and consequently we have only three degrees of freedom in fixing the metric tensor components. One condition has already been imposed by stating that g13 = 0. Hence we still have two conditions for specifying g11 , g12 and g22 . For an orthogonal surface grid {ξ α } we have g12 = 0, then Eq. (2) take the form, with F 2 = g22 /g11 , ξ 1 ≡ ξ and ξ 2 ≡ η, ∂x i ∂ 1 ∂x i ∂ F + = (K1 + K2 )ni . (3) ∂ξ ∂ξ ∂η F ∂η The last condition is employed in specifying the function F (ξ, η), that is the scale-factor ratio. For each specific application the control of the grid is then obtained by acting on F (ξ, η). In the case of conformal mappings we have F = 1. System (3) is formed by three equations, but it is possible to reduce it to a system of two equations by introducing an arbitrary surface coordinate system {uα } such that x i (uα ), with surface metric tensor α . By the chain rule for differentiating composite functions components aαβ , and Christoffel symbols Γβγ 1 2 and considering u and u as dependent variables, from (3) we obtain, with α, β, γ = 1, 2, ∂uα ∂ 1 ∂uα 1 ∂uγ ∂uδ ∂uγ ∂uδ ∂ α F + = −Γ γ δ F + . (4) ∂ξ ∂ξ ∂η F ∂η ∂ξ ∂ξ F ∂η ∂η An alternative way for deriving Eq. (4), starts directly from the constraints on the surface metric tensor components gαβ with respect to the orthogonal coordinates {ξ α }, which are rewritten here g22 = F 2 (ξ, η)g11
and
g12 = g21 = 0,
(5)
where the function F is non-zero and continuously differentiable. If aγ δ are the metric tensor components with respect to the uα -coordinates, by the transformation law for tensor components, ∂uγ ∂uδ , ∂ξ α ∂ξ β Eqs. (5) yield the relations aγ β ∂uγ ∂uα = ε αβ √ , F ∂ξ a ∂η gαβ = aγ δ
(6)
(7)
2 and ε αβ is the permutation symbol (ε 11 = ε 22 = 0, ε 12 = −ε 21 = 1). A similar where a = a11 a22 − a12 expression can be obtained for ∂uα /∂η, and precisely
√ ∂uβ ∂uα = εβγ F aαγ a . (8) ∂η ∂ξ Eqs. (7) and (8) are the generalized Cauchy–Riemann equations. The reason for this terminology is that they reduce to the well known Cauchy–Riemann equations of complex analysis when the surface is flat ({uα } Cartesian coordinates and aαβ = δαβ ) and the coordinates {ξ α } represent conformal coordinates. Differentiating by ξ Eq. (7), by η Eq. (8), and adding them, we obtain γ ∂uγ ∂ aβγ ∂uα ∂ 1 ∂uα aβγ ∂ αβ ∂u ∂ − . (9) F + =ε √ √ ∂ξ ∂ξ ∂η F ∂η ∂η ∂ξ ∂ξ ∂η a a
252
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
Eqs. (4) and (9) are equivalent, but Eq. (9) are more appropriate as mapping system, because for their γ solution it is not necessary to compute explicitly the Christoffel symbols Γαβ . It is worth to note that Eq. (9) represent the extension to a surface domain of the equations proposed by Ryskin and Leal [8], for the numerical generation of two-dimensional orthogonal grids in Euclidean (that is flat) domains. On flat surfaces it is possible to introduce Cartesian coordinates {uα }, in this case the r.h.s of Eq. (9) vanishes. The same is true for the r.h.s of Eq. (4). Hence the system takes the same form as the mapping system proposed in [8]. Following Ryskin and Leal we call function F (ξ, η), the distortion function.
3. Surface conformal mapping Surface grid generation consists in the mapping of a domain D, defined on the surface x i (uα ), onto the rectangular domain D defined on the mathematical plane ξ α . The coordinate transformation uα (ξ β ) is the solution of the boundary-value problem formed by a mapping system obtained from the general differential system (2), by specifying appropriate constraints on the curvilinear metric tensor components gij , and by a set of boundary conditions. In the present case the constraints of orthogonality yield the mapping system (9). We are interested in parameterizing the surface with isothermal coordinates. That is, we want to map conformally D over the plane D; hence, requiring g22 = g11 , F must be constant and equal to unity. The mapping uα (ξ β ) will then be an analytic function. The boundary conditions must ensure that the image of the boundary ∂D coincides with the boundary ∂D on the surface. This correspondence can be obtained by specifying Dirichlet boundary conditions. However the conditions of analyticity (7), (8) must hold through the boundary ∂D, consequently it is not possible to fix an arbitrary distribution of the points along the boundary at the same time, being the problem over determined. It is possible to overcome this problem by representing the boundary ∂D in parametric form, and to enforce a ‘shape correspondence’ of ∂D with the image of ∂D, by moving the boundary grid points along the boundary ∂D, between the corresponding corner points, in order to satisfy Eqs. (7), (8). If the mapping uα = f (ξ β ) represents an allowable mapping, the Jacobian determinant of the function f : D → D must not vanish on D. In the previous section it has been shown that conformal mappings satisfy in each point of D the Cauchy–Riemann equations (7), then it follows that the Jacobian determinant reads aβγ ∂uγ ∂uβ ∂uα ∂uβ =√ . (10) J = εαβ ∂ξ ∂η a ∂η ∂η γ
are zero. This quadratic form is positive definite, and it vanishes only if both the partial derivatives ∂u ∂η But this cannot be true for a regular solution, then the grid will always be unfolded. The existence of surface conformal mappings, solutions of the boundary value problem formed by Eq. (9) and the boundary-point relocation procedure ensuring analytic continuation at the boundaries, can be proved as follows. It can be seen that the problem of finding isothermal parameters is reduced to finding a solution of the Beltrami equation, and that for F (ξ, η) = 1, this solution exists [1]. Isothermal coordinates define a natural conformal structure for an orientable surface, which thus becomes a Riemann surface [5]. By the Koebe uniformization theorem for Riemann surfaces, it follows that an arbitrary region D of the surface bounded by a simple closed curve (Jordan domain), can be conformally mapped
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
253
onto a simply-connected plane domain with one-to-one correspondence. The Koebe theorem states that the surface Jordan domain is conformally equivalent to a plane domain. Moreover the theorem on correspondence of boundaries states that the conformal mapping of a Jordan domain onto another Jordan domain has a homeomorphic extension to the boundary. For such a mapping, the images of three boundary points can be prescribed arbitrarily on the boundary of the image domain. Then, given the Jordan domains D and D, whose positive orientations are given by the points p1 , p2 , p3 on ∂D, and q1 , q2 , q3 on ∂D, then there is a uniquely determined conformal mapping of D onto D, which carries the boundary points pi into the boundary points qi . This result is not only a proof of the existence of a regular solution, globally one-to-one, but it also provides a geometric insight. Being interested in non-degenerate quadrilaterals domains D, it follows that it is not possible to map D on an a priori specified rectangular region D, with the vertices coinciding pairwise. All quadrilaterals are therefore divided into conformal equivalence classes. Every quadrilateral possesses a canonical mapping, transforming four distinguished points to the vertices of a rectangle, and every conformal equivalence class of quadrilaterals contains all similar rectangles characterized by the same ratio of sides M, termed conformal module. For a rectangular region, the conformal module coincides with the side ratio, then it is possible to normalize the region D into the unit square D ! and to treat M(D) as an unknown stretching parameter. The conformal module of the region D is usually unknown. Moreover it is not an easy task to compute it exactly. As pointed out in [7], in many cases a reasonable good estimate is sufficient. A possibility is to update the estimate during the iterative procedure by the relation ! ! g /g22 dξ ! dη! ! , (11) M 2 = D 11 ! ! ! ! D ! g22 /g11 dξ dη obtained by Eqs. (7), (8), where gij! are the metric tensor components with respect to the coordinates {ξ !α } defined on D ! . It is possible to arrive at an expression for M through an extremal problem, by use of a length-area method [6]. Considering the canonical mapping f of the quadrilateral D, with distinguished points (z1 , z2 , z3 , z4 ), z being a complex variable, onto the rectangle D = {u + iv | 0 < u < a, 0 < v < b}. Then 2 f (z) dx dy = ab. (12) D
Let Γ be the family of all locally rectifiable Jordan arcs in D which join the sides (z1 , z2 ) and (z3 , z4 ). Then f (z) |dz| b (13) γ
for every γ ∈ Γ , with equality if γ is the inverse image of a vertical line segment of D joining its horizontal sides. Hence (z)|2 dx dy D |f . (14) M= (infγ ∈Γ γ |f (z)||dz|)2 Eqs. (9), (11) and the boundary point relocation procedure form the boundary-value problem whose solution represents a conformal surface mapping. For its numerical solution, we have adopted an iterative
254
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
Fig. 1. Conformal coordinates on a paraboloid.
procedure, where the basic iteration consists of two steps. In the first one, for a fixed distribution of the boundary points and given conformal module, the finite-difference discretizations of Eq. (9) are solved by an approximate factorization technique. Then the positions of the boundary points are adjusted in order to satisfy Eq. (7). The treatment of the boundary points, explained in [2] for the case of plane domains, can be directly extended to the present case by applying the planar procedure on the local tangent plane of the surface at the points of the boundary line ∂D and individuated by the two tangent vectors of the coordinate lines {ξ α }. Before computing again a new iterate for the interior points, the estimate of the conformal module is updated by using Eq. (11). In Figs. 1 and 2, examples of conformal coordinates on an paraboloid and on a canopy-like configuration are shown.
Fig. 2. Conformal coordinates on a canopy-like surface.
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
255
4. Algebraic surface grids We analyze now the second step of the approach. That is the generation of a grid on the surface by applying a two-dimensional grid algorithm for planar domains, on the domain D, the isothermal parametric rectangle of aspect ratio M(D), and then its projection on the surface. This step can be repeated several times, generating different grids, corresponding to different control functions or different boundary-point distributions, but the parametric rectangle will remain unchanged. The specification of the boundary-point distribution can be made on the surface, where we have a direct control. Being the parametric plane related to the surface domain by the conformal mapping, the corresponding point distribution on D is automatically obtained. Once the grid has been generated, the projection consists in finding the functions {x i }, representing the surface geometry, at the new grid point positions. This operation can be made very efficiently, because the background conformal grid {ξ α } represents a Cartesian net on D, and as a consequence, the procedure for finding the location of the new grid points with respect to the conformal grid points, is trivial. Knowing the new point location on D, the value of the functions x i can be found by interpolation. In Figs. 3 and 4, we display the grids obtained on the paraboloid and canopy-like geometries by a transfinite interpolation, with control of orthogonality at the boundaries by means of Hermite interpolants. It is interesting to note how the orthogonality along the boundaries, obtained on D, has not been destroyed by the conformal projection. Similarly the specified grid clustering is not affected by the surface topology. The projection maps are based on the isothermal coordinates shown in Figs. 1 and 2.
5. Orthogonal surface grid generation In this section we describe a method for generating orthogonal grids, with fixed boundary points distribution. We will apply it to the mapping between the rectangle D, of aspect ratio M(D),
Fig. 3. Algebraic grid on a paraboloid.
256
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
Fig. 4. Algebraic grid on a canopy-like surface.
and a unit square D ! . Being now D ! our computational domain and D the isothermal parametric domain. It must be pointed out that the approach presented in this section, does not apply only to the particular case of rectangular domains but it is valid for general simply-connected non-degenerate quadrilaterals. Then, more generally, we will describe a grid-generator mapping between the twodimensional simply connected planar domain D, representing the physical non-degenerate quadrilateral, and the computational plane D, usually the unit square. 5.1. Quasiconformal mappings From Eq. (9), we have that an orthogonal coordinate transformation is obtained as a solution of the system ∂x i ∂ 1 ∂x i ∂ F + =0 (15) ∂ξ ∂ξ ∂η F ∂η with {x i } Cartesian coordinates, and i = 1, 2. As pointed out in [8], there is a strong connection between orthogonal mappings and quasiconformal mappings. Introducing the complex planes z = x + iy and ζ = ξ + iη, the mapping functions x(ξ, η) and y(ξ, η), can be interpreted as the real and imaginary part of a complex function z = f (ζ ) = x(ζ ) + iy(ζ ), defined on the computational domain D ⊂ ζ , and with image on the physical domain D ⊂ z. Let dso2 = gij dξ i dξ j be a given metric defined on D, it can be expressed in terms of the complex variables ζ and ζ¯ = ξ − iη, as follows 2 2 (16) dso2 = ρ(ξ, η) dζ + µ(ξ, η) dζ¯ , where
√ g11 + g22 + 2 g , |ρ| = 4 2
(17)
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
and g11 − g22 + 2ig12 µ= √ , g11 + g22 + 2 g
√ g11 + g22 − 2 g |µ| = √ 1. g11 + g22 + 2 g 2
Similarly, on D, for the Cartesian metric ds 2 = dx 2 + dy 2 we have, recalling that z = f (ζ ), fζ¯ 2 2 2 2 ¯ ds = |dz| = |fζ | dζ + dζ , fζ
257
(18)
(19)
with fζ = 12 (fξ − ifη ) and fζ¯ = 12 (fξ + ifη ) being the formal complex derivatives. So it appears that if we require, for fζ = 0, fζ¯ = µfζ ,
(20)
then ds 2 = λ(ζ ) dso2 ; that is the mapping z = f (ζ ), is conformal with respect to a given metric gij . The connection between Eq. (15) and the Beltrami equation (20) follows by showing that Eq. (20) corresponds to the Cauchy–Riemann equations. In fact, let µ = µr + iµi , from (20) we have (1 − µr )xξ = (1 + µr )yη − µi (yξ − xη ), (1 + µr )xη = −(1 − µr )yξ + µi (xξ + yη ), let 1−F , µi = 0, (21) 1+F we have the generalized Cauchy–Riemann equations (from now on referred as GCR), which are the plane analog of Eqs. (7) and (8), µr =
F xξ = yη , xη = −Fyξ .
(22) (23)
Hence, solving system (15) amounts to solve the Beltrami equation for µ real. From the theory of quasiconformal mappings [6,1], every homeomorphic function which is a solution of Eq. (20), with K −1 1, (24) K +1 holding at every point of D, is a K-quasiconformal mapping, and µ is termed the complex dilatation. The connection between orthogonal mapping and K-quasiconformal mappings is more evident considering the geometric definition of K-quasiconformal mappings. Geometrically a quadrilateral D is the conformal image of a rectangle D, if the ratio of the sides of D is equal to the modulus of D. If z = f (ζ ) maps D onto D, with modulus m(D) = M(D), then f is said to be K-quasiconformal if the maximal dilatation K is finite and M(D) 1 K. (25) K m(D) |µ| =
In the case of grid generation, D is usually the unit square, then m(D) = 1, and for non-degenerate quadrilaterals M(D) is always finite. Moreover if µ is real and continuous, as in our case, the mapping z = f (ζ ) is a diffeomorphism, and then it represents a coordinate transformation mapping.
258
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
For quasiconformal mappings, it is possible to state a generalized Riemann mapping theorem [1]. In fact for µ regular, the metric gij defines a structure of Riemann surface on the domain D, and the general uniformization theorem asserts the existence of a conformal mapping. Moreover under suitable boundary conditions, the quasiconformal mapping can be extended to the mapping between the closures D and D. 5.2. Distortion function The previous characterization of quasiconformal mappings as conformal mapping between a
In the frame of grid generation,
with a metric given by gij , and D, implies M(D) = M(D). rectangle D, we are interested in finding a transformation between the unit square D and the quadrilateral D. Then we
→ D is can see that the quasiconformal mapping x i = f (ξ j ) will be of the form f = fˆ ◦ h, where fˆ : D conformal. It follows that the choice of a distortion function F (ξ, η) will amount to specify the mapping
And the relation is given by the expressions (18), (21) between the distortion function F , h : D → D.
the complex dilatation µ and the curvilinear metric gij induced by the mapping h on the rectangle D. The following cases can be analyzed: (a) If the function h is the identity map, gij = δij , then F = 1. We are in the case of a conformal transformation between D and the rectangle D. Of course the two domains must belong to the same conformal equivalence class.
(b) If the function h is represented by the transformation, for (ξˆ , η) ˆ defined on D, ξˆ = M(D)ξ, ηˆ = η; 1 . The transformation f is quasiconformal then we have g11 = M(D)2 , g12 = 0, g22 = 1, and F = M(D) with K = M(D). This case corresponds to the calculation of conformal mappings with unknown conformal module, as described in the previous sections (see also [2]). (c) Another possibility is to represent the function h as composed by the one-dimensional stretching
functions, always with (ξˆ , η) ˆ defined on D,
ξˆ = h1 (ξ ), ηˆ = h2 (η); so g11 = (h1ξ )2 , g12 = 0, and g22 = (h2η )2 . From relations (18) and (21) we have that µ= F =
h1ξ − h2η h1ξ + h2η h2η h1ξ
,
.
This case, combined with the conformal module stretching (case (b)), has been adopted for orthogonal mappings with one-dimensional control [3,4]. (d) It is not possible to specify arbitrary functions hi (ξ, η), because it must be g12 = 0, otherwise µ will not be real. However in [3] it has been shown that it is possible to obtain quasiconformal mappings
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
259
with a more effective control than the one-dimensional stretchings of case (c). This control is obtained specifying along each coordinate line a function governing the arc length distribution between the grid points. This function must be normalized with the condition that the grid points should cover the entire grid line. In this way the functions hi (ξ j ) are more general than simple one-dimensional stretchings, because the distribution functions can vary passing from one to another curve of the same family. Once the function F is specified, system (15) can be solved by adding suitable boundary conditions, in order to respect the requirement of quasiconformal continuation, represented by the GCR equations, and computing the conformal module M(D) by Eq. (11). The numerical algorithm is similar to that described for the calculation of the isothermal coordinates. However in most applications, we are interested in fixing the boundary points. In this case, in order to satisfy the requirement of quasiconformal continuation throughout the boundary, it is not possible to specify the distortion function a priori. Both the distortion function along the boundary and the conformal module must be computed from the boundary distribution of the grid points, in such a way that the GCR are still valid along the boundary. In order to devise a consistent procedure, let introduce the new variables {ξ¯ , η}, ¯ defined by the relations 1 dξ, F dη¯ = dη;
dξ¯ =
then the GCR relations become xξ¯ = yη¯ , xη¯ = −yξ¯ . ¯ let now ω be the differential Defining on D a continuous function g(ζ¯ ) = −y(ζ¯ ) + ix(ζ¯ ) with ζ¯ = ξ¯ + i η, form ω = −y dξ¯ − x dη, ¯ and ω∗ its conjugate, given by ¯ ω∗ = x dξ¯ − y dη. From the GCR relations, it follows that the two differential forms are closed, that is dω = −(xξ¯ − yη¯ ) dξ¯ dη¯ = 0, dω∗ = −(yξ¯ + xη¯ ) dξ¯ dη¯ = 0, then ω is a harmonic form [5]. Moreover for the Stokes’ theorem, dω = dω∗ = 0 on D, implies ω = 0, ω∗ = 0 γ
γ
for any closed curve γ in D. From the previous definitions, it follows that
ω + iω∗ = i(x + iy) dξ¯ + i dη¯ = g(ζ¯ ) dζ¯ .
(26)
260
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
Hence
g ζ¯ dζ¯ =
γ
ω+i
γ
ω∗ = 0,
γ
which implies that g(ζ¯ ) is analytic in the region D. This is a classical result (Morera’s theorem). The complex differential dω + i dω∗ , being zero, is closed, then this form may represent the differential of a function f (ζ¯ ), and from (26) we have df (ζ¯ ) = g(ζ¯ ) dζ¯ .
(27)
Comparing Eqs. (26) and (27), we infer that f (ζ¯ ) = x(ζ¯ ) + iy(ζ¯ ). Being its derivative g(ζ¯ ) analytic, f will be analytic in D. Indeed f (ζ¯ ) represents the conformal mapping between the computational plane D and the domain D. The problem in which we are interested, consists into determining f (ζ¯ ), or equivalently ω(ζ¯ ), from boundary data. From the theorems of existence of Abelian differential forms [5], we have that ω(ζ¯ ) on D is uniquely determined by the condition ω∗ = 0. (28) ∂D
This condition, together with the fact that ω is a closed form, and then on ∂D it must be dω = 0, enable us to compute the conformal module M(D) and the values of F along the boundary, knowing the boundary-point distribution. By Eq. (28), and introducing again the coordinates {ξ, η}, we have that
1 x dξ − y dη = 0, (29) x dξ¯ − y dη¯ = F ∂D
∂D
where ∂D = ∂D1 + ∂D2 (∂D1 with ξ = 0 or 1, ∂D2 with η = 0 or 1). The values of the distortion function F along the boundaries are given by the condition dω = 0 along ∂D, which is satisfied by the GCR relation, which in {ξ, η} coordinates is given by the expression F xξ = yη . Along each side of the computational rectangle D, only yη (on ∂D1 ) or xξ (on ∂D2 ) is given by the boundary distribution. The other comes from the computation of the interior points. Then the boundary values of F must be updated iteratively during the computation, and they must satisfy the integral condition (29). Examples of such a method are given in Figs. 5 and 6, where orthogonal grids have been generated, keeping the same distribution of the boundary points employed in the examples of Figs. 3 and 4. Comparing Figs. 3 and 5, it is evident that the orthogonality in the interior points is more pronounced in the last case (Fig. 5), this because the orthogonality constraint has been imposed in all the domain, and not only along the boundaries.
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
261
Fig. 5. Orthogonal grid on a paraboloid.
Fig. 6. Orthogonal grid on a canopy-like surface.
6. Conclusions In this paper we have presented a novel technique for generating surface grids with specified metric characteristics. The technique consists into two main steps: (i) the parameterization of the surface by isothermal coordinates, (ii) the generation of a two-dimensional grid on the conformal parametric plane and its projection on the surface. In this way, the projection is based on a conformal map, preserving angles and scale-length ratios. It has been proved, by numerical examples, that orthogonal surface grids may be generated on quadrilaterals surface domains, employing well established, and computationally efficient, numerical
262
R. Arina / Applied Numerical Mathematics 46 (2003) 249–262
grid generation methods developed for two-dimensional flat domains. In this context we have described the more general problem of two-dimensional orthogonal grid generation, and it has been shown the possibility of generating orthogonal grids specifying the boundary point distribution along the entire domain boundary.
References [1] L. Ahlfors, Quasiconformal mappings and their applications, in: T.L. Saaty (Ed.), Lecture Notes in Modern Mathematics, Vol. II, Wiley, New York, 1964, pp. 151–164. [2] R. Arina, Orthogonal grids with adaptive control, in: J. Hauser, C. Taylor (Eds.), Numerical Grid Generation in CFD, Pineridge Press Ltd., Swansea, 1986. [3] R. Arina, Adaptive orthogonal curvilinear coordinates, in: K. Morton, M.J. Baines (Eds.), Numerical Methods for Fluid Dynamics III, Oxford Univ. Press, 1988, pp. 353–359. [4] R. Arina, Adaptive orthogonal surface coordinates, in: S. Sengupta, J. Hauser, P.R. Eiseman, J.F. Thompson (Eds.), Numerical Grid Generation in CFD ’88, Pineridge Press Ltd., Swansea, 1988, pp. 351–359. [5] H. Cohn, Conformal Mappings on Riemann Surfaces, Dover, New York, 1980. [6] O. Letho, K.I. Virtanen, Quasiconformal Mappings in the Plane, Springer-Verlag, New York, 1973. [7] C.W. Mastin, J.F. Thompson, Quasiconformal mappings and grid generation, SIAM J. Sci. Statist. Comput. 5 (1984) 52–62. [8] G. Ryskin, L.G. Leal, Orthogonal mapping, J. Comput. Phys. 50 (1983) 71–100. [9] Z.U.A. Warsi, Surface grid generation and differential geometry, in: J.E. Castillo (Ed.), Mathematical Aspects of Numerical Grid Generation, SIAM, Philadelphia, PA, 1991, pp. 99–104.