\documentclass[twoside]{article}
\usepackage{amssymb} % used for R in Real numbers
\usepackage{epsf}    % to input ps figures
\pagestyle{myheadings}
\setcounter{page}{55}
\markboth{\hfil Harmonic Parameterization of Geodesic Quadrangles \hfil}%
{\hfil Gennadii A. Chumakov \&  Sergei G. Chumakov \hfil}
\begin{document}
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
{\sc  Differential Equations and Computational Simulations III}\newline
J. Graef, R. Shivaji, B. Soni \& J. Zhu (Editors)\newline
Electronic Journal of Differential Equations, Conference~01, 1997, pp. 55-79. \newline
ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp  147.26.103.110 or 129.120.3.113 (login: ftp)}
 \vspace{\bigskipamount} \\
 Harmonic Parameterization of Geodesic Quadrangles 
 on Surfaces of Constant Curvature  and 2-D Quasi-Isometric Grids 
\thanks{ {\em 1991 Mathematics Subject Classifications:} 65N50, 30C30.
\hfil\break\indent
{\em Key words and phrases:} regular grid generation, quasi-isometric
mappings, geodesic grids, \hfil\break\indent
geodesic quadrangles, surfaces of constant
curvature.
\hfil\break\indent
\copyright 1998 Southwest Texas State University  and University of
North Texas. \hfil\break\indent
Published November 12, 1998.} }

\date{}

\author{Gennadii A. Chumakov \&  Sergei G. Chumakov}
\maketitle

\begin{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 $\cal
  R$ and a given physical domain $\cal 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.
\end{abstract}

\renewcommand{\theequation}{\thesection.\arabic{equation}}
\newtheorem{theorem}{Theorem}
\def \dd#1#2 {\frac{d#1}{d#2}}
\def \DD#1#2 {\frac{\partial {#1}}{\partial {#2}}}
\def \dD { \partial {\cal D} }
\def \dR { \partial {\cal R} }
\def \sign {\mathop {\rm sign}}



\section{Introduction}

\subsection{Quasi-isometric grids}

The generation of a 2-D quasi-isometric grids in a given physical
region ${\cal D}$ (a curvilinear quadrangle with interior angles
$\beta_i$, $ 0 < \beta_i < \pi$, $i=1,\dots,4$, and a conformal
modulus ${\cal M}$) may be considered as a problem of construction of
the mapping
\begin{equation}
X=X(\xi,\eta),\ Y=Y(\xi,\eta)
\label{eq:gen_xy}
\end{equation}
between points $(\xi,\eta)$ of the computational region 
$${\cal R}=\left\{ (\xi,\eta): 0\le\xi\le 1, 0\le\eta\le 1\right\}$$ 
and points
$(X,Y)\in{\cal D}$ such that (\ref{eq:gen_xy}) is the unique solution
of the following boundary value problem (BVP): given a quasi-isometric
mapping between $\dR$ and $\dD$ to extend the mapping inside ${\cal
R}$ as a quasi-isometric solution of the appropriate Beltrami system
\begin{eqnarray}
\sqrt{g_{11}(\xi,\eta,r) g_{22}(\xi,\eta,r) -
        g_{12}^2(\xi,\eta,r)}~~
X_\xi &=&-g_{12}(\xi,\eta,r)~
        Y_\xi+g_{11}(\xi,\eta,r)~Y_\eta, \nonumber \\
\label{eq:belt} \\        
\sqrt{g_{11}(\xi,\eta,r) g_{22}(\xi,\eta,r) -
        g_{12}^2(\xi,\eta,r)}~~
X_\eta &=&-g_{22}(\xi,\eta,r)~
        Y_\xi+g_{12}(\xi,\eta,r)~Y_\eta.
%\\
%g^2(\xi,\eta,r) & = & g_{11}(\xi,\eta,r) g_{22}(\xi,\eta,r) -
%        g_{12}^2(\xi,\eta,r) 
\nonumber
\end{eqnarray}
The coefficients of the system (\ref{eq:belt}) contain one unknown
parameter $r$, restricted to the interval $(r^{\rm min},r^{\rm 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 ${\cal D}$ are not fixed and can move
along $\dD$.

As it is shown in \cite{Bers-John-Schechter-64}, the linear elliptic
 system (\ref{eq:belt}) can be treated as a condition for conformality
 of the mapping (\ref{eq:gen_xy}) with respect to the following
 Riemmanian metric defined on ${\cal R}$:
\begin{equation}
ds^2 = g_{11}(\xi,\eta,r)d\xi^2 + 2 g_{12}(\xi,\eta,r)d\xi d\eta
+ g_{22}(\xi,\eta,r) d\eta^2.
\label{eq:riem_metric}
\end{equation}
Suppose the functions $g_{ik}(\xi,\eta,r)$ satisfy the inequality
\begin{equation}
g_{11}(\xi,\eta,r)+g_{22}(\xi,\eta,r) \leq Q~
\sqrt{g_{11}(\xi,\eta,r) g_{22}(\xi,\eta,r) - g^2_{12}(\xi,\eta,r)},
\label{eq:QQC}
\end{equation}
and are continuously differentiable in ${\cal R}$. Then the mapping
(\ref{eq:gen_xy}) has non-vanishing Jacobian in $\cal R$ and is called
{\sl $Q$-quasiconformal}.  That means that the singular values
$\nu_1=\nu_1(\xi,\eta)$ and $\nu_2=\nu_2(\xi,\eta)$ of the matrix
$$
J= \left[ \begin{array}{cc} X_{\xi} &
X_{\eta} \\ Y_{\xi} & Y_{\eta}  \end{array} \right],
$$
being enumerated in decreasing order, satisfy the condition
$\nu_1(\xi,\eta)/\nu_2(\xi,\eta)\leq Q$ for all $(\xi,\eta) \in {\cal
R}$. The least possible number $Q$ with such property is called {\sl
the coefficient of quasi-conformality} of the mapping
(\ref{eq:gen_xy}) in the domain ${\cal R}$. Clearly, a
$1$-quasiconformal mapping is conformal.

Note that the angle $\alpha$ between the grid lines is defined by the
formula $\cos\alpha=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 (\ref{eq:belt}) generates a metric
conformally equivalent to (\ref{eq:riem_metric}), we are able to
control the grid angle and the cell sides ratio by defining the
coefficients $g_{ik}$ of the metric (\ref{eq:riem_metric}) 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 (\ref{eq:belt}) is
$\mu$-quasi-isometric.

\begin{itemize}
\item A mapping (\ref{eq:gen_xy}) is called {\sl $\mu$-quasi-isometric} if
there exists a constant $\mu\in[1,\infty)$ such that
$$ 
\mu = \max \left\{
\sup_{(\xi,\eta)\in\cal R} \nu_1(\xi,\eta),
\sup_{(\xi,\eta)\in\cal R} \frac{1}{\nu_2(\xi,\eta)}
\right\}.  
$$
The constant $\mu$ is called {\sl the coefficient of
quasi-isometricity} of the mapping (\ref{eq:gen_xy}) in the domain
${\cal R}$.
\end{itemize}

Under a $\mu$-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/\mu$ and $\mu$.  If $\mu=1$ then (\ref{eq:gen_xy}) is an
isometric mapping.  It is clear that $\mu$-quasi-isometric mapping
will also be $\mu^2$-quasiconformal. Note that the relation
(\ref{eq:QQC}) and the system (\ref{eq:belt}) imply that the mapping
(\ref{eq:gen_xy}) is quasi-isometric in any sub-region of $\cal R$.
As a rule a quasiconformal mapping is not quasi-isometric, that is,
function $\nu_1(\xi,\eta)$ or $1/\nu_2(\xi,\eta)$ 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 \cite{GP-1967,GZI-1976}.  The
development of the Godunov-Prokopov method is precisely described in
\cite{GZI-1976} and \cite{TWM-1982}.  In a recent paper Khamayseh and
Mastin \cite{KM-1996} 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 \cite{W-1981}, \cite{Eiseman-1980} and
\cite{TWM-1985}.

The orthogonal mapping technique has been investigated in recent
works \cite{Kang-Leal-1992} and \cite{Dur-Prosp-1992},
\cite{Oh--Kang-1994} from different points of view.

In \cite{Belinskii-1975} Belinsky et.~al.~proposed to
consider the special class of quasi-conformal mappings by defining
$g_{11}=e^{2q(\xi)}$, $g_{22}=e^{2p(\eta)}$,
$g_{12}=\sqrt{g_{11}}\sqrt{g_{22}}\cos [\beta(\eta)-\alpha(\xi)]$. In
other words, the proposed metric contained four arbitrary functions of
either $\xi$ or $\eta$. This metric was obtained as a result of the
Chebyshev mapping ${\cal R}$ onto a plain curvilinear parallelogram.
In a similar way in the paper \cite{GRCh-1990} 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 $\xi$, $\eta$ 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
\cite{Chum-1992,Chum-1993}. In a later paper by Godunov
et.~al.~\cite{GGC-1995} 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 quasi-isometric solution.  The metric coefficients
$g_{ik}$ will be obtained via mapping the computational region ${\cal
R}$ onto a special class of geodesic quadrangles on surfaces of
constant curvature.  We start by investigating some basic properties
of one-parameter families of such geodesic quadrangles with given
angles.

After the basic concepts are introduced, we give several different
quasi-isometric parameterizations of geodesic quadrangles, that is,
quasi-isometric mappings
\begin{equation}
x=x(\xi,\eta,\alpha,{\bf r}),\ y=y(\xi,\eta,\alpha,{\bf r})
\label{eq:gen_xy(r)}
\end{equation} 
between ${\cal R}$ and the geodesic quadrangle ${\cal P}$. These
parameterizations have two vectors of parameters: interior angles
${\bf \alpha}=(\alpha_1,\dots,\alpha_4)$ and so-called Euclidean
lengths of sides ${\bf 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 $\cal
P$
$$
\delta^*=\alpha_1+\alpha_2+\alpha_3+\alpha_4-2\pi
$$ 
is positive or negative.  The mapping (\ref{eq:gen_xy}) generates the
metric tensor $g_{ij}$, coefficients of which are to be used as
coefficients of the Beltrami system (\ref{eq:belt}).

In \cite{Chum-chum-1998} the authors proposed several basic
parameterizations of a geodesic quadrangle ${\cal P}$, and the method
for finding the parameter ${\bf 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 ${\cal P}$ provides a
way to construct a mapping (\ref{eq:gen_xy}) which maps uniform grid
in the unit square on $(\xi,\eta)$-plane into a grid in ${\cal 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 ${\cal P}$ can be treated as a
variant of the Winslow grids \cite{Winslow-1966} with an advantage
that in our case the grid in ${\cal P}$ can be defined explicitly.
These kinds of two-dimensional harmonic mappings have non-vanishing
Jacobian, as it was shown in \cite{MT-1978} and \cite{GZI-1976}.  In
particular, we describe a parameterization (\ref{eq:gen_xy}) which
generates a metric with coefficients $g_{ik}(\xi,\eta,\alpha,{\bf r})$
that have the following properties:
\begin{equation}
\frac
{g_{11}(\xi,\eta,{\bf \alpha},{\bf r})}
{g_{22}(\xi,\eta,{\bf \alpha},{\bf r})}
=
\frac
{g_{22}(1-\eta,\xi,\sigma({\bf \alpha}),\sigma({\bf r}))}
{g_{11}(1-\eta,\xi,\sigma({\bf \alpha}),\sigma({\bf r}))},
\label{eq:metric1}
\end{equation}
\begin{eqnarray}
\lefteqn{ \frac
{g_{12}(\xi,\eta,{\bf \alpha},{\bf r})}
{\sqrt{ g_{11}(\xi,\eta,{\bf \alpha},{\bf r})~ g_{22}(\xi,\eta,{\bf
\alpha},{\bf r}) } }  }\label{eq:metric2}\\
&=&
\frac
{- g_{12}(1-\eta,\xi,\sigma({\bf\alpha}),\sigma({\bf r}))}
{\sqrt{ g_{11}(1-\eta,\xi,\sigma({\bf\alpha}),\sigma({\bf r}))~
g_{22}(1-\eta,\xi,\sigma({\bf\alpha}),\sigma({\bf r})) }}.
\nonumber
\end{eqnarray}
where $\sigma(1,2,3,4)=(4,1,2,3)$.

From equalities (\ref{eq:metric1}) and (\ref{eq:metric2}) it follows
that all elements of vectors $\alpha$ and ${\bf r}$ are equivalent
parameters and consequently any five of these parameters may be used
as characteristic invariants of ${\cal P}$. In particular, we shall
fix angles $\alpha_1,\dots,\alpha_4$ and study a set of geodesic
quadrangles depending on the parameter ${\cal M}$ for the interval
$0<{\cal M}<\infty$. In capacity of the parameter ${\cal M}$ we pick
the conformal module of ${\cal P}$.

We define the conformal module as follows.  For every geodesic
quadrangle ${\cal P}$ with vertices $z_i$, $i=1,\dots,4$, there
exists a unique ${\cal M} \in (0,\infty)$ such that there exists a
conformal mapping of ${\cal P}$ onto the rectangle 
$$\left\{(\xi,\eta):0\le\xi\le 1, 0\le\eta\le {\cal M}\right\}$$ 
under which the vertex
$z_1$ is mapped at the origin, and the other vertices of ${\cal P}$ go
over into vertices of the rectangle.

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

In later sections we describe a procedure of finding the mapping
$f(z)$, $z=x+iy$, by which the geodesic quadrangle ${\cal P}$ is
mapped conformally onto the physical domain ${\cal D}$; such conformal
mapping exists uniquely by virtue of the Riemann Mapping Theorem
\cite{Jen-1958}.  It is important to know how the module of the
derivative of the conformal mapping behaves on the boundary of ${\cal
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 ${\cal D}$.  Then by
virtue of the uniqueness of the conformal mapping of ${\cal P}$ onto
${\cal D}$ there exists a unique proimage $z$ of the point $b$ on
boundary of ${\cal P}$.  Consequently, we allow boundary points of
${\cal P}$ to move; we call such boundary conditions {\sl free} on
${\cal P}$, and {\sl fixed} on ${\cal D}$.  Similarly, if we allow
points of ${\cal D}$ to move along one of the boundaries, and fix
boundary points of ${\cal P}$ on the corresponding boundary, we call
such boundary conditions {\sl free} on ${\cal D}$ and {\sl fixed} on
${\cal P}$.

In order to construct a quasi-isometric grid we should take into
account well-known conditions of the existence of the boundary
derivative (see \cite{Lavrentyev-Shabat}, \cite{Lavrentyev} and
\cite{Gordienko}):
\begin{itemize}
\item If the boundary $\partial {\cal D}$ is twice continuously
  differentiable in some neighborhood of the point $w_0=f(\zeta_0)$
  and $\zeta_0$ is not a corner point of $\cal P$, then the derivative
  $f'$ can be continuously extended in certain part of the boundary
  $\partial {\cal P}$ , containing $\zeta_0$, in such a way that
  $f'(\zeta)\ne0$.
\item If $\zeta_0$ and $w_0=f(\zeta_0)$ are the corner points of
  ${\cal P}$ and ${\cal D}$ respectively and two twice continuously
  differentiable boundary arcs of ${\cal D}$ join in the point $w_0$
  at the same angle as do the corresponding arcs of ${\cal P}$ in the
  point $\zeta_0$, then $f'(\zeta)$ does not equal 0 and is bounded in
  a certain neighborhood of $\zeta_0$.
\end{itemize}
Thus, in order to use the proposed method for construction of a
quasi-isometric grid in ${\cal D}$, the following conditions should be
satisfied:
\begin{enumerate}
\item The boundary $\partial {\cal D}$ of the domain ${\cal D}$ must
  consist of four twice continuously differentiable curves;
\item The sought geodesic quadrangle ${\cal P}$ is to be chosen in such a way
  that the angles $\alpha_1,\dots,\alpha_4$ of ${\cal P}$ should
  coincide with the angles $\beta_1,\dots,\beta_4$ of the domain
  ${\cal D}$, and the conformal modules of ${\cal P}$ and ${\cal D}$
  must be the same.
\end{enumerate}

It was proved in \cite{GGC-1995,Chum-1992} and \cite{Chum-1993},
respectively, that under the restriction 
\begin{equation}
\delta^\ast < 2 \alpha_i, \qquad  i=1,\dots,4, 
\label{eq:cond_on_alpha}
\end{equation}
a geodesic quadrangle ${\cal P}$ satisfying the second condition does
exist uniquely in both cases of the negative and positive angle defect
$\delta^\ast$.

The sought mapping (\ref{eq:gen_xy}) 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
quasi-isometric grids with different boundary conditions are
presented.


% --- Pictures of the plane nose and the corresponding 
% --- geodesic quadrangle
\begin{center}
\begin{tabular}{c}
\epsfxsize=8cm \epsfbox{nose3.eps} \\
{\bf Figure 1a.}
\end{tabular}
\end{center} 

\begin{center}
\begin{tabular}{c}
\epsfxsize=7cm \epsfbox{nose3_gq.eps} \\
{\bf Figure 1b.}
\end{tabular}
\end{center}

In the Figure 1a we show a test domain ${\cal 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 ${\cal D}$ are fixed, then boundary
points of ${\cal 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(\zeta_0)$ split the boundary of ${\cal D}$ into two
$C^2$-curves that intersect at the angle $\beta=b \cdot \pi$,
$b\in[0,1]$. Then for the conformal mapping $f:{\cal
P}\rightarrow{\cal D}$ normed such that $\omega_0=f(\zeta_0) = 0$ in
neighborhood of $\zeta_0$ we have
$$
f(z) = a(z-\zeta_0)^b + o(|z-\zeta_0|^b), \quad
f'(z) = ab(z-\zeta_0)^{b-1} + o(|z-\zeta_0|^{b-1}),
$$
where $a\ne 0$.  In particular, it follows that at the point $\zeta_0$
the derivative $f'$ becomes infinite, and because of that the grid
points next to $f(\zeta_0)$ will not be close to $f(\zeta_0)$ as the
grid is refined.  Examples of quasi-conformal grids with singularities
of such kind can be found in the book \cite{TWM-1985}.

In conclusion we would like to say that from our point of view
Riemannian metrics that have properties (\ref{eq:metric1}) and
(\ref{eq:metric2}), 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 (\ref{eq:belt}) 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.


	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{One-parameter families of convex geodesic \\ quadrangles}
\setcounter{equation}{0}
	%--------------------------------
\subsection{Geodesic lines and bundles on surfaces of constant curvature}

Consider the surface of constant curvature $K=4\delta$ as the $(x,y)$-
parametric plane with the following metric:
\begin{equation}
ds^2=\frac{dx^2 + dy^2}{[1+\delta(x^2+y^2)]^2},
\label{eq:metric}
\end{equation} 
where $\delta$ is a real number \cite{Cartan}. If $\delta$ is
positive, the metric (\ref{eq:metric}) is defined for every $x$ and
$y$ including infinite, and we obtain a representation of spherical
geometry. If $\delta$ is negative, then in the disk $x^2+y^2+\delta<0$
we obtain the Poincare representation of Lobachevsky geometry, or
hyperbolic geometry.  In the case $\delta=0$, the metric
(\ref{eq:metric}) is the Euclidean metric on the plane $(x,y)$.

Let $q=ax+by+c[1-\delta(x^2+y^2)]$. Then geodesics in the metric
(\ref{eq:metric}) are curves defined by the equation $q=0$. Note that
for every $\delta\ne 0$ circles of the form $q=0$ are orthogonal to the
circle $1+\delta(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
\begin{equation}
w(z) = e^{i\omega}\frac{z - \zeta}{1+\delta\overline{\zeta} z},
\label{eq:move}
\end{equation}
where $\omega\in R$. The complex number $\zeta$ must satisfy the
condition $|\zeta|<|\delta|^{-1/2}$ in case $\delta<0$.

Let  $s_1$ and $s_2$ be two distinct geodesics defined by equations $q_1=0$
and $q_2=0$ respectively. The family ${\cal F}$ of geodesics orthogonal to
$s_1$ and $s_2$ is called a geodesic bundle.  The geodesic bundle ${\cal
F}^{\perp}$ in which $s_1$ and $s_2$ can be embedded is called orthogonal
to ${\cal F}$.

	%-------------------------
\subsection{Family of geodesic quadrangles with given angles}

Let ${\cal P}$ be a quadrangle whose sides are geodesics. We assume
further that the vertices $ z_i=(x_i,y_i),~i=1,\dots,4$ of the
quadrangle ${\cal P}$ are enumerated counter-clockwise. Let us denote
sides of ${\cal P}$ as $A_i=z_iz_{i+1}$, angles between $A_{i-1}$ and
$A_i$ as $\alpha_i$, $i=1,\dots,4$, (assume $z_5=z_1$ and $A_0=A_4$).
Let $\varphi_i = \alpha_i - \pi /2$, $i=1,\dots,4$.

It is possible to standardize a geodesic quadrangle ${\cal P}$ by
motions (\ref{eq:move}) 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
${\cal P}$ (say, $z_1$) to the origin by an appropriate transformation
(\ref{eq:move}).  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 (\ref{eq:metric}) under the motions (\ref{eq:move}) implies
that we can associate with every $A_i$ its unique Euclidean length
$r_i$.

In order to construct a geodesic quadrangle ${\cal P}$ it is
sufficient to have five parameters defined.  We shall fix angles
$\alpha_1,\dots,\alpha_4$ and study one-parameter families of
geodesic quadrangles ${\cal P}_{\cal M}$, ${\cal P}_m$ or ${\cal
P}_{r_1}$, where in the capacity of varying parameters we can choose
the conformal module ${\cal M}$, parameter $m=(r_1 r_3)/(r_2 r_4)$, or
a Euclidean side $r_1$.

From the results presented in \cite{Chum-1992,Chum-1993} it follows
that if we consider a one-parametric family ${\cal P}_{\cal M}$ then
invariants $m(\cal M)$ and $r_i(\cal M)$ depend monotonically on the
conformal module $\cal M$, which ranges from 0 to $\infty$.  In other
words, the following theorem holds.

\begin{theorem}
  Let ${\cal P}_{\cal M}$ and ${\cal P}_{\overline{\cal M}}$ be two
  geodesic quadrangles with the same angles and conformal modules
  $\cal M$ and $\overline{\cal M}$, respectively.  If $\cal M =
  \overline{\cal M}$ then ${\cal P}_{\cal M} = {\cal
    P}_{\overline{\cal M}}$; if $\cal M > {\overline{\cal M}}$ then
  $m({\cal M}) > m({\overline{\cal M}})$, $r_i({\cal M}) >
  r_i({\overline{\cal M}})$ and $r_{i+1}({\cal M}) <
  r_{i+1}({\overline{\cal M}})$ for $i=1,3$.
\label{theo:theorem1}
\end{theorem}

Below we provide some examples of convex geodesic quadrangles ${\cal
P}_{\cal M}$ with same angles and different conformal modules.

\begin{tabular}{cc}
\epsfxsize=5.2cm \epsfbox{o1.eps} & \epsfxsize=5.2cm\epsfbox{o2.eps} \\
{\bf a)} & {\bf b)} \\
\epsfxsize=5.2cm \epsfbox{o3.eps} & \epsfxsize=5.2cm\epsfbox{o4.eps} \\
{\bf c)} & {\bf d)} \\
\multicolumn{2}{c}
{{\bf Figure 2.}
Geodesic quadrangle ${\cal P}_{\cal M}$ with}\\
\multicolumn{2}{c}
{same angles  and various conformal module ${\cal M}$.}
\end{tabular}

Note that the parameter $r_1({\cal M})$ varies in the range from
$r_1^{min}=r_1(0)=0$ to $r_1^{\max}= r_1(\infty) < \infty$ as ${\cal
M}$ varies from $0$ to $\infty$ which leads to high parametric
sensitivity of conformal module ${\cal M}$ to $r_1$.  When ${\cal 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 ${\cal P}_{\cal M}$ as ${\cal M} \rightarrow 0$.
In Figure 2d we see that $r_1({\cal M}) \rightarrow r^{max}_1 $ and
$r_2({\cal M}) \rightarrow 0 $ as ${\cal M} \rightarrow \infty $.


\subsection{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 ${\cal P}_{\cal M}$ by given conformal module
$\cal M$, which takes into account the high parametric sensitivity of
$\cal M$ to $r_i$.

First of all for this purpose we shall find positive numbers $r_i^{\rm
min}$ and $r_i^{\rm max}$ such that $r_i(\cal M)$ belongs to the
interval $(r_i^{\rm min},r_i^{\rm max})$ for all ${\cal M}$, and
$r_i^{\rm min}=r_i(0)$, $r_i^{\rm max}=r_i(\infty)$ for
$i=1,\dots,4$.  As it was shown in (\cite{Chum-chum-1998}), we can
calculate $r_i^{\rm min}$ and $r_i^{\rm max}$ using the special
function $B(\psi_1,\psi_2,\psi_3,\psi_4,\delta)$ given by
\begin{equation}
B(\psi_1,\psi_2,\psi_3,\delta) = \sqrt{ \frac{\sin\psi
 \cos(\psi-\psi_3)} {\delta \cos(\psi-\psi_2) \sin(\psi-\psi_2-\psi_3)}},
\label{eq:funcb}
\end{equation}
where $\psi = (\psi_1+\psi_2+\psi_3+\pi/2)/2$.  It is convenient
for us to fix $\delta = \sin \gamma$, where $ \gamma =
(\varphi_1+\varphi_2+\varphi_3+\varphi_4)/2 $.

Let us denote by $l({\cal P})$ the number of sides of ${\cal P}$ the
sum of whose adjacent angles is not less then $\pi$. The value of
$l({\cal P})$ can be 0, 1, 2, 3 and 4. Let us consider each case
separately.  Let $\Sigma_4$ be the cyclic subgroup of the permutation
group ${\cal S}_4$, generated by the element
$$
\bar{\sigma} \in {\cal S}_4, \qquad
\bar{\sigma}=
\left(
\begin{array}{cccc}
1 & 2 & 3 & 4\\
2 & 3 & 4 & 1 \end{array}
\right),
$$

Let $l({\cal P})$ = 0, then for all $\sigma\in\Sigma_4$
\begin{eqnarray}
r_{\sigma(1)}^{\min} = B(\varphi_{\sigma(1)}, \varphi_{\sigma(2)},
                         -\pi/2,\delta), & &
r_{\sigma(1)}^{\max} = B(\varphi_{\sigma(1)}, \pi/2,
                         \varphi_{\sigma(4)}, \delta), \nonumber \\
r_{\sigma(2)}^{\min} = B(\varphi_{\sigma(2)}, \varphi_{\sigma(3)},
                         -\pi/2,\delta), & &
r_{\sigma(2)}^{\max} = B(\varphi_{\sigma(2)}, \pi/2,
                         \varphi_{\sigma(1)}, \delta), \nonumber \\
r_{\sigma(3)}^{\min} = B(\varphi_{\sigma(3)}, \varphi_{\sigma(4)},
                         -\pi/2,\delta), & &
r_{\sigma(3)}^{\max} = B(\varphi_{\sigma(3)}, \pi/2,
                         \varphi_{\sigma(2)}, \delta), \nonumber \\
r_{\sigma(4)}^{\min} = B(\varphi_{\sigma(4)}, \varphi_{\sigma(1)},
                         -\pi/2,\delta), & &
r_{\sigma(4)}^{\max} = B(\varphi_{\sigma(4)}, \pi/2,
                         \varphi_{\sigma(3)}, \delta), \nonumber
\end{eqnarray}

Let $l({\cal P}) = 1$, and $\sigma\in\Sigma_4$ be such that
$\varphi_{\sigma(1)}+\varphi_{\sigma(2)}\ge0$. Then the boundaries will be
as follows.
\begin{eqnarray*}
r_{\sigma(1)}^{\min} &=& 0, \\
r_{\sigma(1)}^{\max} &=& B(\varphi_{\sigma(1)}, \pi/2,
                            \varphi_{\sigma(4)}, \delta), \\
r_{\sigma(2)}^{\min} &=&
   B(\varphi_{\sigma(3)}, \varphi_{\sigma(2)},-\pi/2,\delta), \\
r_{\sigma(2)}^{\max} &=&
   B(\varphi_{\sigma(1)}+\varphi_{\sigma(2)}-\pi/2,
     \varphi_{\sigma(3)}, \varphi_{\sigma(4)},\delta), \\
r_{\sigma(3)}^{\min} &=&
   B(\varphi_{\sigma(3)},\varphi_{\sigma(4)},
     \varphi_{\sigma(1)}+\varphi_{\sigma(2)}-\pi/2,\delta), \\
r_{\sigma(3)}^{\max} &=&
   B(\varphi_{\sigma(3)},-\pi/2,\varphi_{\sigma(2)},\delta), \\
r_{\sigma(4)}^{\min} &=&
   B(\varphi_{\sigma(1)},\varphi_{\sigma(4)},-\pi/2,\delta), \\
 r_{\sigma(4)}^{\max} &=&
   B(\varphi_{\sigma(1)}+\varphi_{\sigma(2)}-\pi/2,
     \varphi_{\sigma(4)},\varphi_{\sigma(3)},\delta). 
\end{eqnarray*}

Consider now the case $l({\cal P})\ge 2$. There always exist such
$\sigma\in\Sigma_4$ that
$$
\varphi_{\sigma(1)} + \varphi_{\sigma(2)}  \ge
  \varphi_{\sigma(3)} + \varphi_{\sigma(4)}, \quad
\varphi_{\sigma(1)} + \varphi_{\sigma(4)}  \ge
  \varphi_{\sigma(2)} + \varphi_{\sigma(3)}
$$
holds. Then the boundaries for $r_j$ will be as follows:
\begin{eqnarray*}
r_{\sigma(1)}^{\min} &=& 0, \\
  r_{\sigma(1)}^{\max} &=&
  B(\varphi_{\sigma(1)}+\varphi_{\sigma(4)}-\pi/2,
      \varphi_{\sigma(2)},\varphi_{\sigma(3)},\delta), \\
r_{\sigma(2)}^{\min} &=&
  B(\varphi_{\sigma(2)},\varphi_{\sigma(3)},
      \varphi_{\sigma(1)}+\varphi_{\sigma(4)}-\pi/2,\delta), \\
 r_{\sigma(2)}^{\max} &=&
  B(\varphi_{\sigma(1)}+\varphi_{\sigma(2)}-\pi/2,
    \varphi_{\sigma(3)},\varphi_{\sigma(4)},\delta), \\
r_{\sigma(3)}^{\min} &=&
  B(\varphi_{\sigma(3)},\varphi_{\sigma(4)},
    \varphi_{\sigma(1)}+\varphi_{\sigma(2)}-\pi/2,\delta), \\
r_{\sigma(3)}^{\max} &=&
  B(\varphi_{\sigma(1)}+\varphi_{\sigma(4)}-\pi/2,
      \varphi_{\sigma(3)},\varphi_{\sigma(2)},\delta),   \\
r_{\sigma(4)}^{\min} &=& 0,\\
r_{\sigma(4)}^{\max} &=&
  B(\varphi_{\sigma(1)}+\varphi_{\sigma(2)}-\pi/2,
    \varphi_{\sigma(4)},\varphi_{\sigma(3)},\delta).
\end{eqnarray*}


\subsection{Relations between Euclidean lengths}

The relations $r_{\sigma(4)}(r_{\sigma(1)})$ and
$dr_{\sigma(4)}/d(r_{\sigma(1)})$ for all $\sigma \in \Sigma_4$ can be
obtained from the following equation which involves parameters
$r_{\sigma(1)}$, $r_{\sigma(4)}$, $\varphi_1$, $\varphi_2$,
$\varphi_3$, $\varphi_4$ and $\delta$:
\begin{eqnarray}
C_{\sigma(3)} S_0 &=& \delta C_{\sigma(2)} S_{\sigma(23)}
r_{\sigma(1)}^2 + \delta C_{\sigma(4)} S_{\sigma(34)}
r_{\sigma(4)}^2 \label{eq:baseq}\\
&& + 2 \delta D_\sigma r_{\sigma(1)} r_{\sigma(4)} - \delta^2
C_{\sigma(1)} S_{\sigma(24)} r_{\sigma(1)}^2 r_{\sigma(4)}^2,
\nonumber
\end{eqnarray}
where
$$
S_0    = \sin \gamma, \quad
D_\sigma = \cos\varphi_{\sigma(2)} \cos\varphi_{\sigma(4)}, \quad
C_i = \cos(\gamma - \varphi_i), \quad
$$
$$
S_{ij} = \sin(\gamma - \varphi_i - \varphi_j), \quad
\sigma(ij) =  \left( \sigma(i), \sigma(j) \right), \quad
i,j=1,\dots,4.
$$

It follows from (\ref{eq:baseq}) that for $\delta \ge 0$:
\begin{equation}
r_{\sigma(4)}(r_{\sigma(1)}) = \frac{1}{B + \sqrt{B^2+C}},
\label{eq:r4(r1)}
\end{equation}
where
$$
B = \frac{r_{\sigma(1)} D_\sigma}
         {C_{\sigma(3)} - r_{\sigma(1)}^2 C_{\sigma(2)} S_{\sigma(23)}},
\qquad
C = \frac{r_{\sigma(1)}^2 C_{\sigma(1)} S_{\sigma(24)} \delta -
         C_{\sigma(4)} S_{\sigma(34)}}
     {C_{\sigma(3)} - r_{\sigma(1)}^2 C_{\sigma(2)} S_{\sigma(23)}}.
$$

In the case $\delta < 0$ we have to take as $r_{\sigma(4)}$ the root
of (\ref{eq:baseq}) that belongs to the interval
$(r_{\sigma(4)}^{\min},r_{\sigma(4)}^{\max})$.

Thus given $r_{\sigma(1)}$ and four angles $\alpha_1,\dots,\alpha_4$
we are able by means of (\ref{eq:baseq}) to determine the function
$r_j=r_j(r_i)$, and its derivative with respect to $r_i$ for
$i,j=1,\dots,4$. In order to do that it is sufficient to find the
derivative of $r_{\sigma(4)}(r_{\sigma(1)})$ with respect to
$r_{\sigma(1)}$:
\begin{equation}
\dd{r_{\sigma(4)}(r_{\sigma(1)})}{r_{\sigma(1)}} = -
\frac{C_{\sigma(2)} S_{\sigma(23)} r_{\sigma(1)} +
      D_\sigma r_{\sigma(4)} - \delta C_{\sigma(1)}
      S_{\sigma(24)} r_{\sigma(1)} r_{\sigma(4)}^2 }
     {C_{\sigma(4)} S_{\sigma(34)} r_{\sigma(4)} +
      D_\sigma r_{\sigma(1)} - \delta C_{\sigma(1)}
      S_{\sigma(24)} r_{\sigma(1)}^2 r_{\sigma(4)} }.
\label{eq:dr4dr1}
\end{equation}
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{ Quasi--isometric parameterizations of \\ geodesic quadrangles}
\setcounter{equation}{0}
It should be noted that any five parameters from vectors
${\alpha}=(\alpha_1,\dots,\alpha_4)$ and ${\bf r}=(r_1,r_2,r_3,r_4)$
may be used as characteristic invariants of a geodesic quadrangle
$\cal P$. We shall now suppose that parameters $\varphi_1, \varphi_2,
\varphi_4, r_1$ and $r_4$ are given.


\subsection{The simplest quasi--isometric parameterization}

Consider $\cal P$ as an intersection of two geodesic bundles, which are
described by the following algebraic equations
\begin{eqnarray}
&(\cos\varphi_1 + \xi a_1) x + (\sin\varphi_1 - \xi b_1) y - \xi
c_1 [1-\delta(x^2+y^2)] =0, \qquad \xi\in[0,1],&\label{eq:vlines}\\
&\eta a_2 x - (1+\eta b_2)y + \eta c_2[1-\delta(x^2+y^2)] = 0,
\qquad \eta\in[0,1],&
\label{eq:hlines}
\end{eqnarray}
where
\begin{eqnarray*}
&a_1  =  \cos\varphi_2 (1 - \delta r_1^2) - \cos\varphi_1, \quad
b_1  =  \sin\varphi_2 (1 + \delta r_1^2) + \sin\varphi_1, &\\
&c_1  =  r_1 \cos \varphi_2, \quad a_2  =  \sin(\varphi_1+\varphi_4) -
          \delta r_4^2 \sin(\varphi_1-\varphi_4), &\\
&b_2  =  \cos(\varphi_1+\varphi_4) -
          \delta r_4^2 \cos(\varphi_1-\varphi_4) - 1, \quad
c_2  =  r_4 \cos \varphi_4.&
\end{eqnarray*}
Let us denote by
\begin{eqnarray}
P(\xi) & = & c_2 \cos \varphi_1 + \xi (a_1 c_2 + a_2 c_1), \nonumber\\
Q(\xi,\eta) &=& \xi c_1 - \eta c_2 \sin\varphi_1 +
                  \xi\eta (b_1 c_2 + b_2 c_1),  \nonumber \\
F(\xi,\eta) & = & \frac{1}{2}
\left[
\cos\varphi_1 + \xi a_1 + \eta (a_2 \sin\varphi_1 +
b_2 \cos\varphi_1) + \xi \eta (a_1 b_2 - a_2 b_1)
\right], \label{eq:pqf} \\
2W(\xi,\eta) & = & F(\xi,\eta) + \sqrt{F(\xi,\eta)^2+
               \delta(Q(\xi,\eta)^2+\eta^2P(\eta)^2)}, \nonumber
\end{eqnarray}
then by virtue of the relations
(\ref{eq:vlines})-(\ref{eq:hlines}) we can obtain the
following mapping of the unit square $\cal R$ onto
the geodesic quadrangle  $\cal P$:
\begin{eqnarray}
x (\xi,\eta) =  \frac{Q(\xi,\eta)}{2W(\xi,\eta)}, \qquad
y (\xi,\eta) = \frac{\eta\cdot P(\xi)}{2W(\xi,\eta)}.
\label{eq:xy}
\end{eqnarray}
Note that if the intersection of geodesic bundles
(\ref{eq:hlines}) and (\ref{eq:vlines}) is a geodesic convex
quadrangle then for all  $\xi\in[0,1]$ the inequality $P(\xi)>0$ holds
and the mapping (\ref{eq:xy}) is a quasi-isometric one.

From (\ref{eq:xy}) we have by differentiation
\begin{eqnarray}
x_\xi = {\cal A} [  (1+\eta b_2) W + \delta \eta^2 c_2 P] z, \quad
x_\eta = P[(b_1 \xi - \sin \varphi_1)W - \delta\xi\eta c_1 P] z,\nonumber \\
y_\xi = \eta \cdot {\cal A} (a_2 W -\delta c_2 Q) z, \quad
y_\eta = P [(\cos \varphi_1+\xi a_1) W + \delta\xi c_1 Q] z,
\label{eq:jacobi}
\end{eqnarray}
where
$P=P(\xi)$, $Q=Q(\xi,\eta)$, $F=F(\xi,\eta)$ and $W=W(\xi,\eta)$
were defined by (\ref{eq:pqf})  and
\begin{eqnarray*}
&{\cal A} = c_1 \cos \varphi_1 +
        \eta \left[ (a_1c_2+a_2c_1) \sin\varphi_1 +
        (b_1c_2+b_2c_1)\cos\varphi_1 \right], &\\
&4z = W^{-2} [F^2+\delta(Q^2+\eta^2P^2)]^{-1/2}.&
\end{eqnarray*}
We should mention that the mapping (\ref{eq:xy}) 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
\begin{eqnarray}
&g_{11}={\cal A}^2 \{ [
   (1+\eta b_2) W + \delta \eta^2 c_2 P ]^2 +
   \eta^2 (  a_2 W -\delta c_2 Q )^2 \},&\label{eq:GGG} \\
&g_{22}=P^2 \{ [
   (b_1 \xi - \sin\varphi_1) W - \delta\xi\eta c_1 P ]^2 +
  [  (\cos \varphi_1+\xi a_1) W +\delta\xi c_1 Q]^2
  \},& \nonumber \\
&g_{12}={\cal A} P \{ [ (1+\eta b_2) W + \delta \eta^2 c_2 P ] [
(b_1 \xi - \sin\varphi_1) W - \delta\xi\eta c_1 P ]& \nonumber\\
&+\eta^2 (a_2 W -\delta c_2 Q )
  [  (\cos \varphi_1+\xi a_1) W +\delta\xi c_1 Q]
  \}. &\nonumber
\end{eqnarray}

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 ${\cal M}$, but
with different enumerations of the sides and angles.

\begin{tabular}{cc}
\epsfxsize=5.2cm \epsfbox{u1.eps} & \epsfxsize=5.2cm \epsfbox{u2.eps} \\
{\bf a)} & {\bf b)} \\
\epsfxsize=5.2cm \epsfbox{u3.eps} & \epsfxsize=5.2cm \epsfbox{u4.eps} \\
{\bf c)} & {\bf d)} \\
\multicolumn{2}{c}
{{\bf Figure 3.} The simplest quasi-isometric parameterization of 
${\cal P}$}
\end{tabular}

\subsection{``Even" quasi-isometric parameterization}

Note that in the transformation (\ref{eq:xy}) we can take in the
capacity of $\xi$ and $\eta$ two arbitrary monotonically increasing
functions $\xi(\xi_1)$ and $\eta(\eta_1)$ satisfying the conditions
$\xi(0)=\eta(0)=0$ and $\xi(1)=\eta(1)=1$. In particular, we can
choose $\xi=\xi(\xi_1)$ and $\eta=\eta(\eta_1)$ in such a way that
under the mapping (\ref{eq:xy}) 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 ${\cal P}$ is also uniform in a sense of the Euclidean
distance.  In order to obtain such a mapping it is sufficient to
choose $\xi(\xi_1)$ and $\eta(\eta_1)$ as follows:
\begin{eqnarray*}
     &\xi(\xi_1)  =  \frac{\xi_1 r_1\cos\varphi_1}
     {c_1(1-\delta\xi_1^2 r_1^2)-a_1 \xi_1 r_1}\,, &\\
&\eta(\eta_1) = \frac{\eta_1 r_4 \cos \varphi_1}
     {c_2 (1-\delta \eta_1^2 r_4^2) -
      \eta_1 r_4 (a_1 \sin \varphi_1 + b_2 \cos \varphi_1)}\,.&
\end{eqnarray*}
Examples of curvilinear coordinate systems in the geodesic quadrangle on
a surface of the constant negative curvature generated by means of
``even" quasi-isometric parameterization are given below.

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\epsfxsize=5.2cm \epsfbox{n1.eps} & \epsfxsize=5.2cm \epsfbox{n2.eps} \\
{\bf a)} & {\bf b)} \\
\epsfxsize=5.2cm \epsfbox{n3.eps} & \epsfxsize=5.2cm \epsfbox{n4.eps} \\
{\bf c)} & {\bf d)} \\
\multicolumn{2}{c}
{{\bf Figure 4.}
The ``even" quasi-isometric parameterization.}
\end{tabular}
\end{center}\end{figure}

Figure 4 represents four geodesic quadrangles with the same angles and
the same "Euclidean" lengths of sides and, consequently, the same
conformal module ${\cal M}$.  In the Figure 4a is shown the geodesic
quadrangle with parameters $(\varphi_1,\varphi_2,\varphi_3,\varphi_4)
= (0.5, -1.0, -0.8, -1.0).$

Note that the grids in the geodesic quadrangle ${\cal P}$ constructed
by the ``even" and the simplest parameterizations are not invariant
under the cyclic permutation vertices of ${\cal 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.

	%----------------------------
\subsection{The harmonic parameterization of $\cal P$ on the Euclidean plane}

We now obtain parameterizations which are based upon geodesic bundles
and have the invariant properties
(\ref{eq:metric1})-(\ref{eq:metric2}).  We first consider the case
$\varphi_1+\varphi_2+\varphi_3+\varphi_4=0$, that is, $\delta=0$, so
that geodesics in the metric (\ref{eq:metric}) are straight lines
defined by the equation $ax+by+c=0$.  Consider ${\cal P}$ as an
intersection of two bundles, which are described by the following
equations:
\begin{eqnarray}
&x\cos[\varphi_1-\xi(\varphi_1+\varphi_2)]+
y\sin[\varphi_1-\xi(\varphi_1+\varphi_2)]-
r_1 \cos(\varphi_2)
\frac{\sin[\xi(\varphi_1+\varphi_2)]}{\sin(\varphi_1+\varphi_2)}= 0& \nonumber\\
& 0\leq\xi\leq 1,& \label{eq:ver} \\
&x\sin[\eta(\varphi_1+\varphi_4)] - y\cos[\eta(\varphi_1+\varphi_4)]
+r_4 \cos\varphi_4
\frac{\sin[\eta(\varphi_1+\varphi_4)]}
     {\sin(\varphi_1+\varphi_4)}=0 &\nonumber\\
& 0\leq\eta\leq 1 & \label{eq:gor}
\end{eqnarray}
In order to solve the system (\ref{eq:ver})--(\ref{eq:gor}) we denote by
\begin{eqnarray}
u(\xi,\eta)=\frac {r_1\cos\varphi_2 }
 {\cos[\varphi_1-\xi(\varphi_1+\varphi_2)-\eta(\varphi_1+\varphi_4)]}
 \frac{\sin[\xi(\varphi_1+\varphi_2)]}{\sin(\varphi_1+\varphi_2)},
\label{eq:u}
\end{eqnarray}
\begin{eqnarray}
v(\xi,\eta)=\frac
{r_4\cos\varphi_4}
 {\cos[\varphi_1-\xi(\varphi_1+\varphi_2)-\eta(\varphi_1+\varphi_4)]}
 \frac{\sin[\eta(\varphi_1+\varphi_4)]}{\sin(\varphi_1+\varphi_4)},
\label{eq:v}
\end{eqnarray}
then
\begin{eqnarray}
 x(\xi,\eta) = u(\xi,\eta) \cos[\eta(\varphi_1+\varphi_4)]-
     v(\xi,\eta) \sin[\varphi_1-\xi(\varphi_1+\varphi_2)],
\label{eq:xuv}\\
 y(\xi,\eta) = u(\xi,\eta) \sin[\eta(\varphi_1+\varphi_4)]+
     v(\xi,\eta) \cos[\varphi_1-\xi(\varphi_1+\varphi_2)],
\label{eq:yuv}
\end{eqnarray}
is the solution.

From (\ref{eq:xuv})--(\ref{eq:yuv})  we have by differentiation
\begin{eqnarray*}
x_{\xi}&=&\frac{V(\eta)\cos[\eta(\varphi_1+\varphi_4)] }
{\cos^2[\varphi_1-\xi(\varphi_1+\varphi_2)-\eta(\varphi_1+\varphi_4)]},\\
x_{\eta}&=&-\frac{U(\xi)\sin[\varphi_1-\xi(\varphi_1+\varphi_2)]}
{\cos^2[\varphi_1-\xi(\varphi_1+\varphi_2)-\eta(\varphi_1+\varphi_4)]},\\
y_{\xi}&=&-\frac{V(\eta)\sin[\eta(\varphi_1+\varphi_4)] }
{\cos^2[\varphi_1-\xi(\varphi_1+\varphi_2)-\eta(\varphi_1+\varphi_4)]},\\
y_{\eta}&=&\frac{U(\xi)\cos[\varphi_1-\xi(\varphi_1+\varphi_2)]}
{\cos^2[\varphi_1-\xi(\varphi_1+\varphi_2)-\eta(\varphi_1+\varphi_4)]},
\end{eqnarray*}
where
\begin{eqnarray}
U(\xi)&=&r_1(\varphi_1+\varphi_4)\cos\varphi_2
\frac{\sin[\xi(\varphi_1+\varphi_2)]}{\sin(\varphi_1+\varphi_2)}\nonumber\\
&&+r_4\cos\varphi_4\cos[\varphi_1-\xi(\varphi_1+\varphi_2)]
\frac{\varphi_1+\varphi_4}{\sin(\varphi_1+\varphi_4)},
\label{eq:1.17} \\
V(\eta)&=&r_1\cos\varphi_2\cos[\varphi_1-\eta(\varphi_1+\varphi_4)]
\frac{\varphi_1+\varphi_2}{sin(\varphi_1+\varphi_2)}\nonumber\\
&&+r_4(\varphi_1+\varphi_2)\cos\varphi_4
\frac{\sin[\eta(\varphi_1+\varphi_4)]}{\sin(\varphi_1+\varphi_4)}.
\label{eq:1.18}
\end{eqnarray}
It can happen that one obtains an indeterminate expression
of the form \newline
$\sin[\xi(\varphi_1+\varphi_2)]/\sin(\varphi_1+\varphi_2)$
as $\varphi_1+\varphi \rightarrow 0$ which is evidently
equal to $\xi$.

It is easily seen that for all $(\xi,\eta) \in {\cal R}$ the
inequalities $U(\xi)>0$ and $V(\eta)>0$  hold and the mapping
(\ref{eq:xuv})--(\ref{eq:yuv}) is a quasi-isometric transformation ${\cal R}$
onto ${\cal 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:
\begin{eqnarray}
&g_{11}=V^2(\eta),\quad
g_{22}=U^2(\xi), & \label{eq:1.19} \\
&g_{12}=U(\xi)V(\eta) \sin[-\varphi_1 +\xi(\varphi_1+\varphi_2)+
       \eta(\varphi_1+\varphi_4)].&\nonumber
\end{eqnarray}
Let us have a uniform grid in $\xi, \eta -$ plane, where $n, m -$ are the
total numbers of grid points to be used in $\xi, \eta $ directions. Then
the harmonic parameterization has a property that under the cyclic
permutation $(m,n) \rightarrow (n,m)$, 
$(\varphi_1,\varphi_2,\varphi_3,\varphi_4) \rightarrow
            (\varphi_4,\varphi_1,\varphi_2,\varphi_3)$,
$(r_1,r_2,r_3,r_4) \rightarrow (r_4,r_1,r_2,r_3)$,
the transformation (\ref{eq:xuv})--(\ref{eq:yuv}) results in the same grid
in ${\cal P}$.

The metric (\ref{eq:1.19}) by virtue of its simplicity can be used for
the generation of quasi-conformal grids in the physical domain ${\cal
D}$ with the angles satisfying the inequality
$\alpha_1+\alpha_2+\alpha_3+\alpha_4 \ne 2\pi$ as well. For this
purpose it is sufficient to substitute in
(\ref{eq:1.17})---(\ref{eq:1.19})
$$
\varphi_i=\alpha_i-\frac{\alpha_1+\alpha_2+\alpha_3+\alpha_4}{4},
\quad i=1,\dots,4.
$$

	%----------------------------
\subsection{The harmonic parameterization of ${\cal P}$ on a surface 
of constant curvature}

Let ${\cal P}$ be a geodesic quadrangle on a surface of the constant
curvature $K = 4 \delta $, where $\delta =
\sin[(\varphi_1+\varphi_2+\varphi_3+\varphi_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 $\cal P$, constants $\Delta_i$
and functions $SH(\xi w_i)$, $CH(\xi w_i)$:
\begin{eqnarray}
&\omega_i=
\left\{
\begin{array}{cc}
  {\rm arccot}(d_i/b_i),\quad{if}\quad l_i \geq 0,\\
  \coth^{-1}(d_i/b_i),\quad{if}\quad l_i <0,
\end{array}
\right.
\quad
\Delta_i=
\left\{ \begin{array}{cc}
    +1,\quad{if}\quad l_i\geq 0,\\
    -1,\quad{if}\quad l_i<0,
\end{array}\right.
&\nonumber\\
&SH(\xi\omega_i) =
\left\{ \begin{array}{cc}
    \sin (\xi\omega_i),\quad{if}\quad l_i\geq 0,\\
    \sinh(\xi\omega_i),\quad{if}\quad l_i<0,
    \end{array}
\right. &\nonumber \\
&CH(\xi\omega_i) =
\left\{ \begin{array}{cc}
    \cos (\xi\omega_i),\quad{if}\quad l_i \geq 0,\\
    \cosh(\xi\omega_i),\quad{if}\quad l_i <0
    \end{array}
\right.&
\label{eq:CH}
\end{eqnarray}
where $l_i=a^2_i+4\delta c_i^2$ and
\begin{eqnarray*}
&a_1=-\sin(\varphi_1+\varphi_2)+\delta r_1^2\sin(\varphi_1-\varphi_2), &\\
&a_2=\sin(\varphi_1+\varphi_4)-\delta r_4^2\sin(\varphi_1-\varphi_4),&\\
&b_1=\sqrt{|l_1|} \sign \delta, \quad
b_2=\sqrt{|l_2|} \sign \delta, &\\
&c_1=r_1\cos\varphi_2, \quad
c_2=r_4\cos\varphi_4, &\\
&d_1=\cos(\varphi_1+\varphi_2)-\delta r_1^2\cos(\varphi_1-\varphi_2),&\\ 
&d_2=\cos(\varphi_1+\varphi_4)-\delta r_4^2\cos(\varphi_1-\varphi_4).&
\end{eqnarray*}
Here we suppose that $\sign \delta=1$ if $\delta=0$.

Consider ${\cal P}$ as an intersection of two geodesic bundles
\begin{eqnarray}
&C(b_1,a_1,\xi\omega_1)x-S(b_1,a_1,\xi\omega_1)y+c_1 SH(\xi\omega_1)
[1-\delta(x^2+y^2)]=0,&\nonumber \\
& 0\leq\xi\leq 1,& \label{eq:vr}\\
&a_2 SH(\eta\omega_2)x - b_2 CH(\eta\omega_2)y +c_2 SH(\eta\omega_2)
[1-\delta(x^2+y^2)]=0,&\nonumber \\
&0\leq\eta\leq 1,&\label{eq:gr}
\end{eqnarray}
where
\begin{eqnarray}
&C(b_1,a_1,\xi\omega_1)=-b_1 \cos \varphi_1 CH(\xi\omega_1)
                       +a_1 \sin \varphi_1 SH(\xi\omega_1), &\label{eq:C}\\
&S(b_1,a_1,\xi\omega_1)= b_1 \sin \varphi_1 CH(\xi\omega_1)
                       +a_1 \cos \varphi_1 SH(\xi\omega_1). &\label{eq:S}
\end{eqnarray}
Introducing functions
\begin{eqnarray}
\bar{C}(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2)&=&
b_2CH(\eta w_2)C(b_1,a_1,\xi w_1)\nonumber \\
&&-a_2 SH(\eta w_2)S(b_1,a_1,\xi w_1), \label{eq:barC}\\
\bar{S}(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2)&=&
b_2CH(\eta w_2)S(b_1,a_1,\xi w_1)\nonumber \\
&&+a_2 SH(\eta w_2)C(b_1,a_1,\xi w_1). \label{eq:barS}
\end{eqnarray}
we can write the desired parameterization of ${\cal P}$ in the form
\begin{eqnarray}
\lefteqn{ x(\xi,\eta) } \label{eq:XX}\\
&=&\frac
{2c_2SH(\eta w_2)S(b_1,a_1,\xi w_1)-2c_1b_2SH(\xi w_1)CH(\eta w_2)}
{\bar{C}(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2)+
\Delta_0
\sqrt{\bar{C}^2(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2)+4\delta h(\xi,\eta)}}, 
\nonumber\\
\lefteqn{ y(\xi,\eta) }\label{eq:YY}\\
&=&\frac
{2 SH(\eta w_2)[c_2C(b_1,a_1,\xi w_1)-c_1a_2SH(\xi w_1)]}
{\bar{C}(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2)+
\Delta_0
\sqrt{\bar{C}^2(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2)+4\delta h(\xi,\eta)}},
\nonumber
\end{eqnarray}
where we have set
$$
\Delta_0={\rm sign}~\left\{SH(\eta w_2)[c_2C(b_1,a_1,\xi w_1)-c_1a_2SH(\xi
w_1)]\right\}, 
$$
\begin{eqnarray*}
 h(\xi,\eta)&=&
c_1^2SH^2(\xi w_1)[a_2^2SH^2(\eta w_2)+b_2^2CH^2(\eta w_2)]\\ 
&&+c_2^2SH^2(\eta w_2)[a_1^2SH^2(\xi w_1)+b_1^2CH^2(\xi w_1)]\\
&&-2c_1c_2SH(\xi w_1)SH(\eta w_2)\bar{S}(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2).
\end{eqnarray*}

\begin{center}
\begin{tabular}{cc}
\epsfxsize=5.2cm \epsfbox{hn1.eps} &  
\epsfxsize=5.2cm \epsfbox{hp1.eps} \\
\multicolumn{2}{c}{
{\bf Figure 5.}
An example of ``harmonic" parameterization.}
\end{tabular}
\end{center}

	%--------------------------
\subsection{Riemannian metric of the harmonic parameterization of ${\cal P}$}

The calculation of $x_\xi$ and $y_\xi$ is based on the relationships:
\begin{eqnarray}
\DD{}{\xi} C(b_1,a_1,\xi w_1)
&=&w_1S(a_1,b_1,\Delta_1\xi w_1), \label{eq:DDCS} \\
\DD{}{\xi} S(b_1,a_1,\xi w_1)
&=&-w_1C(a_1,b_1,\Delta_1\xi w_1), \nonumber
\end{eqnarray}
\begin{eqnarray}
\DD{}{\xi} \bar{C}(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2)&=&
w_1\bar{S}(b_2,a_2,a_1,b_1,\Delta_1\xi w_1,\eta w_2), \nonumber\\
\DD{}{\xi} \bar{S}(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2)&=&
-w_1\bar{C}(b_2,a_2,a_1,b_1,\Delta_1\xi w_1,\eta w_2), \nonumber\\
\DD{}{\eta} \bar{C}(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2)&=&
-w_2\bar{S}(a_2,b_2,b_1,a_1,\xi w_1,\Delta_2\eta w_2).
\label{eq:DDbarC}
\end{eqnarray}
We use the abbreviations
\begin{eqnarray*}
&\bar{C}=\bar{C}(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2), \quad
\bar{S}=\bar{S}(b_2,a_2,b_1,a_1,\xi w_1,\eta w_2), &\\
&h=h(\xi,\eta),\quad
\bar{C}_1=\bar{C}(b_2,a_2,a_1,b_1,\Delta_1\xi w_1,\eta w_2), &\\
&\bar{C}_2=\bar{C}(a_2,b_2,b_1,a_1,\xi w_1,\Delta_2\eta w_2), \quad
z_0=\Delta_0\sqrt{\bar{C}^2+4\delta h},& \\
&\bar{S}_1=\bar{S}(b_2,a_2,a_1,b_1,\Delta_1\xi w_1,\eta w_2), \quad
z_1=\bar{C}+z_0\,, &\\
&\bar{S}_2=\bar{S}(a_2,b_2,b_1,a_1,\xi w_1,\Delta_2\eta w_2),&
\end{eqnarray*}
and find
\begin{eqnarray*}
\lefteqn{ \DD{}{\xi} \left[C(b_1,a_1,\xi w_1)/z_1\right] }\\
 &=& w_1\left[S(a_1,b_1,\Delta_1\xi w_1)\bar{C}-
C(b_1,a_1,\xi w_1)\bar{S}\right]/(z_1 z_0) \\
&&+2\delta\left[2w_1hS(a_1,b_1,\Delta_1\xi w_1)-
h_\xi\bar{C}(b_1,a_1,\xi w_1)\right]/(z_1^2 z_0),
\end{eqnarray*}
\begin{eqnarray*}
\DD{}{\xi} \left[SH(\xi w_1)/z_1\right]&=&
w_1\left[CH(\xi w_1)\bar{C}-
SH(\xi w_1)\bar{S_1}\right]/(z_1 z_0) \\
&&+2\delta\left[2w_1 h CH(\xi w_1)-
h_\xi SH(\xi w_1) \right]/(z_1^2 z_0),
\end{eqnarray*}
\begin{eqnarray*}
\lefteqn{ \DD{}{\xi} \left[S(b_1,a_1,\xi w_1)/z_1\right] }\\
&=&-w_1\left[C(a_1,b_1,\Delta_1\xi w_1)\bar{C}+
S(b_1,a_1,\xi w_1)\bar{S}_1\right]/(z_1 z_0)\\
&& -2\delta\left[2w_1hC(a_1,b_1,\Delta_1\xi w_1)+
h_\xi S(b_1,a_1,\xi w_1)\right]/(z_1^2 z_0).
\end{eqnarray*}
It is easily seen that
\begin{eqnarray*}
&S(a_1,b_1,\Delta_1\xi w_1)\bar{C}-C(b_1,a_1,\xi w_1)\bar{S}=
-a_1 b_1 a_2 SH(\eta w_1),&\\
&CH(\xi w_1)\bar{C}-SH(\xi w_1)\bar{S_1}=
b_1 C(b_2,a_2,-\eta w_2),&\\
&C(a_1,b_1,\Delta_1\xi w_1)\bar{C}+S(b_1,a_1,\xi w_1)\bar{S}_1=
a_1 b_1 b_2 CH(\eta w_2),&
\end{eqnarray*}
and we therefore have
\begin{eqnarray}
\lefteqn{x_\xi}\label{eq:x_xi}\\
&=&-2 b_1 b_2 w_1 CH(\eta w_2)
\left[ c_1 C(b_2,a_2,-\eta w_2) + c_2 a_1 SH(\eta w_2) \right]
/(z_1 z_0)\nonumber\\
&&+4\delta\bigg\{\
2w_1 h [c_1 b_2 CH(\xi w_1) CH(\eta w_2)+
c_2 C(a_1,b_1,\Delta_1\xi w_1) SH(\eta w_2)] \nonumber \\
&&+h_\xi [-c_1 b_2 SH(\xi w_1) CH(\eta w_2)+
c_2 S(b_1,a_1,\xi w_1) SH(\eta w_2)]\bigg\}/(z_1^2 z_0),\nonumber
\end{eqnarray}
\begin{eqnarray}
y_\xi&=&-2 a_2 b_1 w_1 SH(\eta w_2)
\left[c_1 C(b_2,a_2,-\eta w_2) + c_2 a_1 SH(\eta w_2) \right]
/(z_1 z_0) \nonumber\\
&&-4\delta SH(\eta w_2)\left\{\
2 w_1 h [c_1 a_2 CH(\xi w_1)-c_2 S(a_1,b_1,\Delta_1\xi w_1)]-
 \right.\label{eq:y_xi} \\
&& \left.
-h_\xi [c_1 a_2 SH(\xi w_1)-c_2 C(b_1,a_1,\xi w_1)]
\right\} /(z_1^2 z_0),\nonumber
\end{eqnarray}
where
\begin{eqnarray*}
h_\xi&=& 2 w_1 \bigg\{
c_1^2 SH(\xi w_1) CH(\xi w_1)[b_2^2 CH^2(\eta w_2)+a_2^2 SH(\eta w_2)]\\
&&+c_2^2 SH^2(\eta w_2)SH(\xi w_1)CH(\xi w_1)[a_1^2-\Delta_1b_1^2]\\
&&-c_1 c_2 SH(\eta w_2)[CH(\xi w_1) \bar{S}_1-SH(\xi w_1)\bar{C}_1]\bigg\}.
\end{eqnarray*}
Next we calculate $x_\eta$ and $y_\eta$. Using (\ref{eq:DDbarC}) we find
\begin{eqnarray*}
\DD{}{\eta} \left[SH(\eta w_2)/z_1\right]&=&
w_2\left[CH(\eta w_2)\bar{C}+
SH(\eta w_2)\bar{S}_2\right]/(z_1 z_0)\\
&&+2\delta\left[2w_2 h CH(\eta w_2)-
h_\eta SH(\eta w_2) \right]/(z_1^2 z_0),\\
%
\DD{}{\eta} \left[CH(\eta w_2)/z_1\right]&=&
-w_2\left[\Delta_2 SH(\eta w_2)\bar{C}-
CH(\eta w_2)\bar{S}_2\right]/(z_1 z_0)\\ 
&&-2\delta\left[2\Delta_2 w_2 h SH(\eta w_2)+
h_\eta CH(\eta w_2) \right]/(z_1^2 z_0).
\end{eqnarray*}
It is also clear that
\begin{eqnarray*}
&CH(\eta w_2)\bar{C}+SH(\eta w_2)\bar{S}_2=
b_2 C(b_1,a_1,\xi w_1),& \\
&\Delta_2 SH(\eta w_2)\bar{C}-CH(\eta w_2)\bar{S}_2=
-a_2 S(b_1,a_1,\xi w_1),&
\end{eqnarray*}
and we therefore have
\begin{eqnarray}
x_\eta&=&2 b_2 w_2 S(b_1,a_1,\xi w_1)
[-c_1 a_2 SH(\xi w_1) + c_2 C(b_1,a_1,\xi w_1)]
/(z_1 z_0) \label{eq:x_eta} \\
&&+4\delta\{ 2 w_2 h [\Delta_2 c_1 b_2 SH(\xi w_1) SH(\eta w_2)+
c_2 S(b_1,a_1,\xi w_1) CH(\eta w_2)]\nonumber\\
&&+h_\eta [c_1 b_2 SH(\xi w_1) CH(\eta w_2)-
c_2 S(b_1,a_1,\xi w_1) SH(\eta w_2)]\}
/(z_1^2 z_0), \nonumber
\end{eqnarray}
\begin{eqnarray}
y_\eta&=&[c_2 C(b_1,a_1,\xi w_1)-c_1 a_2 SH(\xi w_1)]
\bigg\{ 2 b_2 w_2 C(b_1,a_1,\xi w_1)/(z_1 z_0)\nonumber\\
&&+4\delta [ 2 w_2 h CH(\eta w_2)-h_\eta SH(\eta w_2)]\bigg\}
/(z_1^2 z_0), \label{eq:y_eta}
\end{eqnarray}
where
\begin{eqnarray*}
h_\eta&=& 2 w_2 \bigg\{
c_1^2 SH^2(\xi w_1) SH(\eta w_2) CH(\eta w_2)
[a_2^2-\Delta_2 b_2^2] \\
&&+c_2^2 SH^2(\eta w_2) CH(\eta w_2)
[b_1^2 CH^2(\xi w_1)+a_1^2 SH^2(\xi w_1)]\\
&&-c_1 c_2 SH(\xi w_1)[CH(\eta w_2) \bar{S}+SH(\eta w_2)\bar{C}_2]\bigg\}.
\end{eqnarray*}
Thus, introducing functions (\ref{eq:CH})--(\ref{eq:S})
and (\ref{eq:barC}), (\ref{eq:barS}) we have found
$x_\xi, x_\eta, y_\xi, y_\eta$ and therefore coefficients of the metric
tensor
\begin{eqnarray}
g_{11}=x_\xi^2+y_\xi^2,\quad
g_{22}=x_\eta^2+y_\eta^2,\quad
g_{12}=x_\xi x_\eta+y_\xi y_\eta
\label{eq:ggg}
\end{eqnarray}
for all $(\xi, \eta) \in {\cal R}$.


\subsection{Fixed boundary points}

Consider a geodesic quadrangle ${\cal P}$ with sides $A_i$,
$i=1,\dots,4$ and its four boundary points
$(x_d,y_d)\in A_1$, $(x_r,y_r)\in A_2$, $(x_u,y_u)\in A_3$ and
$(x_l,y_l)\in 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:
\begin{eqnarray}
a_v x + b_v y - [1-\delta (x^2 + y^2)] & = & 0, \label{eq:vline1} \\
a_h x + b_h y + [1-\delta (x^2 + y^2)] & = & 0. \label{eq:hline1}
\end{eqnarray}

Now we will find the point $x=x_{vh},~y=y_{vh}$ of intersection of
(\ref{eq:vline1}) and (\ref{eq:hline1}). Denoting

$$ a = \frac{a_hb_v-a_vb_h}{2}, \qquad b =
(a_h+a_v)^2+(b_h+b_v)^2,
$$
we find:
\begin{equation}
x_{vh} = \frac{-(b_h+b_v)}{a +
\sqrt{a^2+\delta b}},
\qquad
y_{vh} = \frac{a_h+a_v}{a + \sqrt{a^2+\delta b}}.
\label{eq:xy_even}
\end{equation}
In order to calculate the coefficients $g_{ik}$ of the metric tensor
of the mentioned quasi-isometric transformation (\ref{eq:xy_even})
that maps $\cal R$ onto $\cal 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 (\ref{eq:move})
and use formulas (\ref{eq:GGG}) or (\ref{eq:x_xi})---(\ref{eq:y_xi}),
(\ref{eq:x_eta})---(\ref{eq:ggg}).
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Variational method for the generation of \\ quasi-isometric grids}
\setcounter{equation}{0}

\subsection{A certain class of Riemannian manifolds}

Consider a geodesic quadrangle ${\cal P}$ with given angles
$\alpha_1,\dots,\alpha_4$ and sides of Euclidean lengths
$r_1,\dots,r_4$. Let $\sigma\in\Sigma_4$,
$\alpha=(\alpha_{\sigma(1)},\dots,\alpha_{\sigma(4)})$, and
$r=r_{\sigma(1)}$.  Consider one-parameter family of geodesic quadrangles
${\cal P}_{r}$ with angles
$\alpha_{\sigma(1)},\dots,\alpha_{\sigma(4)}$ and parameter
$r\in(r_{\sigma(1)}^{\min},r_{\sigma(1)}^{\max})$.
Every mapping of the form (\ref{eq:xy}), (\ref{eq:XX})--(\ref{eq:YY}) or
(\ref{eq:xy_even}) of the unit square ${\cal R}$ onto a quadrangle
${\cal P}_{r}$ has the following form:
\begin{equation}
x=x(\xi,\eta,r), \qquad y=y(\xi,\eta,r),
\label{eq:ri-vid}
\end{equation}
and the mapping generates the class of conformally equivalent metrics. Let
the following metric be a representative of the class:
\begin{equation}
ds^2=g_{11}(\xi,\eta,r)
d\xi^2 + 2 g_{12}(\xi,\eta,r) d\xi d\eta +
g_{22}(\xi,\eta,r) d\eta^2.
\label{eq:metric_genr}
\end{equation}

So we can define the Riemannian manifold ${\cal N}(g_{ij}(\xi,\eta,r),
{\cal R})$ with the coordinate domain ${\cal R}$, the metric tensor
$g_{ij}$ and the parameter $r$.  Then the Riemann mapping theorem
implies the unique existence of the parameter $r^{\ast}$ such that the
Riemannian manifold ${\cal N}(g_{ij}(\xi,\eta,r^{\ast}),{\cal R})$ can
be mapped onto given curvilinear quadrangle ${\cal D}$ conformally
with respect to the metric (\ref{eq:metric_genr}). The mapping is
quasi-isometric if all sides of ${\cal D}$ are smooth enough (belong
to $C^2$) and in addition ${\cal P}_{r^{\ast}}$ and ${\cal D}$ have
the same angles \cite{Lavrentyev-Shabat,Gordienko}.

Thus the main problem consists of finding the parameter $r^{\ast}$ and
a mapping $X^{\ast}(\xi,\eta)$, $Y^{\ast}(\xi,\eta)$ such that this
mapping is conformal with respect to the metric (\ref{eq:metric_genr}) with
metric tensor $g_{ij}(\xi,\eta,r^{\ast})$.


\subsection{Functional $\Phi$}

Using the the class of functions $g_{ij}=g_{ij}(\xi,\eta,r)$ we can define
the class of admitted functions $A=A(\xi,\eta,r)$, $B=B(\xi,\eta,r)$,
$C=C(\xi,\eta,r)$, where
\begin{eqnarray*}
&A(\xi,\eta,r) = \frac{g_{22}(\xi,\eta,r)}{g}, \quad
B(\xi,\eta,r) = \frac{g_{12}(\xi,\eta,r)}{g}, & \\
&C(\xi,\eta,r) = \frac{g_{11}(\xi,\eta,r)}{g}, \quad
g^2 = g_{11}g_{22}-g_{12}^2\,.&
\label{eq:method_abc}
\end{eqnarray*}

Let us also introduce the class of admitted mappings $X=X(\xi,\eta),\
Y=Y(\xi,\eta)$ of the computational region ${\cal R}$ onto ${\cal D}$
that has the following properties:
\begin{enumerate}
\item $X(\xi,\eta),\ Y(\xi,\eta)$ define a quasi-isometric correspondence
between $\dR$ and $\dD$;
\item $X(\xi,\eta),\ Y(\xi,\eta)$ can be continued inside ${\cal R}$
in such a way that the functional
\begin{eqnarray}
\Phi(X,Y,r)&=&
\int\limits_0^1 \!  \int\limits_0^1
A(\xi,\eta,r) [X_\xi^2+Y_\xi^2] -
2B (\xi,\eta,r) [X_\xi X_\eta + Y_\xi Y_\eta] \nonumber\\
&&\qquad +C(\xi,\eta,r) [X_\eta^2+Y_\eta^2]\,d\xi\,d\eta\,,
\label{eq:functional}
\end{eqnarray}
is bounded.
\end{enumerate}

The minimum value of the functional is equal to the area $S_{\cal D}$
of the domain ${\cal D}$, as it was shown in \cite{GGC-1995}.
Functions $X^{\ast}(\xi,\eta)$, $Y^{\ast}(\xi,\eta)$ from the
described class and the number $r^{\ast}$ that provide the minimum of
the functional $\Phi$, give us the desired mapping of the Riemannian
manifold ${\cal N}(g_{ij}(\xi,\eta,r),{\cal R})$ onto ${\cal D}$.

\subsection{The variational principle}

In order to find $X^{\ast}(\xi,\eta)$, $Y^{\ast}(\xi,\eta)$ and
$r^{\ast}$ we will construct a minimizing sequence $\left\{ X^n, Y^n,
r^n \right\}$ that has the following property:
$$
\Phi(X^{n+1},Y^{n+1},r^{n+1}) < \Phi(X^n,Y^n,r^n), \quad
\lim_{n\to\infty} {\cal M} ({\cal P}_{r^n}) =
{\cal M} ({\cal P}_{r^{\ast}}) = {\cal M} ({\cal D}).
$$
To begin the minimization process we assume that the
functions $X^n(\xi,\eta)$, $Y^n(\xi,\eta)$ 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
$$
\Phi(X^{n},Y^{n},r^{n+1}) < \Phi(X^n,Y^n,r^n).
$$
In particular,
\begin{equation}
r^{n+1} = r^n +
\frac{\Lambda(X^n,Y^n,r^n) - V(X^n,Y^n,r^n)}
     {\nu \Lambda(X^n,Y^n,r^n) + \mu V(X^n,Y^n,r^n)},
\label{eq:r_new}
\end{equation}
where
\begin{eqnarray*}
\Lambda^2(X,Y,r) &=& \int\limits_0^1 \!  \int\limits_0^1
   A(\xi,\eta,r) [X_\xi^2+Y_\xi^2] \;d\xi\, d\eta, \\
V^2(X,Y,r) &=& \int\limits_0^1 \!  \int\limits_0^1
   C(\xi,\eta,r) [X_\eta^2+Y_\eta^2] \;d\xi\, d\eta.
\end{eqnarray*}
Then we calculate the new coefficients of the functional
$$
A = A(\xi,\eta,r^{n+1}), \quad
B = B(\xi,\eta,r^{n+1}), \quad
C = C(\xi,\eta,r^{n+1}).
$$

In (\ref{eq:r_new}) we can use $\mu = 1/\sqrt{g_{11}(1,0,r)}$ and $\nu
= -\dd{r_{\sigma(4)}}{r} /\sqrt{g_{22}(0,1,r)}$ \cite{Chum-chum-1998}.

On the second step for computation of new approximation
$X^{n+1}$, $Y^{n+1}$ such that
\begin{eqnarray}
\Phi(X^{n+1},Y^{n+1},r^{n+1}) < \Phi(X^n,Y^n,r^{n+1}), \nonumber
\end{eqnarray}
we use obtained $A$, $B$, $C$ as coefficients of the elliptic
equations that represent the variational Euler -- Lagrange equations
for the functional (\ref{eq:functional}) being minimized on $X$ and
$Y$:
\begin{eqnarray}
-\DD{}{\xi} A \DD{X}{\xi} - \DD{}{\eta} C \DD{X}{\eta} +
        \left(
        \DD{}{\xi} B \DD{X}{\eta} + \DD{}{\eta} B \DD{X}{\xi}
        \right) & = & 0, \label{eq:Euler_x} \\
-\DD{}{\xi} A \DD{Y}{\xi} - \DD{}{\eta} C \DD{Y}{\eta} +
        \left(
        \DD{}{\xi} B \DD{Y}{\eta} + \DD{}{\eta} B \DD{Y}{\xi}
        \right) & = & 0. \label{eq:Euler_y}
\end{eqnarray}

The solution of the system (\ref{eq:Euler_x})--(\ref{eq:Euler_y})
with appropriate boundary conditions can be used as a new
approximation $X^{n+1},Y^{n+1}$

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 (\ref{eq:r_new}) 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
\cite{GGC-1995,Chum-1992}.


\begin{figure}\begin{center}
 \begin{tabular}{cc}
  \epsfxsize=5cm \epsfclipon \epsfbox{wing1-10.eps} 
& \epsfxsize=6cm \epsfclipon \epsfbox{star30.eps} \\
{\bf Figure 6.} & {\bf Figure 7.} \end{tabular}\end{center}  \end{figure}



\begin{figure}\begin{center}
\begin{tabular}{cc}  
  \epsfxsize=5cm \epsfclipon \epsfbox{w1-10-2.eps} 
& \epsfxsize=5cm \epsfclipon \epsfbox{w1-10-1.eps}\\
{\bf Figure 8. }& {\bf Figure 9.}  \end{tabular} \end{center}\end{figure}



\section{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.


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


	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{thebibliography}{99}

\bibitem{Belinskii-1975}
P.~P.~Belinskii, S.~K.~Godunov, Y.~B.~Ivanov and I.~K.~Yanenko,
{\it USSR Comp.  Math.  Math.  Phys.}  15 (1975), 133

\bibitem{Berger}
M.~Berger: {\it G\'eom\'etrie}, CEDIC/Nathan, Paris (1977)

\bibitem{Bers-John-Schechter-64}
L.~Bers, F.~John and M.~Schechter: {\it Partial Differential Equations},
Inter-science Publishers, New York (1964)

\bibitem{Cartan}
A.~Cartan: {\it Geometrie des espaces de Riemann}, Paris. (1928)

\bibitem{Chum-1992}
G.~A.~Chumakov: Riemannian Metric of the Harmonic Parameterization of
Geodesic Quadrangles on Surfaces of Constant Curvature. {\it Proc. of
Inst. of Math. SB RAS,} Vol. 22, pp.133--151. (1992) (in Russian)

\bibitem{Chum-1993}
G.~A.~Chumakov: Conformal Parameterization of Curvilinear Quadrangles by
Means of Geodesic Quadrangles on Surfaces of Constant Positive Curvature.
{\it Siberian Math. J., Vol. 34(1)}, pp.~172--180. (1993)

\bibitem{Chum-chum-1998} 
G.~A.~Chumakov, S.~G.~Chumakov: A Method for the 2-D Quasi-Isometric
Regular Grid Generation, {\it J. Comput. Phys.} {\bf 143,} 1-28 (1998)

\bibitem{Dur-Prosp-1992}
R.~Duraiswami, A.~Prosperetti, {\it Orthogonal Mapping in two Dimensions},
J. Comput. Phys. 98 (1992), 254-268

\bibitem{Eiseman-1980}
P.~R.~Eiseman, {\it Geometric Methods in Computational Fluid Dynamics},
ICASE 80-11, NASA Langley Research Center (1980)

\bibitem{GGC-1995}
S.~K.~Godunov, V.~M.~Gordienko and G.~A.~Chumakov: Variational Principle
for 2-D Regular Quasi-isometric Grid Generation. {\it Comp.~Fluid Dyn.},
Vol.~5, p.99--118. (1995)

\bibitem{GP-1967} S.~K.~Godunov and G.~P.~Prokopov, {\it USSR
Comp. Math. Math. Phys.} 7 (1967),89

\bibitem{GRCh-1990}
S.~K.~Godunov, E.~I.~Romenskii and G.~A.~Chumakov: Grid generation in
complex domains by means of quasi-conformal mappings, {\it Proc. of
Institute of Math, Vol 18, pp.~75--84}, Novosibirsk, Nauka (1990)

\bibitem{GZI-1976}
S.~K.~Godunov, L.~V.~Zabrodin, M.~Ya.~Ivanov, et.~al:
{\it Numerical Solution of the Multidimensional Problems of Gas Dynamics.}
(in Russian), Nauka, Moscow (1976)

\bibitem{Gordienko}
V.M. Gordienko: Boundary properties of conformal and quasi-conformal
mappings of polygons. {\it Siberian Math. J.}, v.33, N 6, pp. 966--972.
(1992)

\bibitem{Jen-1958}
J.~A.~Jenkins: {\it Univalent functions and conformal mappings},
Springer-Verlag (1958)

\bibitem{Kang-Leal-1992}
I.~S.~Kang, I.~G.~Leal,
Orthogonal Grid Generation in 2D Domain via the Boundary Integral 
Technique, {\it J. Comput. Phys.} 102 (1992), 78-87

\bibitem{KM-1996} A.~K.~Khamayseh, C.~W.~Mastin, {\it
J. Comput. Phys.} 123, 394-401, (1996)

\bibitem{Lavrentyev}
M.~A.~Lavrentyev: {\it Variational method for BVP for elliptic systems},
Moscow (1962)

\bibitem{Lavrentyev-Shabat}
M.A. Lavrentyev and B.V. Shabat,
{\it Hydrodynamic problems and their mathematical models.} (in Russian),
Nauka, Moscow (1973)

\bibitem{MT-1978}
C.~W.~Mastin and J.~F.~Thompson, {\it Math. Anal. Appl.} 62, 52 (1978)

\bibitem{Oh--Kang-1994}
H.~J.~ Oh, I.~S.~Kang,
A non-iterative Scheme for Orthogonal Grid Generation with
Control Function and Specified Boundary Correspondence on three Sides,
{\it J. Comput. Phys.} 112 (1994), 138-148

\bibitem{Ward-1971}
B.~L.~van der Waerden: {\it Algebra}, Springer-Verlag (1971)

\bibitem{TWM-1985}
J.~F.~Thompson, Z.~U.~A.~Warsi, C.~W.~Mastin: {\it Numerical Grid
Generation},
North-Holland, New York (1985)

\bibitem{TWM-1982} J.~F.~Thompson, Z.~U.~A.~Warsi, C.~W.~Mastin:
Boundary-Fitted Coordinate Systems for Numerical Solution of Partial
Differential Equations -- A Review, {\it J. Comput. Phys.} 47, 1
(1982)

\bibitem{W-1981}
Z.~U.~A.~Warsi,  {\it Tensors and Differential Geometry Applied
to Numerical Coordinate Generation},
MSUU-EIRS 8-1, Mississippi State Univ. (1981)

\bibitem{Winslow-1966} A.~M.~Winslow: Numerical Solution of the
Quasilinear Poisson Equation in a Nonuniform Triangle Mesh, {\it
J. Comput. Phys.} 1, 149-172 (1966)

\end{thebibliography}

\bigskip

\noindent{\sc Gennadii A. Chumakov}\\
Sobolev Institute of Mathematics,  pr.~acad.~Koptyuga, 4 \\
Novosibirsk 630090, Russia \\
Email address: chumakov@math.nsc.ru \medskip
 
\noindent{\sc Sergei G. Chumakov}\\
Mathematics Department, University of Wisconsin -- Madison\\
Madison, WI 53703, U.S.A. \\
E-mail address: chumakov@math.wisc.edu
\end{document}