Differential Equations and Computational Simulations III J. Graef, R. Shivaji, B. Soni & J. Zhu (Editors) Electronic Journal of Differential Equations, Conference 01, 1997, pp. 55-79. ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp 147.26.103.110 or 129.120.3.113 (login: ftp)

Harmonic Parameterization of Geodesic Quadrangles on Surfaces of Constant Curvature

and 2-D Quasi-Isometric Grids ^{∗}

Gennadii A. Chumakov & Sergei G. Chumakov

Abstract

A method for the generation of quasi-isometric boundary-fitted curvilinear coordinates for arbitrary domains is developed on the basis of the quasi-isometric mappings theory and conformal representation of spherical and hyperbolic geometries. A one-parameter family of Riemannian metrics with some attractive invariant properties is analytically described. We construct the quasi-isometric mapping between the regular computation domain R and a given physical domain D that is conformal with respect to the unique metric from the proposed one-parameter class. The identification process of the unknown parameter takes into account the high parametric sensitivity of metrics to the parameter. For this purpose we use a new technique for finding the geodesic quadrangle with given angles and a conformal module on the surface of constant curvature, which makes the method more robust. The method allows more direct control of the grid cells size and angle over the field as the grid is refined. Illustrations of this technique are presented for the case of one-element airfoil and several test domains.

1 Introduction

1.1 Quasi-isometric grids

The generation of a 2-D quasi-isometric grids in a given physical region D (a curvilinear quadrangle with interior angles β_{i}, 0 < β_{i} < π, i = 1, …, 4, and a conformal modulus M) may be considered as a problem of construction of the mapping

X = X(ξ, η), Y = Y (ξ, η)  (1.1)

between points (ξ, η) of the computational region

R = {(ξ, η) : 0 ≤ ξ ≤ 1, 0 ≤ η ≤ 1}

∗1991 Mathematics Subject Classifications: 65N50, 30C30.

Key words and phrases: regular grid generation, quasi-isometric mappings, geodesic grids, geodesic quadrangles, surfaces of constant curvature.

\begin{matrix}c \\ ◯ \end{matrix}1998 Southwest Texas State University and University of North Texas.

Published November 12, 1998.

and points (X, Y) ∈ D such that (1.1) is the unique solution of the following boundary value problem (BVP): given a quasi-isometric mapping between ∂R and ∂D to extend the mapping inside R as a quasi-isometric solution of the appropriate Beltrami system

\sqrt{g_{11}(ξ, η, r)g_{22}(ξ, η, r) − g_{12}^{2}(ξ, η, r)} X_{ξ} = −g_{12}(ξ, η, r) Y_{ξ} + g_{11}(ξ, η, r) Y_{η},

(1.2)

\sqrt{g_{11}(ξ, η, r)g_{22}(ξ, η, r) − g_{12}^{2}(ξ, η, r)} X_{η} = −g_{22}(ξ, η, r) Y_{ξ} + g_{12}(ξ, η, r) Y_{η}.

The coefficients of the system (1.2) contain one unknown parameter r, restricted to the interval (r^{min}, r^{max}), which is to be found in the process of solving the BVP.

In the capacity of the boundary conditions in this BVP we can also choose so-called ”free” conditions, under which grid points on the boundary of the physical region D are not fixed and can move along ∂D.

As it is shown in [3], the linear elliptic system (1.2) can be treated as a condition for conformality of the mapping (1.1) with respect to the following Riemmanian metric defined on R:

ds^{2} = g_{11}(ξ, η, r)dξ^{2} + 2g_{12}(ξ, η, r)dξdη + g_{22}(ξ, η, r)dη^{2}.  (1.3)

Suppose the functions g_{ik}(ξ, η, r) satisfy the inequality

g_{11}(ξ, η, r) + g_{22}(ξ, η, r) ≤ Q \sqrt{g_{11}(ξ, η, r)g_{22}(ξ, η, r) − g_{12}^{2}(ξ, η, r)}, (1.4)

and are continuously differentiable in R. Then the mapping (1.1) has nonvanishing Jacobian in R and is called Q-quasiconformal. That means that the singular values ν_{1} = ν_{1}(ξ, η) and ν_{2} = ν_{2}(ξ, η) of the matrix

J = [ \begin{matrix}X_{ξ} & X_{η} \\ Y_{ξ} & Y_{η} \end{matrix} ] ,

being enumerated in decreasing order, satisfy the condition ν_{1}(ξ, η)/ν_{2}(ξ, η) ≤ Q for all (ξ, η) ∈ R. The least possible number Q with such property is called the coefficient of quasi-conformality of the mapping (1.1) in the domain R. Clearly, a 1-quasiconformal mapping is conformal.

Note that the angle α between the grid lines is defined by the formula cos α = g_{12}(g_{11}g_{22})^{−1/2}, and the ratio of lengths of cell sides is (g_{22}/g_{11})^{1/2}. Thus, since every solution of the Beltrami system (1.2) generates a metric conformally equivalent to (1.3), we are able to control the grid angle and the cell sides ratio by defining the coefficients g_{ik} of the metric (1.3) manually. In order to have better control of the size of cells as the grid is refined, we want to define g_{ik} in such a way that the solution of the corresponding Beltrami system (1.2) is μ-quasi-isometric.

• A mapping (1.1) is called μ-quasi-isometric if there exists a constant μ ∈ [1, ∞) such that

μ = max \{ sup_{(ξ,η)∈R} ν_{1}(ξ, η), sup_{(ξ,η)∈R} \frac{1}{ν_{2}(ξ, η)}\} .

The constant μ is called the coefficient of quasi-isometricity of the mapping (1.1) in the domain R.

Under a μ-quasi-isometric mapping an infinitesimal square will go over into a parallelogram and the ratio of the lengths of any side of the parallelogram and a corresponding side of the square will be bounded by 1/μ and μ. If μ = 1 then (1.1) is an isometric mapping. It is clear that μ-quasi-isometric mapping will also be μ^{2}-quasiconformal. Note that the relation (1.4) and the system (1.2) imply that the mapping (1.1) is quasi-isometric in any sub-region of R. As a rule a quasiconformal mapping is not quasi-isometric, that is, function ν_{1}(ξ, η) or 1/ν_{2}(ξ, η) is not bounded.

Nevertheless use of quasiconformal mappings combined with conditions of existence of bounded derivative of a holomorphic function provides the key to the generation of quasi-isometric coordinate systems.

Several attempts were made in that direction. If we put g_{11} = 1/r, g_{12} = 0, g_{22} = r, we get the quasi-conformal grid generation problem posed by Godunov and Prokopov in [11, 13]. The development of the Godunov-Prokopov method is precisely described in [13] and [24]. In a recent paper Khamayseh and Mastin [17] have extended this method to surface grid generation. For various concepts from tensor analysis and differential geometry applicable to the generation of curvilinear coordinate systems we refer to [25], [9] and [23].

The orthogonal mapping technique has been investigated in recent works [16] and [8], [21] from different points of view.

In [1] Belinsky et. al. proposed to consider the special class of quasi-conformal mappings by defining g_{11} = e^{2q(ξ)}, g_{22} = e^{2p(η)}, g_{12} = \sqrt{g_{11}}\sqrt{g_{22}} cos[β(η)−α(ξ)]. In other words, the proposed metric contained four arbitrary functions of either ξ or η. This metric was obtained as a result of the Chebyshev mapping R onto a plain curvilinear parallelogram. In a similar way in the paper [12] for the purpose of construction of structured multi-block quasi-conformal grids in complex domains it was proposed to use a certain class of functions g_{ij} depending on ξ, η and an unknown vector of parameters r. The metric was obtained by mapping a computational domain that consisted of several squares onto a figure made of several parallelograms.

In order to obtain the unique quasi-isometric solution to the grid generation problem, a special one-parameter family of metrics was closely studied by one of author in papers [5, 6]. In a later paper by Godunov et. al. [10] it was proposed to study a five-parameter family of metrics.

The main goal of the present work is to describe analytically several types of one-parametric families of metrics for which the posed BVP has a unique quasiisometric solution. The metric coefficients g_{ik} will be obtained via mapping the computational region R onto a special class of geodesic quadrangles on surfaces of constant curvature. We start by investigating some basic properties of oneparameter families of such geodesic quadrangles with given angles.

After the basic concepts are introduced, we give several different quasiisometric parameterizations of geodesic quadrangles, that is, quasi-isometric mappings

x = x(ξ, η, α, r), y = y(ξ, η, α, r)  (1.5)

between R and the geodesic quadrangle P. These parameterizations have two vectors of parameters: interior angles α = (α_{1}, …, α_{4}) and so-called Euclidean lengths of sides r = (r_{1}, r_{2}, r_{3}, r_{4}) in the parametric plane (x, y). The plane itself which conformally represents spherical or hyperbolic geometries according to whether the angle defect of P

δ^{∗} = α_{1} + α_{2} + α_{3} + α_{4} − 2π

is positive or negative. The mapping (1.1) generates the metric tensor g_{ij}, coefficients of which are to be used as coefficients of the Beltrami system (1.2).

In [7] the authors proposed several basic parameterizations of a geodesic quadrangle P, and the method for finding the parameter r was described. All parameterizations were obtained by means of using geodesic bundles in order to build vertical and horizontal families of grid lines. In this paper we continue with the subject.

Use of geodesic bundles for parameterizations of P provides a way to construct a mapping (1.1) which maps uniform grid in the unit square on (ξ, η)- plane into a grid in P that is invariant under rotations. By rotation we mean the motion that places another vertex of the geodesic quadrangle at the origin of the parametric plane, and “invariant under rotations” means that rotations and the grid generation process commute. In this case the grid lines appear to be level lines of some harmonic functions, consequently, the geodesic grids in P can be treated as a variant of the Winslow grids [26] with an advantage that in our case the grid in P can be defined explicitly. These kinds of two-dimensional harmonic mappings have non-vanishing Jacobian, as it was shown in [20] and [13]. In particular, we describe a parameterization (1.1) which generates a metric with coefficients g_{ik}(ξ, η, α, r) that have the following properties:

\frac{g_{11}(ξ, η, α, r)}{g_{22}(ξ, η, α, r)} = \frac{g_{22}(1 − η, ξ, σ(α), σ(r))}{g_{11}(1 − η, ξ, σ(α), σ(r))},  (1.6)

\frac{g_{12}(ξ, η, α, r)}{\sqrt{g_{11}(ξ, η, α, r) g_{22}(ξ, η, α, r)}}  (1.7)
= \frac{−g_{12}(1 − η, ξ, σ(α), σ(r))}{\sqrt{g_{11}(1 − η, ξ, σ(α), σ(r)) g_{22}(1 − η, ξ, σ(α), σ(r))}}.

where σ(1, 2, 3, 4) = (4, 1, 2, 3).

From equalities (1.6) and (1.7) it follows that all elements of vectors α and r are equivalent parameters and consequently any five of these parameters may be used as characteristic invariants of P. In particular, we shall fix angles α_{1}, …, α_{4} and study a set of geodesic quadrangles depending on the parameter M for the interval 0 < M < ∞. In capacity of the parameter M we pick the conformal module of P.

We define the conformal module as follows. For every geodesic quadrangle P with vertices z_{i}, i = 1, …, 4, there exists a unique M ∈ (0, ∞) such that there exists a conformal mapping of P onto the rectangle

{(ξ, η) : 0 ≤ ξ ≤ 1, 0 ≤ η ≤ M}

under which the vertex z_{1} is mapped at the origin, and the other vertices of P go over into vertices of the rectangle.

The additional parameters r_{1}, …, r_{4} are then to be determined as monotonic functions r_{i}(M) of the fundamental parameter M. Later we are going to use one of r_{i} as or main parameter for finding the geodesic quadrangle P with given angles α_{1}, …, α_{4} and given a conformal module M. Thus we want to develop the technique for finding such geodesic quadrangle P that takes into account the high parametric sensitivity of M to changes of the parameters r_{i}. For this purpose we shall introduce positive numbers r_{i}^{min}, r_{i}^{max} such that for all possible values of M parameters r_{i} stay within the boundaries r_{i}^{min}, and r_{i}^{max} for i = 1, …, 4. We also will need the derivatives dr_{j}(M)/dr_{i}(M).

In later sections we describe a procedure of finding the mapping f(z), z = x + iy, by which the geodesic quadrangle P is mapped conformally onto the physical domain D; such conformal mapping exists uniquely by virtue of the Riemann Mapping Theorem [15]. It is important to know how the module of the derivative of the conformal mapping behaves on the boundary of P, because it can affect the grid cell size near the boundary of the domain when the grid is refined.

Suppose we fix a point b on a boundary of D. Then by virtue of the uniqueness of the conformal mapping of P onto D there exists a unique proimage z of the point b on boundary of P. Consequently, we allow boundary points of P to move; we call such boundary conditions free on P, and fixed on D. Similarly, if we allow points of D to move along one of the boundaries, and fix boundary points of P on the corresponding boundary, we call such boundary conditions free on D and fixed on P.

In order to construct a quasi-isometric grid we should take into account wellknown conditions of the existence of the boundary derivative (see [19], [18] and [14]):

• If the boundary ∂D is twice continuously differentiable in some neighborhood of the point w_{0} = f(ζ_{0}) and ζ_{0} is not a corner point of P, then the derivative f′ can be continuously extended in certain part of the boundary ∂P , containing ζ_{0}, in such a way that f′(ζ) ≠ 0.

• If ζ_{0} and w_{0} = f(ζ_{0}) are the corner points of P and D respectively and two twice continuously differentiable boundary arcs of D join in the point w_{0} at the same angle as do the corresponding arcs of P in the point ζ_{0}, then f′(ζ) does not equal 0 and is bounded in a certain neighborhood of ζ_{0}.

Thus, in order to use the proposed method for construction of a quasi-isometric grid in D, the following conditions should be satisfied:

1. The boundary ∂D of the domain D must consist of four twice continuously differentiable curves;

2. The sought geodesic quadrangle P is to be chosen in such a way that the angles α_{1}, …, α_{4} of P should coincide with the angles β_{1}, …, β_{4} of the domain D, and the conformal modules of P and D must be the same.

It was proved in [10, 5] and [6], respectively, that under the restriction

δ^{∗} < 2α_{i},  i = 1, …, 4,  (1.8)

a geodesic quadrangle P satisfying the second condition does exist uniquely in both cases of the negative and positive angle defect δ^{∗}.

The sought mapping (1.1) is the superposition of two quasi-isometric mappings mentioned above, and it is to be found as the solution of a variational problem of minimizing a functional of Dirichlet type.

To illustrate the effectiveness of the method, numerical examples of quasiisometric grids with different boundary conditions are presented.

Figure 1a.

Figure 1b.

In the Figure 1a we show a test domain D with a quasi-conformal grid in it, and Figure 1b represents the geodesic quadrangle with same angles and conformal module as of the test domain. If boundary points on D are fixed, then boundary points of P are free, and vice versa.

Note that the Figure 1 illustrates the case when the conditions 1 and 2 are not satisfied, that is, one side of the physical domain is not a C^{2}-curve. We can treat it as a violation of the condition 2. Let the point f(ζ_{0}) split the boundary of D into two C^{2}-curves that intersect at the angle β = b · π, b ∈ [0, 1]. Then for the conformal mapping f : P → D normed such that ω_{0} = f(ζ_{0}) = 0 in neighborhood of ζ_{0} we have

f(z) = a(z − ζ_{0})^{b} + o(|z − ζ_{0}|^{b}), f′(z) = ab(z − ζ_{0})^{b−1} + o(|z − ζ_{0}|^{b−1}),

where a ≠ 0. In particular, it follows that at the point ζ_{0} the derivative f′ becomes infinite, and because of that the grid points next to f(ζ_{0}) will not be close to f(ζ_{0}) as the grid is refined. Examples of quasi-conformal grids with singularities of such kind can be found in the book [23].

In conclusion we would like to say that from our point of view Riemannian metrics that have properties (1.6) and (1.7), and in particular, the metrics induced by the harmonic parameterization, are more efficient for the generation of grids through the solution of the Beltrami system (1.2) with ”free” boundary conditions. We used the harmonic parameterization to generate the distribution of boundary points on the right side of the physical domain in the figure 1.

2 One-parameter families of convex geodesic quadrangles

2.1 Geodesic lines and bundles on surfaces of constant curvature

Consider the surface of constant curvature K = 4δ as the (x, y)- parametric plane with the following metric:

ds^{2} = \frac{dx^{2} + dy^{2}}{[1 + δ(x^{2} + y^{2})]^{2}} ,  (2.1)

where δ is a real number [4]. If δ is positive, the metric (2.1) is defined for every x and y including infinite, and we obtain a representation of spherical geometry. If δ is negative, then in the disk x^{2} + y^{2} + δ < 0 we obtain the Poincare representation of Lobachevsky geometry, or hyperbolic geometry. In the case δ = 0, the metric (2.1) is the Euclidean metric on the plane (x, y).

Let q = ax + by + c[1 − δ(x^{2} + y^{2})]. Then geodesics in the metric (2.1) are curves defined by the equation q = 0. Note that for every δ ≠ 0 circles of the form q = 0 are orthogonal to the circle 1 + δ(x^{2} + y^{2}) = 0 which is called the absolute.

Each of three types of non-Euclidean spaces of a constant curvature indicated above admits the continuous group of isometric mappings, that is, motions. If we denote x + iy by z, then every motion has the form

w(z) = e^{iω} \frac{z − ζ}{1 + δ\overline{ζ}z} ,  (2.2)

where ω ∈ R. The complex number ζ must satisfy the condition |ζ| < |δ|^{−1/2} in case δ < 0.

Let s_{1} and s_{2} be two distinct geodesics defined by equations q_{1} = 0 and q_{2} = 0 respectively. The family F of geodesics orthogonal to s_{1} and s_{2} is called a geodesic bundle. The geodesic bundle F^{⊥} in which s_{1} and s_{2} can be embedded is called orthogonal to F.

2.2 Family of geodesic quadrangles with given angles

Let P be a quadrangle whose sides are geodesics. We assume further that the vertices z_{i} = (x_{i}, y_{i}), i = 1, …, 4 of the quadrangle P are enumerated counterclockwise. Let us denote sides of P as A_{i} = z_{i}z_{i+1}, angles between A_{i−1} and A_{i} as α_{i}, i = 1, …, 4, (assume z_{5} = z_{1} and A_{0} = A_{4}). Let φ_{i} = α_{i} − π/2, i = 1, …, 4.

It is possible to standardize a geodesic quadrangle P by motions (2.2) in such a way that two of its sides on the (x, y)-parametric plane will be linear segments, as shown on the figure 1. In order to do that it is sufficient to move one vertex of P (say, z_{1}) to the origin by an appropriate transformation (2.2). Then by rotation we can place the side A_{1} on the x axis so that the vertex z_{2} has positive x-coordinate. Denote by r_{1} the Euclidean length of the segment A_{1}. The invariance of the metric (2.1) under the motions (2.2) implies that we can associate with every A_{i} its unique Euclidean length r_{i}.

In order to construct a geodesic quadrangle P it is sufficient to have five parameters defined. We shall fix angles α_{1}, …, α_{4} and study one-parameter families of geodesic quadrangles P_{M}, P_{m} or P_{r_{1}}, where in the capacity of varying parameters we can choose the conformal module M, parameter m = (r_{1}r_{3})/(r_{2}r_{4}), or a Euclidean side r_{1}.

From the results presented in [5, 6] it follows that if we consider a oneparametric family P_{M} then invariants m(M) and r_{i}(M) depend monotonically on the conformal module M, which ranges from 0 to ∞. In other words, the following theorem holds.

Theorem 1 Let P_{M} and P\overline{M} be two geodesic quadrangles with the same angles and conformal modules M and \overline{M}, respectively. If M = \overline{M} then P_{M} = P\overline{M}; if M > \overline{M} then m(M) > m(\overline{M}), r_{i}(M) > r_{i}(\overline{M}) and r_{i+1}(M) < r_{i+1}(\overline{M}) for i = 1, 3.

Below we provide some examples of convex geodesic quadrangles P_{M} with same angles and different conformal modules.

a)  b)

c)  d)

Figure 2. Geodesic quadrangle P_{M} with same angles and various conformal module M.

Note that the parameter r_{1}(M) varies in the range from r_{1}^{min} = r_{1}(0) = 0 to r_{1}^{max} = r_{1}(∞) < ∞ as M varies from 0 to ∞ which leads to high parametric sensitivity of conformal module M to r_{1}. When M is close to zero, practically we have a geodesic triangle with a polar coordinate system (Figure 2a ). It is the topological limit of the geodesic quadrangle P_{M} as M → 0. In Figure 2d we see that r_{1}(M) → r_{1}^{max} and r_{2}(M) → 0 as M → ∞.

2.3 Boundaries for Euclidean lengths

Since we are going to consider Euclidean lengths as our main parameters, one of our tasks is to develop a technique for finding the geodesic quadrangle P_{M} by given conformal module M, which takes into account the high parametric sensitivity of M to r_{i}.

First of all for this purpose we shall find positive numbers r_{i}^{min} and r_{i}^{max} such that r_{i}(M) belongs to the interval (r_{i}^{min}, r_{i}^{max}) for all M, and r_{i}^{min} = r_{i}(0), r_{i}^{max} = r_{i}(∞) for i = 1, …, 4. As it was shown in ([7]), we can calculate r_{i}^{min} and r_{i}^{max} using the special function B(ψ_{1}, ψ_{2}, ψ_{3}, ψ_{4}, δ) given by

B(ψ_{1}, ψ_{2}, ψ_{3}, δ) = \sqrt{\frac{sin ψ cos(ψ − ψ_{3})}{δ cos(ψ − ψ_{2}) sin(ψ − ψ_{2} − ψ_{3})}},  (2.3)

where ψ = (ψ_{1} + ψ_{2} + ψ_{3} + π/2)/2. It is convenient for us to fix δ = sin γ, where γ = (φ_{1} + φ_{2} + φ_{3} + φ_{4})/2.

Let us denote by l(P) the number of sides of P the sum of whose adjacent angles is not less then π. The value of l(P) can be 0, 1, 2, 3 and 4. Let us consider each case separately. Let Σ_{4} be the cyclic subgroup of the permutation group S_{4}, generated by the element

σ̄ ∈ S_{4},  σ̄ = ( \begin{matrix}1 & 2 & 3 & 4 \\ 2 & 3 & 4 & 1 \end{matrix} ) ,

Let l(P) = 0, then for all σ ∈ Σ_{4}

r_{σ(1)}^{min} = B(φ_{σ(1)}, φ_{σ(2)}, −π/2, δ),  r_{σ(1)}^{max} = B(φ_{σ(1)}, π/2, φ_{σ(4)}, δ),
r_{σ(2)}^{min} = B(φ_{σ(2)}, φ_{σ(3)}, −π/2, δ),  r_{σ(2)}^{max} = B(φ_{σ(2)}, π/2, φ_{σ(1)}, δ),
r_{σ(3)}^{min} = B(φ_{σ(3)}, φ_{σ(4)}, −π/2, δ),  r_{σ(3)}^{max} = B(φ_{σ(3)}, π/2, φ_{σ(2)}, δ),
r_{σ(4)}^{min} = B(φ_{σ(4)}, φ_{σ(1)}, −π/2, δ),  r_{σ(4)}^{max} = B(φ_{σ(4)}, π/2, φ_{σ(3)}, δ),

Let l(P) = 1, and σ ∈ Σ_{4} be such that φ_{σ(1)} + φ_{σ(2)} ≥ 0. Then the boundaries will be as follows.

r_{σ(1)}^{min} = 0,
r_{σ(1)}^{max} = B(φ_{σ(1)}, π/2, φ_{σ(4)}, δ),
r_{σ(2)}^{min} = B(φ_{σ(3)}, φ_{σ(2)}, −π/2, δ),
r_{σ(2)}^{max} = B(φ_{σ(1)} + φ_{σ(2)} − π/2, φ_{σ(3)}, φ_{σ(4)}, δ),
r_{σ(3)}^{min} = B(φ_{σ(3)}, φ_{σ(4)}, φ_{σ(1)} + φ_{σ(2)} − π/2, δ),
r_{σ(3)}^{max} = B(φ_{σ(3)}, −π/2, φ_{σ(2)}, δ),
r_{σ(4)}^{min} = B(φ_{σ(1)}, φ_{σ(4)}, −π/2, δ),
r_{σ(4)}^{max} = B(φ_{σ(1)} + φ_{σ(2)} − π/2, φ_{σ(4)}, φ_{σ(3)}, δ).

Consider now the case l(P) ≥ 2. There always exist such σ ∈ Σ_{4} that

φ_{σ(1)} + φ_{σ(2)} ≥ φ_{σ(3)} + φ_{σ(4)}, φ_{σ(1)} + φ_{σ(4)} ≥ φ_{σ(2)} + φ_{σ(3)}

holds. Then the boundaries for r_{j} will be as follows:

r_{σ(1)}^{min} = 0,
r_{σ(1)}^{max} = B(φ_{σ(1)} + φ_{σ(4)} − π/2, φ_{σ(2)}, φ_{σ(3)}, δ),
r_{σ(2)}^{min} = B(φ_{σ(2)}, φ_{σ(3)}, φ_{σ(1)} + φ_{σ(4)} − π/2, δ),
r_{σ(2)}^{max} = B(φ_{σ(1)} + φ_{σ(2)} − π/2, φ_{σ(3)}, φ_{σ(4)}, δ),
r_{σ(3)}^{min} = B(φ_{σ(3)}, φ_{σ(4)}, φ_{σ(1)} + φ_{σ(2)} − π/2, δ),
r_{σ(3)}^{max} = B(φ_{σ(1)} + φ_{σ(4)} − π/2, φ_{σ(3)}, φ_{σ(2)}, δ),
r_{σ(4)}^{min} = 0,
r_{σ(4)}^{max} = B(φ_{σ(1)} + φ_{σ(2)} − π/2, φ_{σ(4)}, φ_{σ(3)}, δ).

2.4 Relations between Euclidean lengths

The relations r_{σ(4)}(r_{σ(1)}) and dr_{σ(4)}/d(r_{σ(1)}) for all σ ∈ Σ_{4} can be obtained from the following equation which involves parameters r_{σ(1)}, r_{σ(4)}, φ_{1}, φ_{2}, φ_{3}, φ_{4} and δ:

C_{σ(3)}S_{0} = δC_{σ(2)}S_{σ(23)}r_{σ(1)}^{2} + δC_{σ(4)}S_{σ(34)}r_{σ(4)}^{2}  (2.4)
+2δD_{σ}r_{σ(1)}r_{σ(4)} − δ^{2}C_{σ(1)}S_{σ(24)}r_{σ(1)}^{2}r_{σ(4)}^{2},

where

S_{0} = sin γ, D_{σ} = cos φ_{σ(2)} cos φ_{σ(4)}, C_{i} = cos(γ − φ_{i}),
S_{ij} = sin(γ − φ_{i} − φ_{j}), σ(ij) = (σ(i), σ(j)) , i, j = 1, …, 4.

It follows from (2.4) that for δ ≥ 0:

r_{σ(4)}(r_{σ(1)}) = \frac{1}{B + \sqrt{B^{2} + C}} ,  (2.5)

where

B = \frac{r_{σ(1)}D_{σ}}{C_{σ(3)} − r_{σ(1)}^{2}C_{σ(2)}S_{σ(23)}} , C = \frac{r_{σ(1)}^{2}C_{σ(1)}S_{σ(24)}δ − C_{σ(4)}S_{σ(34)}}{C_{σ(3)} − r_{σ(1)}^{2}C_{σ(2)}S_{σ(23)}} .

In the case δ < 0 we have to take as r_{σ(4)} the root of (2.4) that belongs to the interval (r_{σ(4)}^{min}, r_{σ(4)}^{max}).

Thus given r_{σ(1)} and four angles α_{1}, …, α_{4} we are able by means of (2.4) to determine the function r_{j} = r_{j}(r_{i}), and its derivative with respect to r_{i} for i, j = 1, …, 4. In order to do that it is sufficient to find the derivative of r_{σ(4)}(r_{σ(1)}) with respect to r_{σ(1)}:

\frac{dr_{σ(4)}(r_{σ(1)})}{dr_{σ(1)}} = −\frac{C_{σ(2)}S_{σ(23)}r_{σ(1)} + D_{σ}r_{σ(4)} − δC_{σ(1)}S_{σ(24)}r_{σ(1)}r_{σ(4)}^{2}}{C_{σ(4)}S_{σ(34)}r_{σ(4)} + D_{σ}r_{σ(1)} − δC_{σ(1)}S_{σ(24)}r_{σ(1)}^{2}r_{σ(4)}} . (2.6)

3 Quasi–isometric parameterizations of geodesic quadrangles

It should be noted that any five parameters from vectors α = (α_{1}, …, α_{4}) and r = (r_{1}, r_{2}, r_{3}, r_{4}) may be used as characteristic invariants of a geodesic quadrangle P. We shall now suppose that parameters φ_{1}, φ_{2}, φ_{4}, r_{1} and r_{4} are given.

3.1 The simplest quasi–isometric parameterization

Consider P as an intersection of two geodesic bundles, which are described by the following algebraic equations

(cos φ_{1} + ξa_{1})x + (sin φ_{1} − ξb_{1})y − ξc_{1}[1 − δ(x^{2} + y^{2})] = 0, ξ ∈ [0, 1], (3.1)

ηa_{2}x − (1 + ηb_{2})y + ηc_{2}[1 − δ(x^{2} + y^{2})] = 0,  η ∈ [0, 1],  (3.2)

where

a_{1} = cos φ_{2}(1 − δr_{1}^{2}) − cos φ_{1}, b_{1} = sin φ_{2}(1 + δr_{1}^{2}) + sin φ_{1},
c_{1} = r_{1} cos φ_{2}, a_{2} = sin(φ_{1} + φ_{4}) − δr_{4}^{2} sin(φ_{1} − φ_{4}),
b_{2} = cos(φ_{1} + φ_{4}) − δr_{4}^{2} cos(φ_{1} − φ_{4}) − 1, c_{2} = r_{4} cos φ_{4}.

Let us denote by

P (ξ) = c_{2} cos φ_{1} + ξ(a_{1}c_{2} + a_{2}c_{1}),
Q(ξ, η) = ξc_{1} − ηc_{2} sin φ_{1} + ξη(b_{1}c_{2} + b_{2}c_{1}),
F (ξ, η) = \frac{1}{2} [cos φ_{1} + ξa_{1} + η(a_{2} sin φ_{1} + b_{2} cos φ_{1}) + ξη(a_{1}b_{2} − a_{2}b_{1})] , (3.3)
2W (ξ, η) = F (ξ, η) + \sqrt{F (ξ, η)^{2} + δ(Q(ξ, η)^{2} + η^{2}P (η)^{2})},

then by virtue of the relations (3.1)-(3.2) we can obtain the following mapping of the unit square R onto the geodesic quadrangle P:

x(ξ, η) = \frac{Q(ξ, η)}{2W (ξ, η)},  y(ξ, η) = \frac{η · P (ξ)}{2W (ξ, η)}.  (3.4)

Note that if the intersection of geodesic bundles (3.2) and (3.1) is a geodesic convex quadrangle then for all ξ ∈ [0, 1] the inequality P (ξ) > 0 holds and the mapping (3.4) is a quasi-isometric one.

From (3.4) we have by differentiation

x_{ξ} = A[(1 + ηb_{2})W + δη^{2}c_{2}P]z, x_{η} = P [(b_{1}ξ − sin φ_{1})W − δξηc_{1}P]z,

y_{ξ} = η · A(a_{2}W − δc_{2}Q)z, y_{η} = P [(cos φ_{1} + ξa_{1})W + δξc_{1}Q]z, (3.5)

where P = P (ξ), Q = Q(ξ, η), F = F (ξ, η) and W = W (ξ, η) were defined by (3.3) and

A = c_{1} cos φ_{1} + η [(a_{1}c_{2} + a_{2}c_{1}) sin φ_{1} + (b_{1}c_{2} + b_{2}c_{1}) cos φ_{1}] ,

4z = W^{−2}[F^{2} + δ(Q^{2} + η^{2}P^{2})]^{−1/2}.

We should mention that the mapping (3.4) generates a class of conformally equivalent Riemannian metrics. Since any two conformally equivalent metrics lead to the same Beltrami system, we can take in the capacity of a representative the metric with following coefficients

g_{11} = A^{2}{[(1 + ηb_{2})W + δη^{2}c_{2}P]^{2} + η^{2}(a_{2}W − δc_{2}Q)^{2}},  (3.6)

g_{22} = P^{2}{[(b_{1}ξ − sin φ_{1})W − δξηc_{1}P]^{2} + [(cos φ_{1} + ξa_{1})W + δξc_{1}Q]^{2}},

g_{12} = AP {[(1 + ηb_{2})W + δη^{2}c_{2}P][(b_{1}ξ − sin φ_{1})W − δξηc_{1}P]
+η^{2}(a_{2}W − δc_{2}Q)[(cos φ_{1} + ξa_{1})W + δξc_{1}Q]}.

Below we provide the examples where we used the simplest quasi-isometric parameterization in order to construct a grid inside four geodesic quadrangles on the surface of the positive constant curvature with the same angles and conformal module M, but with different enumerations of the sides and angles.

a)  b)

c)  d)

Figure 3. The simplest quasi-isometric parameterization of P

3.2 “Even” quasi-isometric parameterization

Note that in the transformation (3.4) we can take in the capacity of ξ and η two arbitrary monotonically increasing functions ξ(ξ_{1}) and η(η_{1}) satisfying the conditions ξ(0) = η(0) = 0 and ξ(1) = η(1) = 1. In particular, we can choose ξ = ξ(ξ_{1}) and η = η(η_{1}) in such a way that under the mapping (3.4) the uniform distribution of points on the lower and the left sides of the unit square holds, i.e., the distribution of points on sides A_{1} and A_{4} of the geodesic quadrangle P

a)  b)

c)  d)

Figure 4. The “even” quasi-isometric parameterization.

is also uniform in a sense of the Euclidean distance. In order to obtain such a mapping it is sufficient to choose ξ(ξ_{1}) and η(η_{1}) as follows:

ξ(ξ_{1}) = \frac{ξ_{1}r_{1} cos φ_{1}}{c_{1}(1−δξ_{1}^{2}r_{1}^{2})−a_{1}ξ_{1}r_{1}} ,
η(η_{1}) = \frac{η_{1}r_{4} cos φ_{1}}{c_{2}(1−δη_{1}^{2}r_{4}^{2})−η_{1}r_{4}(a_{1} sin φ_{1}+b_{2} cos φ_{1})} .

Examples of curvilinear coordinate systems in the geodesic quadrangle on a surface of the constant negative curvature generated by means of “even” quasiisometric parameterization are given below.

Figure 4 represents four geodesic quadrangles with the same angles and the same ”Euclidean” lengths of sides and, consequently, the same conformal module M. In the Figure 4a is shown the geodesic quadrangle with parameters (φ_{1}, φ_{2}, φ_{3}, φ_{4}) = (0.5, −1.0, −0.8, −1.0).

Note that the grids in the geodesic quadrangle P constructed by the “even” and the simplest parameterizations are not invariant under the cyclic permutation vertices of P in the origin by means of the motion group. If one wants to obtain the grid that is invariant under the motion group, one can look below in the next two sections where we describe the harmonic parameterization.

3.3 The harmonic parameterization of P on the Euclidean plane

We now obtain parameterizations which are based upon geodesic bundles and have the invariant properties (1.6)-(1.7). We first consider the case φ_{1} + φ_{2} + φ_{3} + φ_{4} = 0, that is, δ = 0, so that geodesics in the metric (2.1) are straight lines defined by the equation ax + by + c = 0. Consider P as an intersection of two bundles, which are described by the following equations:

x cos[φ_{1} − ξ(φ_{1} + φ_{2})] + y sin[φ_{1} − ξ(φ_{1} + φ_{2})] − r_{1} cos(φ_{2})\frac{sin[ξ(φ_{1}+φ_{2})]}{sin(φ_{1}+φ_{2})} = 0

0 ≤ ξ ≤ 1,  (3.7)
x sin[η(φ_{1} + φ_{4})] − y cos[η(φ_{1} + φ_{4})] + r_{4} cos φ_{4} \frac{sin[η(φ_{1}+φ_{4})]}{sin(φ_{1}+φ_{4})} = 0
0 ≤ η ≤ 1  (3.8)

In order to solve the system (3.7)–(3.8) we denote by

u(ξ, η) = \frac{r_{1} cos φ_{2}}{cos[φ_{1} − ξ(φ_{1} + φ_{2}) − η(φ_{1} + φ_{4})]} \frac{sin[ξ(φ_{1} + φ_{2})]}{sin(φ_{1} + φ_{2})},  (3.9)

v(ξ, η) = \frac{r_{4} cos φ_{4}}{cos[φ_{1} − ξ(φ_{1} + φ_{2}) − η(φ_{1} + φ_{4})]} \frac{sin[η(φ_{1} + φ_{4})]}{sin(φ_{1} + φ_{4})},  (3.10)

then

x(ξ, η) = u(ξ, η) cos[η(φ_{1} + φ_{4})] − v(ξ, η) sin[φ_{1} − ξ(φ_{1} + φ_{2})],  (3.11)
y(ξ, η) = u(ξ, η) sin[η(φ_{1} + φ_{4})] + v(ξ, η) cos[φ_{1} − ξ(φ_{1} + φ_{2})],  (3.12)

is the solution.

From (3.11)–(3.12) we have by differentiation

x_{ξ} = \frac{V (η) cos[η(φ_{1} + φ_{4})]}{cos^{2}[φ_{1} − ξ(φ_{1} + φ_{2}) − η(φ_{1} + φ_{4})]},
x_{η} = −\frac{U(ξ) sin[φ_{1} − ξ(φ_{1} + φ_{2})]}{cos^{2}[φ_{1} − ξ(φ_{1} + φ_{2}) − η(φ_{1} + φ_{4})]},
y_{ξ} = −\frac{V (η) sin[η(φ_{1} + φ_{4})]}{cos^{2}[φ_{1} − ξ(φ_{1} + φ_{2}) − η(φ_{1} + φ_{4})]},
y_{η} = \frac{U(ξ) cos[φ_{1} − ξ(φ_{1} + φ_{2})]}{cos^{2}[φ_{1} − ξ(φ_{1} + φ_{2}) − η(φ_{1} + φ_{4})]},

where

U(ξ) = r_{1}(φ_{1} + φ_{4}) cos φ_{2} \frac{sin[ξ(φ_{1} + φ_{2})]}{sin(φ_{1} + φ_{2})}
+r_{4} cos φ_{4} cos[φ_{1} − ξ(φ_{1} + φ_{2})]\frac{φ_{1} + φ_{4}}{sin(φ_{1} + φ_{4})},  (3.13)
V (η) = r_{1} cos φ_{2} cos[φ_{1} − η(φ_{1} + φ_{4})]\frac{φ_{1} + φ_{2}}{sin(φ_{1} + φ_{2})}
+r_{4}(φ_{1} + φ_{2}) cos φ_{4} \frac{sin[η(φ_{1} + φ_{4})]}{sin(φ_{1} + φ_{4})}.  (3.14)

It can happen that one obtains an indeterminate expression of the form sin[ξ(φ_{1} + φ_{2})]/ sin(φ_{1} + φ_{2}) as φ_{1} + φ → 0 which is evidently equal to ξ.

It is easily seen that for all (ξ, η) ∈ R the inequalities U(ξ) > 0 and V (η) > 0 hold and the mapping (3.11)–(3.12) is a quasi-isometric transformation R onto P. This mapping generates a class conformally equivalent Riemannian metrics. In the capacity of a representative of the class of metrics we can take the metric with the following coefficients:

g_{11} = V^{2}(η), g_{22} = U^{2}(ξ),  (3.15)
g_{12} = U(ξ)V (η) sin[−φ_{1} + ξ(φ_{1} + φ_{2}) + η(φ_{1} + φ_{4})].

Let us have a uniform grid in ξ, η− plane, where n, m− are the total numbers of grid points to be used in ξ, η directions. Then the harmonic parameterization has a property that under the cyclic permutation (m, n) → (n, m), (φ_{1}, φ_{2}, φ_{3}, φ_{4}) → (φ_{4}, φ_{1}, φ_{2}, φ_{3}), (r_{1}, r_{2}, r_{3}, r_{4}) → (r_{4}, r_{1}, r_{2}, r_{3}), the transformation (3.11)–(3.12) results in the same grid in P.

The metric (3.15) by virtue of its simplicity can be used for the generation of quasi-conformal grids in the physical domain D with the angles satisfying the inequality α_{1} + α_{2} + α_{3} + α_{4} ≠ 2π as well. For this purpose it is sufficient to substitute in (3.13)—(3.15)

φ_{i} = α_{i} − \frac{α_{1} + α_{2} + α_{3} + α_{4}}{4} , i = 1, …, 4.

3.4 The harmonic parameterization of P on a surface of constant curvature

Let P be a geodesic quadrangle on a surface of the constant curvature K = 4δ, where δ = sin[(φ_{1} + φ_{2} + φ_{3} + φ_{4})/2]. Let us define for i = 1, 2 the angles w_{i} of the intersection of sides A_{i} and A_{i+2} of the geodesic quadrangle P, constants Δ_{i} and functions SH(ξw_{i}), CH(ξw_{i}):

ω_{i} = \{ \begin{matrix}arccot(d_{i}/b_{i}), if l_{i} ≥ 0, \\ coth^{−1}(d_{i}/b_{i}), if l_{i} < 0, \end{matrix} Δ_{i} = \{ \begin{matrix}+1, if l_{i} ≥ 0, \\ −1, if l_{i} < 0, \end{matrix}

SH(ξω_{i}) = \{ \begin{matrix}sin(ξω_{i}), if l_{i} ≥ 0, \\ sinh(ξω_{i}), if l_{i} < 0, \end{matrix}

CH(ξω_{i}) = \{ \begin{matrix}cos(ξω_{i}), if l_{i} ≥ 0, \\ cosh(ξω_{i}), if l_{i} < 0 \end{matrix}  (3.16)

where l_{i} = a_{i}^{2} + 4δc_{i}^{2} and

a_{1} = − sin(φ_{1} + φ_{2}) + δr_{1}^{2} sin(φ_{1} − φ_{2}),
a_{2} = sin(φ_{1} + φ_{4}) − δr_{4}^{2} sin(φ_{1} − φ_{4}),
b_{1} = \sqrt{|l_{1}|} sign δ, b_{2} = \sqrt{|l_{2}|} sign δ,
c_{1} = r_{1} cos φ_{2}, c_{2} = r_{4} cos φ_{4},
d_{1} = cos(φ_{1} + φ_{2}) − δr_{1}^{2} cos(φ_{1} − φ_{2}),
d_{2} = cos(φ_{1} + φ_{4}) − δr_{4}^{2} cos(φ_{1} − φ_{4}).

Here we suppose that sign δ = 1 if δ = 0.

Consider P as an intersection of two geodesic bundles

C(b_{1}, a_{1}, ξω_{1})x − S(b_{1}, a_{1}, ξω_{1})y + c_{1}SH(ξω_{1})[1 − δ(x^{2} + y^{2})] = 0,

0 ≤ ξ ≤ 1,  (3.17)
a_{2}SH(ηω_{2})x − b_{2}CH(ηω_{2})y + c_{2}SH(ηω_{2})[1 − δ(x^{2} + y^{2})] = 0,
0 ≤ η ≤ 1,  (3.18)

where

C(b_{1}, a_{1}, ξω_{1}) = −b_{1} cos φ_{1}CH(ξω_{1}) + a_{1} sin φ_{1}SH(ξω_{1}),  (3.19)
S(b_{1}, a_{1}, ξω_{1}) = b_{1} sin φ_{1}CH(ξω_{1}) + a_{1} cos φ_{1}SH(ξω_{1}).  (3.20)

Introducing functions

C̄(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}) = b_{2}CH(ηw_{2})C(b_{1}, a_{1}, ξw_{1})
−a_{2}SH(ηw_{2})S(b_{1}, a_{1}, ξw_{1}),  (3.21)
S̄(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}) = b_{2}CH(ηw_{2})S(b_{1}, a_{1}, ξw_{1})
+a_{2}SH(ηw_{2})C(b_{1}, a_{1}, ξw_{1}). (3.22)

we can write the desired parameterization of P in the form

x(ξ, η)  (3.23)
= \frac{2c_{2}SH(ηw_{2})S(b_{1}, a_{1}, ξw_{1}) − 2c_{1}b_{2}SH(ξw_{1})CH(ηw_{2})}{C̄(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}) + Δ_{0}\sqrt{C̄^{2}(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}) + 4δh(ξ, η)}},
y(ξ, η)  (3.24)
= \frac{2SH(ηw_{2})[c_{2}C(b_{1}, a_{1}, ξw_{1}) − c_{1}a_{2}SH(ξw_{1})]}{C̄(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}) + Δ_{0}\sqrt{C̄^{2}(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}) + 4δh(ξ, η)}},

where we have set Δ_{0} = sign {SH(ηw_{2})[c_{2}C(b_{1}, a_{1}, ξw_{1}) − c_{1}a_{2}SH(ξw_{1})]} ,

h(ξ, η) = c_{1}^{2}SH^{2}(ξw_{1})[a_{2}^{2}SH^{2}(ηw_{2}) + b_{2}^{2}CH^{2}(ηw_{2})]
+c_{2}^{2}SH^{2}(ηw_{2})[a_{1}^{2}SH^{2}(ξw_{1}) + b_{1}^{2}CH^{2}(ξw_{1})]
−2c_{1}c_{2}SH(ξw_{1})SH(ηw_{2})S̄(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}).

Figure 5. An example of “harmonic” parameterization.

72 Harmonic Parameterization of Geodesic Quadrangles

3.5 Riemannian metric of the harmonic parameterization

of P

The calculation of x_{ξ} and y_{ξ} is based on the relationships:

\frac{∂}{∂ξ} C(b_{1}, a_{1}, ξw_{1}) = w_{1}S(a_{1}, b_{1}, Δ_{1}ξw_{1}),  (3.25)
\frac{∂}{∂ξ} S(b_{1}, a_{1}, ξw_{1}) = −w_{1}C(a_{1}, b_{1}, Δ_{1}ξw_{1}),
\frac{∂}{∂ξ} C̄(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}) = w_{1}S̄(b_{2}, a_{2}, a_{1}, b_{1}, Δ_{1}ξw_{1}, ηw_{2}),
\frac{∂}{∂ξ} S̄(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}) = −w_{1}C̄(b_{2}, a_{2}, a_{1}, b_{1}, Δ_{1}ξw_{1}, ηw_{2}),
\frac{∂}{∂η} C̄(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}) = −w_{2}S̄(a_{2}, b_{2}, b_{1}, a_{1}, ξw_{1}, Δ_{2}ηw_{2}). (3.26)

We use the abbreviations

C̄ = C̄(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}), S̄ = S̄(b_{2}, a_{2}, b_{1}, a_{1}, ξw_{1}, ηw_{2}),

h = h(ξ, η), C̄_{1} = C̄(b_{2}, a_{2}, a_{1}, b_{1}, Δ_{1}ξw_{1}, ηw_{2}),

C̄_{2} = C̄(a_{2}, b_{2}, b_{1}, a_{1}, ξw_{1}, Δ_{2}ηw_{2}), z_{0} = Δ_{0}\sqrt{C̄^{2} + 4δh},

S̄_{1} = S̄(b_{2}, a_{2}, a_{1}, b_{1}, Δ_{1}ξw_{1}, ηw_{2}), z_{1} = C̄ + z_{0} ,
S̄_{2} = S̄(a_{2}, b_{2}, b_{1}, a_{1}, ξw_{1}, Δ_{2}ηw_{2}),

and find

\frac{∂}{∂ξ} [C(b_{1}, a_{1}, ξw_{1})/z_{1}]
= w_{1} [S(a_{1}, b_{1}, Δ_{1}ξw_{1})C̄ − C(b_{1}, a_{1}, ξw_{1})S̄] /(z_{1}z_{0})
+2δ [2w_{1}hS(a_{1}, b_{1}, Δ_{1}ξw_{1}) − h_{ξ}C̄(b_{1}, a_{1}, ξw_{1})] /(z_{1}^{2}z_{0}),
\frac{∂}{∂ξ} [SH(ξw_{1})/z_{1}] = w_{1} [CH(ξw_{1})C̄ − SH(ξw_{1})S̄_{1}] /(z_{1}z_{0})
+2δ [2w_{1}hCH(ξw_{1}) − h_{ξ}SH(ξw_{1})] /(z_{1}^{2}z_{0}),
\frac{∂}{∂ξ} [S(b_{1}, a_{1}, ξw_{1})/z_{1}]
= −w_{1} [C(a_{1}, b_{1}, Δ_{1}ξw_{1})C̄ + S(b_{1}, a_{1}, ξw_{1})S̄_{1}] /(z_{1}z_{0})
−2δ [2w_{1}hC(a_{1}, b_{1}, Δ_{1}ξw_{1}) + h_{ξ}S(b_{1}, a_{1}, ξw_{1})] /(z_{1}^{2}z_{0}).

It is easily seen that

S(a_{1}, b_{1}, Δ_{1}ξw_{1})C̄ − C(b_{1}, a_{1}, ξw_{1})S̄ = −a_{1}b_{1}a_{2}SH(ηw_{1}),

CH(ξw_{1})C̄ − SH(ξw_{1})S̄_{1} = b_{1}C(b_{2}, a_{2}, −ηw_{2}),

C(a_{1}, b_{1}, Δ_{1}ξw_{1})C̄ + S(b_{1}, a_{1}, ξw_{1})S̄_{1} = a_{1}b_{1}b_{2}CH(ηw_{2}),

Gennadii A. Chumakov & Sergei G. Chumakov 73

and we therefore have

x_{ξ}  (3.27)
= −2b_{1}b_{2}w_{1}CH(ηw_{2}) [c_{1}C(b_{2}, a_{2}, −ηw_{2}) + c_{2}a_{1}SH(ηw_{2})] /(z_{1}z_{0})

+4δ{ 2w_{1}h[c_{1}b_{2}CH(ξw_{1})CH(ηw_{2}) + c_{2}C(a_{1}, b_{1}, Δ_{1}ξw_{1})SH(ηw_{2})] +h_{ξ}[−c_{1}b_{2}SH(ξw_{1})CH(ηw_{2}) + c_{2}S(b_{1}, a_{1}, ξw_{1})SH(ηw_{2})]}/(z_{1}^{2}z_{0}),

y_{ξ} = −2a_{2}b_{1}w_{1}SH(ηw_{2}) [c_{1}C(b_{2}, a_{2}, −ηw_{2}) + c_{2}a_{1}SH(ηw_{2})] /(z_{1}z_{0})
−4δSH(ηw_{2}) { 2w_{1}h[c_{1}a_{2}CH(ξw_{1}) − c_{2}S(a_{1}, b_{1}, Δ_{1}ξw_{1})]− (3.28)
−h_{ξ}[c_{1}a_{2}SH(ξw_{1}) − c_{2}C(b_{1}, a_{1}, ξw_{1})]} /(z_{1}^{2}z_{0}),

where

h_{ξ} = 2w_{1}{c_{1}^{2}SH(ξw_{1})CH(ξw_{1})[b_{2}^{2}CH^{2}(ηw_{2}) + a_{2}^{2}SH(ηw_{2})]
+c_{2}^{2}SH^{2}(ηw_{2})SH(ξw_{1})CH(ξw_{1})[a_{1}^{2} − Δ_{1}b_{1}^{2}]
−c_{1}c_{2}SH(ηw_{2})[CH(ξw_{1})S̄_{1} − SH(ξw_{1})C̄_{1}]}.

Next we calculate x_{η} and y_{η}. Using (3.26) we find

\frac{∂}{∂η} [SH(ηw_{2})/z_{1}] = w_{2} [CH(ηw_{2})C̄ + SH(ηw_{2})S̄_{2}] /(z_{1}z_{0})
+2δ [2w_{2}hCH(ηw_{2}) − h_{η}SH(ηw_{2})] /(z_{1}^{2}z_{0}),
\frac{∂}{∂η} [CH(ηw_{2})/z_{1}] = −w_{2} [Δ_{2}SH(ηw_{2})C̄ − CH(ηw_{2})S̄_{2}] /(z_{1}z_{0})
−2δ [2Δ_{2}w_{2}hSH(ηw_{2}) + h_{η}CH(ηw_{2})] /(z_{1}^{2}z_{0}).

It is also clear that

CH(ηw_{2})C̄ + SH(ηw_{2})S̄_{2} = b_{2}C(b_{1}, a_{1}, ξw_{1}),

Δ_{2}SH(ηw_{2})C̄ − CH(ηw_{2})S̄_{2} = −a_{2}S(b_{1}, a_{1}, ξw_{1}),

and we therefore have

x_{η} = 2b_{2}w_{2}S(b_{1}, a_{1}, ξw_{1})[−c_{1}a_{2}SH(ξw_{1}) + c_{2}C(b_{1}, a_{1}, ξw_{1})]/(z_{1}z_{0}) (3.29)

+4δ{2w_{2}h[Δ_{2}c_{1}b_{2}SH(ξw_{1})SH(ηw_{2}) + c_{2}S(b_{1}, a_{1}, ξw_{1})CH(ηw_{2})] +h_{η}[c_{1}b_{2}SH(ξw_{1})CH(ηw_{2}) − c_{2}S(b_{1}, a_{1}, ξw_{1})SH(ηw_{2})]}/(z_{1}^{2}z_{0}),

y_{η} = [c_{2}C(b_{1}, a_{1}, ξw_{1}) − c_{1}a_{2}SH(ξw_{1})]{2b_{2}w_{2}C(b_{1}, a_{1}, ξw_{1})/(z_{1}z_{0})
+4δ[2w_{2}hCH(ηw_{2}) − h_{η}SH(ηw_{2})]}/(z_{1}^{2}z_{0}),  (3.30)

where

h_{η} = 2w_{2}{c_{1}^{2}SH^{2}(ξw_{1})SH(ηw_{2})CH(ηw_{2})[a_{2}^{2} − Δ_{2}b_{2}^{2}]
+c_{2}^{2}SH^{2}(ηw_{2})CH(ηw_{2})[b_{1}^{2}CH^{2}(ξw_{1}) + a_{1}^{2}SH^{2}(ξw_{1})]
−c_{1}c_{2}SH(ξw_{1})[CH(ηw_{2})S̄ + SH(ηw_{2})C̄_{2}]}.

Thus, introducing functions (3.16)–(3.20) and (3.21), (3.22) we have found x_{ξ}, x_{η}, y_{ξ}, y_{η} and therefore coefficients of the metric tensor

g_{11} = x_{ξ}^{2} + y_{ξ}^{2}, g_{22} = x_{η}^{2} + y_{η}^{2}, g_{12} = x_{ξ}x_{η} + y_{ξ}y_{η}  (3.31)

for all (ξ, η) ∈ R.

3.6 Fixed boundary points

Consider a geodesic quadrangle P with sides A_{i}, i = 1, …, 4 and its four boundary points (x_{d}, y_{d}) ∈ A_{1}, (x_{r}, y_{r}) ∈ A_{2}, (x_{u}, y_{u}) ∈ A_{3} and (x_{l}, y_{l}) ∈ A_{4} such that A_{4} or A_{1} will not contain both points (x_{d}, y_{d}) and (x_{l}, y_{l}).

Then equations of the geodesics passing through (x_{d}, y_{d}) and (x_{u}, y_{u}), and through (x_{l}, y_{l}) and (x_{r}, y_{r}) respectively, are as follows:

a_{v}x + b_{v}y − [1 − δ(x^{2} + y^{2})] = 0,  (3.32)
a_{h}x + b_{h}y + [1 − δ(x^{2} + y^{2})] = 0.  (3.33)

Now we will find the point x = x_{vh}, y = y_{vh} of intersection of (3.32) and (3.33). Denoting

a = \frac{a_{h}b_{v} − a_{v}b_{h}}{2} , b = (a_{h} + a_{v})^{2} + (b_{h} + b_{v})^{2},

we find:

x_{vh} = \frac{−(b_{h} + b_{v})}{a + \sqrt{a^{2} + δb}},  y_{vh} = \frac{a_{h} + a_{v}}{a + \sqrt{a^{2} + δb}}.  (3.34)

In order to calculate the coefficients g_{ik} of the metric tensor of the mentioned quasi-isometric transformation (3.34) that maps R onto P we can use different explicit formulas. For this purpose we place one of the vertices of a cell of the grid in the origin by means of the motion group (2.2) and use formulas (3.6) or (3.27)—(3.28), (3.29)—(3.31).

4 Variational method for the generation of quasi-isometric grids

4.1 A certain class of Riemannian manifolds

Consider a geodesic quadrangle P with given angles α_{1}, …, α_{4} and sides of Euclidean lengths r_{1}, …, r_{4}. Let σ ∈ Σ_{4}, α = (α_{σ(1)}, …, α_{σ(4)}), and r = r_{σ(1)}. Consider one-parameter family of geodesic quadrangles P_{r} with angles α_{σ(1)}, …, α_{σ(4)} and parameter r ∈ (r_{σ(1)}^{min}, r_{σ(1)}^{max}). Every mapping of the form (3.4), (3.23)–(3.24) or (3.34) of the unit square R onto a quadrangle P_{r} has the following form:

x = x(ξ, η, r),  y = y(ξ, η, r),  (4.1)

and the mapping generates the class of conformally equivalent metrics. Let the following metric be a representative of the class:

ds^{2} = g_{11}(ξ, η, r)dξ^{2} + 2g_{12}(ξ, η, r)dξdη + g_{22}(ξ, η, r)dη^{2}.  (4.2)

So we can define the Riemannian manifold N (g_{ij}(ξ, η, r), R) with the coordinate domain R, the metric tensor g_{ij} and the parameter r. Then the Riemann mapping theorem implies the unique existence of the parameter r^{∗} such that the Riemannian manifold N (g_{ij}(ξ, η, r^{∗}), R) can be mapped onto given curvilinear quadrangle D conformally with respect to the metric (4.2). The mapping is quasi-isometric if all sides of D are smooth enough (belong to C^{2}) and in addition P_{r^{∗}} and D have the same angles [19, 14].

Thus the main problem consists of finding the parameter r^{∗} and a mapping X^{∗}(ξ, η), Y^{∗}(ξ, η) such that this mapping is conformal with respect to the metric (4.2) with metric tensor g_{ij}(ξ, η, r^{∗}).

4.2 Functional Φ

Using the the class of functions g_{ij} = g_{ij}(ξ, η, r) we can define the class of admitted functions A = A(ξ, η, r), B = B(ξ, η, r), C = C(ξ, η, r), where

A(ξ, η, r) = \frac{g_{22}(ξ,η,r)}{g}, B(ξ, η, r) = \frac{g_{12}(ξ,η,r)}{g},
C(ξ, η, r) = \frac{g_{11}(ξ,η,r)}{g}, g^{2} = g_{11}g_{22} − g_{12}^{2} .

Let us also introduce the class of admitted mappings X = X(ξ, η), Y = Y (ξ, η) of the computational region R onto D that has the following properties:

1. X(ξ, η), Y (ξ, η) define a quasi-isometric correspondence between ∂R and ∂D;

2. X(ξ, η), Y (ξ, η) can be continued inside R in such a way that the functional

Φ(X, Y, r) = ∫_{0}^{1}∫_{0}^{1} A(ξ, η, r)[X_{ξ}^{2} + Y_{ξ}^{2}] − 2B(ξ, η, r)[X_{ξ}X_{η} + Y_{ξ}Y_{η}]
+C(ξ, η, r)[X_{η}^{2} + Y_{η}^{2}] dξ dη ,  (4.3)

is bounded.

The minimum value of the functional is equal to the area S_{D} of the domain D, as it was shown in [10]. Functions X^{∗}(ξ, η), Y^{∗}(ξ, η) from the described class and the number r^{∗} that provide the minimum of the functional Φ, give us the desired mapping of the Riemannian manifold N (g_{ij}(ξ, η, r), R) onto D.

4.3 The variational principle

In order to find X^{∗}(ξ, η), Y^{∗}(ξ, η) and r^{∗} we will construct a minimizing sequence {X^{n}, Y^{n}, r^{n}} that has the following property:

Φ(X^{n+1}, Y^{n+1}, r^{n+1}) < Φ(X^{n}, Y^{n}, r^{n}), lim_{n→∞} M(P_{r^{n}}) = M(P_{r^{∗}}) = M(D).

To begin the minimization process we assume that the functions X^{n}(ξ, η), Y^{n}(ξ, η) and r^{n} are known to us.

On the first step, using the Montel variational principle from the conformal mapping theory, we obtain r^{n+1} such that

Φ(X^{n}, Y^{n}, r^{n+1}) < Φ(X^{n}, Y^{n}, r^{n}).

In particular,

r^{n+1} = r^{n} + \frac{Λ(X^{n}, Y^{n}, r^{n}) − V (X^{n}, Y^{n}, r^{n})}{νΛ(X^{n}, Y^{n}, r^{n}) + μV (X^{n}, Y^{n}, r^{n})},  (4.4)

where

Λ^{2}(X, Y, r) = ∫_{0}^{1}∫_{0}^{1} A(ξ, η, r)[X_{ξ}^{2} + Y_{ξ}^{2}] dξ dη,

V^{2}(X, Y, r) = ∫_{0}^{1}∫_{0}^{1} C(ξ, η, r)[X_{η}^{2} + Y_{η}^{2}] dξ dη.

Then we calculate the new coefficients of the functional

A = A(ξ, η, r^{n+1}), B = B(ξ, η, r^{n+1}), C = C(ξ, η, r^{n+1}).

In (4.4) we can use μ = 1/\sqrt{g_{11}(1, 0, r)} and ν = −\frac{dr_{σ(4)}}{dr} /\sqrt{g_{22}(0, 1, r)} [7]. On the second step for computation of new approximation X^{n+1}, Y^{n+1} such

that

Φ(X^{n+1}, Y^{n+1}, r^{n+1}) < Φ(X^{n}, Y^{n}, r^{n+1}),

we use obtained A, B, C as coefficients of the elliptic equations that represent the variational Euler – Lagrange equations for the functional (4.3) being minimized on X and Y:

−\frac{∂}{∂ξ} A\frac{∂X}{∂ξ} − \frac{∂}{∂η} C \frac{∂X}{∂η} + (\frac{∂}{∂ξ} B \frac{∂X}{∂η} + \frac{∂}{∂η} B \frac{∂X}{∂ξ} ) = 0,  (4.5)
−\frac{∂}{∂ξ} A\frac{∂Y}{∂ξ} − \frac{∂}{∂η} C \frac{∂Y}{∂η} + (\frac{∂}{∂ξ} B \frac{∂Y}{∂η} + \frac{∂}{∂η} B \frac{∂Y}{∂ξ} ) = 0.  (4.6)

The solution of the system (4.5)–(4.6) with appropriate boundary conditions can be used as a new approximation X^{n+1}, Y^{n+1}

Figure 6.  Figure 7.

Figure 8.  Figure 9.

Steps 1 and 2 are to be repeated till the desired accuracy of determining the solution of the variational problem is achieved.

Note that the formula (4.4) is invariant with respect to the metric chosen, i.e., the iteration process does not depend on the way the metric depend on the varying parameter. This feature is one of major developments in our paper comparing with works [10, 5].

5 Grid examples

In Figure 6, the boundary points are free on the straight boundaries, fixed on the outer half-circle and ”same as opposite” on the wing surface. In figure 7, all boundary points are fixed. In figure 8, the grid is nearly orthogonal. And figure 9 demonstrates that the grid cells remain parallelograms near corners as the grid is refined.

Acknowledgments This research was supported by the Russian Fundamental Research Foundation under the grant 96-15-96290.

References

[1] P. P. Belinskii, S. K. Godunov, Y. B. Ivanov and I. K. Yanenko, USSR Comp. Math. Math. Phys. 15 (1975), 133

[2] M. Berger: Géométrie, CEDIC/Nathan, Paris (1977)

[3] L. Bers, F. John and M. Schechter: Partial Differential Equations, Interscience Publishers, New York (1964)

[4] A. Cartan: Geometrie des espaces de Riemann, Paris. (1928)

[5] G. A. Chumakov: Riemannian Metric of the Harmonic Parameterization of Geodesic Quadrangles on Surfaces of Constant Curvature. Proc. of Inst. of Math. SB RAS, Vol. 22, pp.133–151. (1992) (in Russian)

[6] G. A. Chumakov: Conformal Parameterization of Curvilinear Quadrangles by Means of Geodesic Quadrangles on Surfaces of Constant Positive Curvature. Siberian Math. J., Vol. 34(1), pp. 172–180. (1993)

[7] G. A. Chumakov, S. G. Chumakov: A Method for the 2-D Quasi-Isometric Regular Grid Generation, J. Comput. Phys. 143, 1-28 (1998)

[8] R. Duraiswami, A. Prosperetti, Orthogonal Mapping in two Dimensions, J. Comput. Phys. 98 (1992), 254-268

[9] P. R. Eiseman, Geometric Methods in Computational Fluid Dynamics, ICASE 80-11, NASA Langley Research Center (1980)

[10] S. K. Godunov, V. M. Gordienko and G. A. Chumakov: Variational Principle for 2-D Regular Quasi-isometric Grid Generation. Comp. Fluid Dyn., Vol. 5, p.99–118. (1995)

[11] S. K. Godunov and G. P. Prokopov, USSR Comp. Math. Math. Phys. 7 (1967),89

[12] S. K. Godunov, E. I. Romenskii and G. A. Chumakov: Grid generation in complex domains by means of quasi-conformal mappings, Proc. of Institute of Math, Vol 18, pp. 75–84, Novosibirsk, Nauka (1990)

[13] S. K. Godunov, L. V. Zabrodin, M. Ya. Ivanov, et. al: Numerical Solution of the Multidimensional Problems of Gas Dynamics. (in Russian), Nauka, Moscow (1976)

[14] V.M. Gordienko: Boundary properties of conformal and quasi-conformal mappings of polygons. Siberian Math. J., v.33, N 6, pp. 966–972. (1992)

[15] J. A. Jenkins: Univalent functions and conformal mappings, Springer-Verlag (1958)

[16] I. S. Kang, I. G. Leal, Orthogonal Grid Generation in 2D Domain via the Boundary Integral Technique, J. Comput. Phys. 102 (1992), 78-87

[17] A. K. Khamayseh, C. W. Mastin, J. Comput. Phys. 123, 394-401, (1996)

[18] M. A. Lavrentyev: Variational method for BVP for elliptic systems, Moscow (1962)

[19] M.A. Lavrentyev and B.V. Shabat, Hydrodynamic problems and their mathematical models. (in Russian), Nauka, Moscow (1973)

[20] C. W. Mastin and J. F. Thompson, Math. Anal. Appl. 62, 52 (1978)

[21] H. J. Oh, I. S. Kang, A non-iterative Scheme for Orthogonal Grid Generation with Control Function and Specified Boundary Correspondence on three Sides, J. Comput. Phys. 112 (1994), 138-148

[22] B. L. van der Waerden: Algebra, Springer-Verlag (1971)

[23] J. F. Thompson, Z. U. A. Warsi, C. W. Mastin: Numerical Grid Generation, North-Holland, New York (1985)

[24] J. F. Thompson, Z. U. A. Warsi, C. W. Mastin: Boundary-Fitted Coordinate Systems for Numerical Solution of Partial Differential Equations – A Review, J. Comput. Phys. 47, 1 (1982)

[25] Z. U. A. Warsi, Tensors and Differential Geometry Applied to Numerical Coordinate Generation, MSUU-EIRS 8-1, Mississippi State Univ. (1981)

[26] A. M. Winslow: Numerical Solution of the Quasilinear Poisson Equation in a Nonuniform Triangle Mesh, J. Comput. Phys. 1, 149-172 (1966)

Gennadii A. Chumakov Sobolev Institute of Mathematics, pr. acad. Koptyuga, 4 Novosibirsk 630090, Russia Email address: chumakov@math.nsc.ru

Sergei G. Chumakov Mathematics Department, University of Wisconsin – Madison Madison, WI 53703, U.S.A. E-mail address: chumakov@math.wisc.edu