Adaptive Discrete Harmonic Grid Generation - Semantic Scholar
Adaptive Discrete Harmonic Grid Generation Pablo Barrera-S´anchez a, Longina Castellanos Noda b, Francisco J. Dom´ınguez-Mota c,∗, Guilmer F. Gonz´alez Flores a, Angel P´erez Dom´ınguez b a U.N.A.M.,
Facultad de Ciencias, Mexico City, Mexico. b ICIMAF,
c U.M.S.N.H.,
Havana City, Cuba.
Facultad de Ciencias F´ısico-Matem´ aticas, Morelia, Mexico.
Abstract We present a new adaptive-harmonic structured grid generation method. It is based on a functional that shares a common set of minimizers with Ivanenko’s harmonic functional [9]. An unconstrained optimization process related to a continuation parameter is used to guarantee the convexity of the grid cells. Several numerical examples of grids generated on quite irregular 2D regions, show the effective performance of the proposed method and its robustness when preset monitor functions are given. Key words: Variational grid generation, adaptive grids, quasi-harmonic grids, structured grids.
1
Introduction
Variational grid generation has been used quite successfully in the last years. It started with the work of Brackbill and Saltzman in 1982 [5], and the contributions of Steinberg and Roache in 1986 [19], who developed the basis for the continuous problem. Afterwards, the discrete version of the problem was first developed by Castillo [6].
∗ Corresponding author. Email addresses: [email protected] (Pablo Barrera-S´ anchez), [email protected] (Longina Castellanos Noda), [email protected] (Francisco J. Dom´ınguez-Mota), [email protected] (Guilmer F. Gonz´ alez Flores), [email protected] (Angel P´erez Dom´ınguez).
Preprint submitted to Mathematics and Computers in Simulations 24 April 2007
In the classical variational setting, the strategy was to obtain the optimal grid by solving -numerically in many cases- the Euler-Lagrange equations associated to a functional. There is a wealth of literature available in the subject, mainly devoted to the study of the theoretical properties of the functionals in order to guarantee that their minimizers are convex grids. It has been successfully applied to many regions with a number of symmetries or with an elongated geometry. Comprehensive treatments of the problem can be found in Castillo [6], de Almeida [10], and Knupp and Robidoux [15]. There are also some other approaches to the grid generation problem, for instance, by subdividing the region into simpler subregions, which transforms the problem of gridding an irregular complex domain into the efficient computational choice of a useful subdivision [1,18]. However, from our point of view, the structured gridding of highly irregular planar regions has not been satisfactorily explored. The need of solving modeling problems in domains as realistic as possible was indeed what led to Ivanenko and Charakhch’yan [13] in 1988 to solve the grid generation problem, using the discretization of the variational principle rather than the Euler Lagrange equations. This allows to control directly the interior nodes of the grid and obtain a solution through a large scale optimization method, known as the Direct Optimization Method [16]. Besides, in order to be able to work with the discretized functionals it was necessary to transform the requirement of the positiveness of the Jacobian of the mapping between the regions, that guarantees the convexity of the grid. They proposed a very simple check [8,9,13,14]: the convexity of each grid cell, which represents a major advance from the practical point of view. Ivanenko and Charakhch’yan’s method allowed to generate smooth grids with a great number of nodes, on very irregular regions, in a very efficient and practical way: by solving a large scale unconstrained optimization problem. But they needed an already convex initial grid to guarantee the convergence of their method. To overcome the latter problem, Tinoco and Barrera [4] presented in 1997 a new quasi-harmonic parameter-dependent functional, and used a continuation approach that permitted to avoid the need for an initial convex grid. Later, in 2003, Barrera and Dom´ınguez-Mota proposed a general regularization procedure for a family of area functionals ([2], [12]) which was suitable for the development of numerically stable continuation approaches. The procedure was shown to be very effective on highly irregular regions. In this paper, we present the use of a regularization scheme, that follows the ideas of Barrera and Dom´ınguez-Mota applied to the harmonic functional developed by Ivanenko and Charakhch’yan. Also, we consider their adaptive discrete functional [9] and derive a matrix form for it. We present a practical implementation of our proposal of a continuation algorithm, to generate convex grids using our regularized harmonic and adaptive functionals. The last one is tested using monitor functions that represent a given surface to show that its performance is the one expected. 2
The paper is organized as follows: The continuous formulation of the general grid generation problem is made in Section 2. Sections 3 and 4 describe the discrete approach and some discrete functionals. In section 5 we introduce the new functional Hω and its properties regarding convergence are presented in section 6. Ivanenko’s adaptive harmonic functional is described in section 7 and also our matrix formulation of the continuous functional. In section 8 the adaptive functional HAd,w is presented. Implementation details and the algorithm are described in section 9, and the results of the numerical tests in sections 10 and 11. As expected, conclusions about the work end this paper.
2
Problem Formulation
The regions to be considered for the grid generation problem of our interest, are simple connected domains Ω in the plane whose boundaries are closed polygonal positively oriented Jordan curves. The grid generation problem can be described as the construction of continuous almost everywhere smooth functions x(ξ, η), y(ξ, η) that define a one-to-one mapping x : R 7→ Ω with x = (x(ξ, η), y(ξ, η)) from the unit square R = {(ξ, η)|0 ≤ ξ ≤ 1, 0 ≤ η ≤ 1} onto the physical region Ω to be gridded (Fig. 1). The mapping x(ξ, η) is required to have a full rank Jacobian matrix of positive determinant to preserve the orientation on R and Ω defined by the boundaries. Thus, the continuous grid generation problem is posed in the following way: Problem 1. To find a one to one smooth mapping (or continuous grid) x from the unit square onto the domain Ω that satisfies J = xξ yη − xη yξ > 0
(1)
on the whole unit square. A grid is admissible if x(∂R) = ∂Ω.
3
Fig. 1. Mapping between R and Ω.
In order to solve this problem, Winslow [20] developed a method based on the calculation of the component functions ξ(x, y), η(x, y) of the inverse mapping x−1 defined over Ω and that satisfy the set of Laplace equations ξxx + ξyy = 0, ηxx + ηyy = 0
(2)
with Dirichlet boundary conditions given by x−1 (∂Ω) = ∂R. The usual approach is to make a change of variables to transform (2) into a couple of equations for x(ξ, η) and y(ξ, η): αxξξ − 2βxξη + γxηη = 0 αyξξ − 2βyξη + γyηη = 0
(3)
where α = x2η + yη2 , β = xξ xη + yξ yη , γ = x2ξ + yξ2 . Nonetheless, it is important to notice that the numerical solution of equations (3) using a finite difference scheme does not guarantee positiveness of the Jacobian (i.e., the convexity) of the transformation (ξ, η) 7→ (x, y) either ([9], [16]). A variational formulation for problem 1, related to the set of equations (3) for the inverse mapping x−1 , was first considered in 1982 by Brackbill and Saltzman [5]. They proposed the functional IS =
Z Z
Ω
[(∇ξ)2 + (∇η)2 ] dxdy =
Z Z
x2ξ + x2η + yξ2 + yη2 dξ dη. J R
(4)
Equations (3) are the Euler-Lagrange expressions for the functional (4).
3
Discrete Grid Formulation
One of the most appealing features of Ivanenko’s Direct Optimization Method is the simplicity of its theoretical framework based on a few geometrical and 4
intuitive concepts. This will be clear in the next sections, where we will pose a minimal set of requirements in order to propose an efficient solution to the problem we are interested in. Let us consider a region Ω of the plane, defined by a simple, closed and counterclockwise oriented polygonal curve γ of vertices V = {v1 , v2 , · · · , vq } (Fig. 2). V11 V12
V10
V9
V13
V8
v14 Ω
V7 ∂Ω V6
V1 V2 V4
V5
V3
Fig. 2. Example of a region defined by a simple closed polygonal curve.
Definition. Let m, n be natural numbers with m, n > 2. A set G of points of the plane with boundaries L1 (G) = {Pi,1 |i = 1, ..., m} L2 (G) = {Pm,j |j = 1, ..., n} L3 (G) = {Pi,n |i = 1, ..., m} L4 (G) = {P1,j |j = 1, ..., n} is called a structured, admissible and discrete grid V ⊆
4 [
1
for Ω, of order m × n, if
Li (G).
i=1
Besides, we will say that G is convex if each one of the (m − 1)(n − 1) quadrilaterals (or cells) ci,j of vertices {Pi,j , Pi+1,j , Pi,j+1, Pi+1,j+1}, with 1 ≤ i < m y 1 ≤ j < n, is convex. The sets L1 (G), L2 (G), L3 (G) y L4 (G) will be considered the sides of the grid boundary or the grid sides, and appear in the definition to emphasize our interest in having the same boundary for the region and for the grid. 1
With quadrilateral elements.
5
4
Discrete functionals
Ivanenko and Charakhch’yan’s used quadrilateral isoperimetric finite elements to discretize the continuous harmonic functional (4) and obtained H(G), the Discrete Harmonic Functional: N λ(△ki,j ) X λ(△q ) H(G) = = k q=1 α(△q ) i=1 j=1 k=1 α(△i,j ) m X n X 4 X
(5)
where m and n are the number of horizontal and vertical points of a discrete grid G respectively, k is used to sequentially denote the four triangles formed with the vertices of each quadrilateral grid cell, λ(△q ) is the length functional and α(△q ) is twice the oriented area for such triangles, as described below. In this formula, the variable N stands for the total number of triangles, equal to four times the number of grid cells. This is the basic functional on which we are going to develop the new ones. Hereinafter, it will be of the greatest relevance to consider every grid cell as S
S
R
(4)
(3)
∆
∆
(1)
(2)
∆
∆
Q
P
R
P
Q
Fig. 3. The four oriented triangles defined by a quadrilateral grid cell.
divided into four oriented triangles in order to control the convexity of the cells (Fig. 3). The orientation of the boundary induces an orientation of the cells and triangles of the cells, that induces a sign on the triangle areas which is the key to solve the grid generation problem as an unconstrained optimization one. Thus, for the oriented triangle with vertices Q, P, R ∈ IR2 , the functions λ and α are defined as λ(△(Q, P, R)) = kP − Qk2 + kP − Rk2
α(△(Q, P, R)) = (P − Q)t J2 (P − R)
where k · k denotes the Euclidean norm and J2 is the matrix
J2 =
0 1 . −1 0 6
(6)
A grid G is convex if min{α(△q ) > 0|q = 1, ..., N}. In what follows we will use α and α(△) indistinctly. In this context, the discrete grid generation problem can be posed, in general, as a large scale optimization one: Problem 2. To solve G∗ = arg min G
N X
f (△q )
q=1
over the set of admissible grids for a region. The adequate selection of f is the keystone to get some important properties on those minimizers, such as convexity and smoothness. It is even possible to assert that for numerical purposes, the main goal is to propose functions f for which such minimizers are convex grids. Thus, we will also refer to a numerical grid that minimizes some discrete functional as an optimal grid. For the particular case of IS , the minimizers are approximated by minimizing H(G), and it is possible to show that it has at least a minimizer in the set of all the convex grids for Ω [13]. However, the procedure has the drawback that an initial convex grid is required for the optimization. Tinoco and Barrera [4] proposed a new quasi-harmonic functional that does not require an initial convex grid: Fω (G) =
N X
λ(△q ) q=1 ω + α(△q )
(7)
where ω > 0 is chosen in such a way that in the optimization process, the value in the denominator of (7) in always positive. They also used a continuation approach and proved, under minimal requirements, the existence of optimal convex grids for Fω on plane irregular regions and that for ω sufficiently small, those grids are essentially harmonic. Most of the times, the functional Fω is very good for irregular regions. This means that is very fast due to the update of ω of the form ω ← k1 α− (G) + k2 for appropriate values of k1 > 1 and k2 ≪ 1. However, for quite difficult irregular regions, this strategy may fail, making necessary to find, frequently by trial and error, new values of k1 and k2 to succeed. An example of this fact and grid comparisons between Fω and Hω are given in section 12. 7
Another drawback of Fω is the discontinuity due to the pole at −ω. This can be avoided with some kind of regularization of Fω to guarantee continuity and differentiability. A natural question that arises immediately is whether it is possible to develop an alternative theoretical framework for C 2 functions to generate convex grids. In this direction, Barrera and Dom´ınguez-Mota made a comprehensive analysis of the properties of a family of continuous discrete functionals of area (i.e., with no explicit dependence on λ) to solve the convex grid generation problem ([2], [12]) and proved that if f is a C 2 strictly decreasing convex and bounded below function such that f (α) → 0 as α → ∞, then Sω (G) =
N X
f (ωα(△q ))
(8)
q=1
considered as the objective function in the optimization problem 2, is minimized by convex grids for ω > 0 large enough [2]. It is very important to emphasize that, as a straight consequence, for numerical purposes some very economical choices for f like
f (α) =
1/α,
α≥1
(α − 1)(α − 2) + 1,
α 0, with the same optimal grids as H(G), for ω small enough. To accomplish this goal, the main idea is to replace the area function α−1 in H(G) by a function capable of acting as a continuous barrier and defining f (△q ) = λ(△q )ϕω (α). A natural choice for ϕω (α) is a continuous, convex and strictly decreasing C 2 function identical to 1/α for α ≥ ω (see [2] for details). For instance, we can consider for a positive real number ω
ϕω (α) =
(2ω
− α)/ω 2 , α < ω
1/α,
(9)
α ≥ ω. 8
Thus, we propose the functional Hω defined by Hω (G) =
N X
λ(△q )ϕω (△q ),
(10)
q=1
which is always positive.
6
Properties of Hω
Our immediate goal is to show that the minimizers of Hω (G) are harmonic grids for ω sufficiently small. In our results we are going to use the minimum and average oriented triangle area: α (G) = min {α(△q )}, 1≤q≤N
α(G) = 2A/N, respectively, where A is the area of the region Ω. If α (G) > 0, then G is a convex grid. We will also denote the set of admissible grids on Ω as M(Ω) and the set of convex grids as M0 (Ω) = {G ∈ Ω|α (G) > 0}. The following theorem provides upper and lower bounds for the length functional on a grid. Theorem 1 Let Ω be a polygonal region in the plane such that M0 (Ω) 6= φ. Then there exists a positive real number K such that for every G in M0 (Ω) L(G) =
N X
q=1
λ(△q ) ≤ K
(11)
and 4A/K ≤ 1, where A is the positive oriented area of Ω. PROOF. Let us note that if G is a convex grid, then the vertices of every cell in G lie within the region defined by the boundary of Ω. Hence, the first 9
inequality (11) is true since the continuous function L(G) is evaluated in the compact closure of Ω whenever G ∈ M0 (Ω). The second is a consequence of the definitions (6) and the orthogonality of the matrix J2 when applied to the △(S, P, R) to transform 0 ≤ k(R − P ) ± J2 (S − P ))k22 into the equivalent form 2|α(△(S, P, R))| ≤ λ(△(S, P, R)) and therefore
4A = 2(2A) = 2
N X
q=1
α(△q ) ≤
N X
q=1
2α(△q ) ≤
N X
q=1
λ(△q ) ≤ K.
2
The next theorem provides a lower bound on α (G∗ ) for the minimizers G∗ of Hω (G) for ω fixed. Theorem 2 Let Ω be a polygonal region on the plane with positive oriented area A for which M0 (Ω) 6= φ. Let G0 ∈ M0 (Ω) and ω such that 0 < ω ≤ α (G0 ), then for G∗ = arg min {Hω (G)} G∈M (Ω)
√ the inequality −2 Aω < α (G∗ ) holds. ˆ q the triangles in the cells of the grids PROOF. Let us denote as △q and △ ˆ a triangle of minimum oriented area of G∗ . G0 and G∗ respectively, and by △ n By choosing ω such that 0 ω. Also the functionals 13
Hω (G) and HAd,ω (G) share similar bounds and properties regarding convergence. Thus, for the sake of brevity, the analogous proofs of theorems 1 and 2 for the functional HAd,ω (G) will not be included.
9
Implementation details
Theorem 3 can be readily seen as a procedure to generate a grid sequence that converges to a harmonic one. The details of our implementation are the following: (1) Convex grid. For a practical implementation the definition of convexity, i.e. α (G∗n ) > 0, should be tested using a tolerance value ǫc and we will say that a grid is ǫc -convex if α (G∗n ) > ǫc . To make it scale invariant we set ǫc = ǫc · α(G). Thus, an ǫc -convex grid satisfies that α (G∗n ) > ǫc · α(G) = ǫc . NOTE: A large value of ǫc is sometimes the reason that the procedure do not converge, but might do for a lower one. (2) Initial value for ω1 . For an initial non convex grid, a reasonable choice for this value is half the average area: ω1 = α(G)/2. If the initial grid is already convex, it is convenient to set ω1 = ǫc . (3) Optimization loop. In order to optimize the grid for each value of ωn , the tolerances for the stopping criteria, tolf for the relative error between successive function values and tolg for the gradient norm, should be set with care, see [17] for details. In the implementation presented below we write the values we proved in our codes. (4) Refinement. Once a convex grid is obtained, the optimization should be executed again with ω = ωmin , where ωmin ≤ ǫc , in order to guarantee that the minimizer is also harmonic. Thus, the algorithm for the generation of harmonic grids using the functional Hω and HAd,ω follows:
Algorithm 1. (Adaptive) Harmonic grid generation. Step 0. Given ǫc , ωmin and an initial grid G1 . Set ǫc = ǫc · α(G), n = 1 if (α (G∗n ) < ǫc ) then (non-convex) Set ω1 = α(G)/2, tolf= 10−4 , tolg= 10−2 else (convex) Set ω1 = ωmin, tolf= (epsmach)1/2 , tolg= 10−5 end 14
Step 1. Solve G∗n = arg min {Hωn (G)} G∈M (Ω)
Step 2.
or
G∗n = arg min {HAd,ωn (G)} G∈M (Ω)
if (α (G∗n ) > ǫc ) then (finish) ‘‘an ǫc -convex grid has been obtained for the current value of ǫc ” else (non-convex) Set ωn+1 = max{ωn /2, ωmin} if(ωn+1 = ωmin ) (finish)‘‘NO ǫc -convex grid has been obtained Try with a lower value for ǫc " else Set tolf= max{(epsmach)1/2 , tolf · 10−1 }, tolg= max{10−5, tolf · 10−1 }, n = n + 1 go back to step 1 end end
It is important to note that if the initial grid for HAd,ω is not convex, then one might begin using Hω and once the convex harmonic grid is obtained, continue minimizing HAd,ω .
10
Numerical tests for Hω
The algorithm of the previous section was applied to four regions defined by quite irregular boundaries, three of them representing the shapes of Havana bay, lake Ucha in Russia and Peru, and the fourth one is just a strongly nonconvex region that we will call the Swan. For the examples presented, a 40 × 40 grid was generated by transfinite interpolation as the initial estimation. The optimization was carried out using a Trust Region Newton Method (TRON) that solves bound constrained nonlinear optimization problems [17]. The final harmonic discrete grids are shown in figure 5. Even though each one of the selected boundaries has a great inherent complexity, the harmonic grids were obtained in at most 5 updates of ω.In order to show this last assertion, we present in figure 6 the grids obtained at the end of the optimization of the harmonic functional at the first four values for the parameter ω in the Swan region.
15
Havana Bay.
Peru.
Ucha lake.
Swan.
Fig. 5. Harmonic grids generated with Hω (G)
As it is well known (see [14], [9]) , the harmonic functional produces coordinate lines moving away from convex portions of the domain. This can be easily controlled using a convex combination of the Harmonic functional with the classical discrete area functional, i.e.
A(G) =
N X
[α(∆k )]2 ,
(14)
k=1
that allows us to control the cell areas without losing the convexity of the optimal grids. Such control is possible due to the continuity of A(G), even though there is no continuation approach for it (Some other useful combinations are addressed in [2]). An example of the effect of this combination of functionals is presented in fig. 7, just taking σ = 0.3 for the convex combination F (G) = σHω (G) + (1 − σ)A(G)
(15)
16
Initial grid
ω=1
ω =0.5
ω =0.25
Fig. 6. Optimal grids for the Region Swan generated by optimizing Hω (G) for the ω values shown. The optimal convex grid for ω =0.125 was already shown in fig. 5
11
Numerical tests for HAd,ω
The algorithm for HAd,ω (G) was used to generate adaptive grids using the following expressions for the monitor functions g(x, y): a) g1 (x, y) = exp(−2(4(2x − 1)2 + 9(2y − 1)2 − 1)2 ). b) g2 (x, y) = sin(2π(x − 2y)). Adaptive grids generated for the already convex grid of fig. 5 obtained with the quasi-harmonic functional Hω (G), for monitor functions g1 and g2 , are presented in figures 8 and 9 respectively.
12
Comparison between Fω and Hω
In order to complete our analysis, the remaining task is to compare the performance of Fω and Hω . The algorithm for the continuation approach for Fω as described in [4] is the following: 17
Fig. 7. Grid for region Swan, using the convex combination of Hω and A(G), with σ =0.3
Algorithm 2. (1) Generate an initial grid G0 (2) Let ω= (3) Solve
−1.05α −α
(G0 ),
α (G0 ) ≤ α(G0 )
(G0 ) + 0.01α(G0 ), otherwise
G∗ = arg min {Fω (G)} G∈M (Ω)
(4) Taking G∗ as a new initial grid, repeat steps 2 and 3 until a convex grid is obtained. To produce the test, we use the regions presented in previous sections, as well as two extra complex regions, Great Britain and Mexico, shown in figure 10. Initial grids were generated by Transfinite Interpolation with different numbers of points by side. The time required to generated a convex grid, usc processor running at ing a personal computer with 512Mb RAM and Intel 1.73GHz, is shown in table 1. These results are quite self explanatory. In general, the algorithm for Fω is very 18
Havana Bay.
Peru.
Ucha lake.
Swan.
Fig. 8. Adaptive grids generated with g1 = exp(−2(4(2x − 1)2 + 9(2y − 1)2 − 1)2 )
HAd,ω (G),
monitor
function
fast, but the step 2 in algorithm 2 is not suitable to guarantee that convexity can be achieved for regions where max{α (G)|G ∈ M(Ω)} is very close to zero. The regions Great Britain and Mexico have this characteristic, causing the algorithm to fail in a number of the test cases. On the other hand, even though the running time was greater, the algorithm for Hω is safer to generate the sought convex grids in all the tests. We are testing a much faster version that will be reported in a future article. One example of harmonic grids with 60 points per side is shown in figure 11.
19
Havana Bay.
Peru.
Ucha lake. Fig. 9. Adaptive g2 = sin(2π(x − 2y))
grids
Swan. generated
with
HAd,ω (G),
Great Britain
monitor
function
Mexico
Fig. 10. Two examples of highly complex regions for the grid generation problem
13
Conclusions and future developments
The cornerstone of an efficient and robust method to generate adaptive harmonic grids was presented. The high reliability of the proposed approach is based on the fact that our functionals HAd,ω and Hω use a very effective general regularization in terms of α, that makes it possible to optimize on the 20
Grid
Hω
Fω
Swan 40x40
0:28.04
0:12.26
Havana 40x40
0:10.02
0:10.83
Ucha 40x40
0:7.47
0:5.62
Swan 60x60
0:46.83
0:30.34
Mexico 60x60
7:50.14
FGCG
Mexico 40x40
1:43.08
FGCG
Havana 60x60
0:32.34
0:28.19
Mexico 30x30
0:33.86
FGCG
Great Britain 40x40
0:22.49
0:13.69
Great Britain 80x80
4:57.46
FGCG
Great Britain 70x70
2:53.29
FGCG
Great Britain 60x60
1:37.50
FGCG
Havana 70x70
0:50.12
0:40.48
Swan 80x80
5:54.35
1:47.00
Ucha 80x80
2:37.08
1:34.50
Table 1 Performance comparison on time between Hω and Fω . Times are measured in minutes. F.G.C.G. stands for “Fail to Get a Convex Grid”.
Great Britain
Mexico
Fig. 11. Examples of harmonic grids generated with Hω for two highly complex regions
whole set of admissible grids for a region, and allows us in consequence to define a safe and robust algorithm. Examples of the numerical solution of partial differential equations using the generated grids are reserved for a future article.
21
Acknowledgements We want to thank to Intercambio Acad´emico of UNAM and the “Macroproyecto: Tecnolog´ıas para la Universidad de la Informaci´on y la Computaci´on”, CIC-UMSNH Grant 9.16 for the financial supporting for this work, and also the Institute of Cybernetics, Mathematics and Physics of Cuba. Many thanks due to the reviewers of this paper for their valuable suggestions and comments.
References [1] C. Armstrong, S. Bridgett, R. Donaghy, T. Li, C. McGleenan, R. McKeag and D. Robinson, Medials for Meshing and More, Proceedings of 4th International Meshing Roundtable, Sandia National Laboratories (1995), p.277-288. [2] P. Barrera-S´ anchez, F.J. Dom´ınguez-Mota and G.F. Gonz´ alez-Flores, Robust Discrete Grid Generation on Plane Irregular Regions, Computational Mathematics and Mathematical Physics, Vol.43, No. 6, (2003) pp. 845-854. [3] P. Barrera-S´ anchez, L. Castellanos and A. P´erez Dom´ınguez, M´etodos variacionales discretos para la generaci´ on de mallas, DGAPA-UNAM, M´exico, (1994). [4] P. Barrera-S´ anchez and J.G. Tinoco Ruiz, Smooth and convex grid generation over general plane regions, Mathematics and computers in simulation 46, (1998) pp. 87-102 . [5] J.U. Brackbill and J.S. Saltzman, Adaptive zoning for singular problems in two dimensions, J. Comput. Phys., 46(3), (1982) p. 342-368. [6] J.E. Castillo, Variational Grid Generation. Ph.D. Thesis, University of New Mexico, (1987). [7] J.E. Castillo, Mathematical Aspects of Numerical Grid Generation. SIAM, (1991). [8] A.A. Charakhch’yan and S.A. Ivanenko, Curvilinear Grids of Convex Quadrilaterals, Comput. Maths. Math. Phys. 28, No.2, (1988) 126-133 . [9] A.A. Charakhch’yan and S.A. Ivanenko, A variational form of the Winslow grid generator, Journal of Computational Physics 135, No.5, (1997) 385-398. [10] V.F. De Almeida, Domain Deformation Mapping: Application to Variational Mesh Generation, SIAM J. Sci. Comput., 20, No. 4, (1999) pp. 1252-1275. [11] A. de la Cruz Uribe, Generaci´ on Num´erica de Mallas Arm´ onicas-Adaptivas y su Aplicaci´ on a la soluci´ on de algunas EDP’s (in spanish) M.Sc. Thesis. Facultad de Ciencias, U.N.A.M, (2006).
22
[12] F.J. Dom´ınguez-Mota, Sobre la generaci´ on variacional discreta de mallas casiortogonales en el plano (in spanish), Ph.D. Thesis., Facultad de Ciencias, U.N.A.M, (2005). [13] S.A. Ivanenko, Adaptive Grids and Grids on Surfaces, Comput. Maths. Math. Phys. 33, No.9, (1993) 1179-1193. [14] S.A. Ivanenko, Harmonic Mappings, Handbook of Grid Generation, CRC Press, pp. 8.1-8.41, (1999). [15] P. Knupp and N. Rovideaux, A Framework for Variational Grid Generation: Conditioning the Jacobian Matrix with Matrix Norms, SIAM J. Sci. Comput., 21, No. 6, (1999) pp. 2029-2047. [16] P. Knupp and S. Steinberg, Fundamentals of grid generation. CRC Press, (1992). [17] C.J. Lin and J. J. Mor´e, Newton’s method for large-scale bound constrained optimization problems, SIAM Journal on Opt. 4, (1999) pp. 1100-1127. [18] P. Sampl, Semi-structured mesh generation based on medial axis, Proceedings of 9th International Meshing Roundtable, Sandia National Laboratories (2000), p.21-32. [19] S. Steinberg and P.J. Roache, Variational Grid Generation, Numerical Methods for P.D.E.s 2, (1992) pp. 71-96. [20] A.M. Winslow, Numerical solution of quasilinear Poisson equation in nonuniform triangle mesh, J. Comput. Phys. 1, (1967) 149. [21] UNAMALLA. An automatic package for numerical grid generation. Web site: http://www.matematicas.unam.mx/unamalla
23
Facultad de Ciencias, Mexico City, Mexico. b ICIMAF,
c U.M.S.N.H.,
Havana City, Cuba.
Facultad de Ciencias F´ısico-Matem´ aticas, Morelia, Mexico.
Abstract We present a new adaptive-harmonic structured grid generation method. It is based on a functional that shares a common set of minimizers with Ivanenko’s harmonic functional [9]. An unconstrained optimization process related to a continuation parameter is used to guarantee the convexity of the grid cells. Several numerical examples of grids generated on quite irregular 2D regions, show the effective performance of the proposed method and its robustness when preset monitor functions are given. Key words: Variational grid generation, adaptive grids, quasi-harmonic grids, structured grids.
1
Introduction
Variational grid generation has been used quite successfully in the last years. It started with the work of Brackbill and Saltzman in 1982 [5], and the contributions of Steinberg and Roache in 1986 [19], who developed the basis for the continuous problem. Afterwards, the discrete version of the problem was first developed by Castillo [6].
∗ Corresponding author. Email addresses: [email protected] (Pablo Barrera-S´ anchez), [email protected] (Longina Castellanos Noda), [email protected] (Francisco J. Dom´ınguez-Mota), [email protected] (Guilmer F. Gonz´ alez Flores), [email protected] (Angel P´erez Dom´ınguez).
Preprint submitted to Mathematics and Computers in Simulations 24 April 2007
In the classical variational setting, the strategy was to obtain the optimal grid by solving -numerically in many cases- the Euler-Lagrange equations associated to a functional. There is a wealth of literature available in the subject, mainly devoted to the study of the theoretical properties of the functionals in order to guarantee that their minimizers are convex grids. It has been successfully applied to many regions with a number of symmetries or with an elongated geometry. Comprehensive treatments of the problem can be found in Castillo [6], de Almeida [10], and Knupp and Robidoux [15]. There are also some other approaches to the grid generation problem, for instance, by subdividing the region into simpler subregions, which transforms the problem of gridding an irregular complex domain into the efficient computational choice of a useful subdivision [1,18]. However, from our point of view, the structured gridding of highly irregular planar regions has not been satisfactorily explored. The need of solving modeling problems in domains as realistic as possible was indeed what led to Ivanenko and Charakhch’yan [13] in 1988 to solve the grid generation problem, using the discretization of the variational principle rather than the Euler Lagrange equations. This allows to control directly the interior nodes of the grid and obtain a solution through a large scale optimization method, known as the Direct Optimization Method [16]. Besides, in order to be able to work with the discretized functionals it was necessary to transform the requirement of the positiveness of the Jacobian of the mapping between the regions, that guarantees the convexity of the grid. They proposed a very simple check [8,9,13,14]: the convexity of each grid cell, which represents a major advance from the practical point of view. Ivanenko and Charakhch’yan’s method allowed to generate smooth grids with a great number of nodes, on very irregular regions, in a very efficient and practical way: by solving a large scale unconstrained optimization problem. But they needed an already convex initial grid to guarantee the convergence of their method. To overcome the latter problem, Tinoco and Barrera [4] presented in 1997 a new quasi-harmonic parameter-dependent functional, and used a continuation approach that permitted to avoid the need for an initial convex grid. Later, in 2003, Barrera and Dom´ınguez-Mota proposed a general regularization procedure for a family of area functionals ([2], [12]) which was suitable for the development of numerically stable continuation approaches. The procedure was shown to be very effective on highly irregular regions. In this paper, we present the use of a regularization scheme, that follows the ideas of Barrera and Dom´ınguez-Mota applied to the harmonic functional developed by Ivanenko and Charakhch’yan. Also, we consider their adaptive discrete functional [9] and derive a matrix form for it. We present a practical implementation of our proposal of a continuation algorithm, to generate convex grids using our regularized harmonic and adaptive functionals. The last one is tested using monitor functions that represent a given surface to show that its performance is the one expected. 2
The paper is organized as follows: The continuous formulation of the general grid generation problem is made in Section 2. Sections 3 and 4 describe the discrete approach and some discrete functionals. In section 5 we introduce the new functional Hω and its properties regarding convergence are presented in section 6. Ivanenko’s adaptive harmonic functional is described in section 7 and also our matrix formulation of the continuous functional. In section 8 the adaptive functional HAd,w is presented. Implementation details and the algorithm are described in section 9, and the results of the numerical tests in sections 10 and 11. As expected, conclusions about the work end this paper.
2
Problem Formulation
The regions to be considered for the grid generation problem of our interest, are simple connected domains Ω in the plane whose boundaries are closed polygonal positively oriented Jordan curves. The grid generation problem can be described as the construction of continuous almost everywhere smooth functions x(ξ, η), y(ξ, η) that define a one-to-one mapping x : R 7→ Ω with x = (x(ξ, η), y(ξ, η)) from the unit square R = {(ξ, η)|0 ≤ ξ ≤ 1, 0 ≤ η ≤ 1} onto the physical region Ω to be gridded (Fig. 1). The mapping x(ξ, η) is required to have a full rank Jacobian matrix of positive determinant to preserve the orientation on R and Ω defined by the boundaries. Thus, the continuous grid generation problem is posed in the following way: Problem 1. To find a one to one smooth mapping (or continuous grid) x from the unit square onto the domain Ω that satisfies J = xξ yη − xη yξ > 0
(1)
on the whole unit square. A grid is admissible if x(∂R) = ∂Ω.
3
Fig. 1. Mapping between R and Ω.
In order to solve this problem, Winslow [20] developed a method based on the calculation of the component functions ξ(x, y), η(x, y) of the inverse mapping x−1 defined over Ω and that satisfy the set of Laplace equations ξxx + ξyy = 0, ηxx + ηyy = 0
(2)
with Dirichlet boundary conditions given by x−1 (∂Ω) = ∂R. The usual approach is to make a change of variables to transform (2) into a couple of equations for x(ξ, η) and y(ξ, η): αxξξ − 2βxξη + γxηη = 0 αyξξ − 2βyξη + γyηη = 0
(3)
where α = x2η + yη2 , β = xξ xη + yξ yη , γ = x2ξ + yξ2 . Nonetheless, it is important to notice that the numerical solution of equations (3) using a finite difference scheme does not guarantee positiveness of the Jacobian (i.e., the convexity) of the transformation (ξ, η) 7→ (x, y) either ([9], [16]). A variational formulation for problem 1, related to the set of equations (3) for the inverse mapping x−1 , was first considered in 1982 by Brackbill and Saltzman [5]. They proposed the functional IS =
Z Z
Ω
[(∇ξ)2 + (∇η)2 ] dxdy =
Z Z
x2ξ + x2η + yξ2 + yη2 dξ dη. J R
(4)
Equations (3) are the Euler-Lagrange expressions for the functional (4).
3
Discrete Grid Formulation
One of the most appealing features of Ivanenko’s Direct Optimization Method is the simplicity of its theoretical framework based on a few geometrical and 4
intuitive concepts. This will be clear in the next sections, where we will pose a minimal set of requirements in order to propose an efficient solution to the problem we are interested in. Let us consider a region Ω of the plane, defined by a simple, closed and counterclockwise oriented polygonal curve γ of vertices V = {v1 , v2 , · · · , vq } (Fig. 2). V11 V12
V10
V9
V13
V8
v14 Ω
V7 ∂Ω V6
V1 V2 V4
V5
V3
Fig. 2. Example of a region defined by a simple closed polygonal curve.
Definition. Let m, n be natural numbers with m, n > 2. A set G of points of the plane with boundaries L1 (G) = {Pi,1 |i = 1, ..., m} L2 (G) = {Pm,j |j = 1, ..., n} L3 (G) = {Pi,n |i = 1, ..., m} L4 (G) = {P1,j |j = 1, ..., n} is called a structured, admissible and discrete grid V ⊆
4 [
1
for Ω, of order m × n, if
Li (G).
i=1
Besides, we will say that G is convex if each one of the (m − 1)(n − 1) quadrilaterals (or cells) ci,j of vertices {Pi,j , Pi+1,j , Pi,j+1, Pi+1,j+1}, with 1 ≤ i < m y 1 ≤ j < n, is convex. The sets L1 (G), L2 (G), L3 (G) y L4 (G) will be considered the sides of the grid boundary or the grid sides, and appear in the definition to emphasize our interest in having the same boundary for the region and for the grid. 1
With quadrilateral elements.
5
4
Discrete functionals
Ivanenko and Charakhch’yan’s used quadrilateral isoperimetric finite elements to discretize the continuous harmonic functional (4) and obtained H(G), the Discrete Harmonic Functional: N λ(△ki,j ) X λ(△q ) H(G) = = k q=1 α(△q ) i=1 j=1 k=1 α(△i,j ) m X n X 4 X
(5)
where m and n are the number of horizontal and vertical points of a discrete grid G respectively, k is used to sequentially denote the four triangles formed with the vertices of each quadrilateral grid cell, λ(△q ) is the length functional and α(△q ) is twice the oriented area for such triangles, as described below. In this formula, the variable N stands for the total number of triangles, equal to four times the number of grid cells. This is the basic functional on which we are going to develop the new ones. Hereinafter, it will be of the greatest relevance to consider every grid cell as S
S
R
(4)
(3)
∆
∆
(1)
(2)
∆
∆
Q
P
R
P
Q
Fig. 3. The four oriented triangles defined by a quadrilateral grid cell.
divided into four oriented triangles in order to control the convexity of the cells (Fig. 3). The orientation of the boundary induces an orientation of the cells and triangles of the cells, that induces a sign on the triangle areas which is the key to solve the grid generation problem as an unconstrained optimization one. Thus, for the oriented triangle with vertices Q, P, R ∈ IR2 , the functions λ and α are defined as λ(△(Q, P, R)) = kP − Qk2 + kP − Rk2
α(△(Q, P, R)) = (P − Q)t J2 (P − R)
where k · k denotes the Euclidean norm and J2 is the matrix
J2 =
0 1 . −1 0 6
(6)
A grid G is convex if min{α(△q ) > 0|q = 1, ..., N}. In what follows we will use α and α(△) indistinctly. In this context, the discrete grid generation problem can be posed, in general, as a large scale optimization one: Problem 2. To solve G∗ = arg min G
N X
f (△q )
q=1
over the set of admissible grids for a region. The adequate selection of f is the keystone to get some important properties on those minimizers, such as convexity and smoothness. It is even possible to assert that for numerical purposes, the main goal is to propose functions f for which such minimizers are convex grids. Thus, we will also refer to a numerical grid that minimizes some discrete functional as an optimal grid. For the particular case of IS , the minimizers are approximated by minimizing H(G), and it is possible to show that it has at least a minimizer in the set of all the convex grids for Ω [13]. However, the procedure has the drawback that an initial convex grid is required for the optimization. Tinoco and Barrera [4] proposed a new quasi-harmonic functional that does not require an initial convex grid: Fω (G) =
N X
λ(△q ) q=1 ω + α(△q )
(7)
where ω > 0 is chosen in such a way that in the optimization process, the value in the denominator of (7) in always positive. They also used a continuation approach and proved, under minimal requirements, the existence of optimal convex grids for Fω on plane irregular regions and that for ω sufficiently small, those grids are essentially harmonic. Most of the times, the functional Fω is very good for irregular regions. This means that is very fast due to the update of ω of the form ω ← k1 α− (G) + k2 for appropriate values of k1 > 1 and k2 ≪ 1. However, for quite difficult irregular regions, this strategy may fail, making necessary to find, frequently by trial and error, new values of k1 and k2 to succeed. An example of this fact and grid comparisons between Fω and Hω are given in section 12. 7
Another drawback of Fω is the discontinuity due to the pole at −ω. This can be avoided with some kind of regularization of Fω to guarantee continuity and differentiability. A natural question that arises immediately is whether it is possible to develop an alternative theoretical framework for C 2 functions to generate convex grids. In this direction, Barrera and Dom´ınguez-Mota made a comprehensive analysis of the properties of a family of continuous discrete functionals of area (i.e., with no explicit dependence on λ) to solve the convex grid generation problem ([2], [12]) and proved that if f is a C 2 strictly decreasing convex and bounded below function such that f (α) → 0 as α → ∞, then Sω (G) =
N X
f (ωα(△q ))
(8)
q=1
considered as the objective function in the optimization problem 2, is minimized by convex grids for ω > 0 large enough [2]. It is very important to emphasize that, as a straight consequence, for numerical purposes some very economical choices for f like
f (α) =
1/α,
α≥1
(α − 1)(α − 2) + 1,
α 0, with the same optimal grids as H(G), for ω small enough. To accomplish this goal, the main idea is to replace the area function α−1 in H(G) by a function capable of acting as a continuous barrier and defining f (△q ) = λ(△q )ϕω (α). A natural choice for ϕω (α) is a continuous, convex and strictly decreasing C 2 function identical to 1/α for α ≥ ω (see [2] for details). For instance, we can consider for a positive real number ω
ϕω (α) =
(2ω
− α)/ω 2 , α < ω
1/α,
(9)
α ≥ ω. 8
Thus, we propose the functional Hω defined by Hω (G) =
N X
λ(△q )ϕω (△q ),
(10)
q=1
which is always positive.
6
Properties of Hω
Our immediate goal is to show that the minimizers of Hω (G) are harmonic grids for ω sufficiently small. In our results we are going to use the minimum and average oriented triangle area: α (G) = min {α(△q )}, 1≤q≤N
α(G) = 2A/N, respectively, where A is the area of the region Ω. If α (G) > 0, then G is a convex grid. We will also denote the set of admissible grids on Ω as M(Ω) and the set of convex grids as M0 (Ω) = {G ∈ Ω|α (G) > 0}. The following theorem provides upper and lower bounds for the length functional on a grid. Theorem 1 Let Ω be a polygonal region in the plane such that M0 (Ω) 6= φ. Then there exists a positive real number K such that for every G in M0 (Ω) L(G) =
N X
q=1
λ(△q ) ≤ K
(11)
and 4A/K ≤ 1, where A is the positive oriented area of Ω. PROOF. Let us note that if G is a convex grid, then the vertices of every cell in G lie within the region defined by the boundary of Ω. Hence, the first 9
inequality (11) is true since the continuous function L(G) is evaluated in the compact closure of Ω whenever G ∈ M0 (Ω). The second is a consequence of the definitions (6) and the orthogonality of the matrix J2 when applied to the △(S, P, R) to transform 0 ≤ k(R − P ) ± J2 (S − P ))k22 into the equivalent form 2|α(△(S, P, R))| ≤ λ(△(S, P, R)) and therefore
4A = 2(2A) = 2
N X
q=1
α(△q ) ≤
N X
q=1
2α(△q ) ≤
N X
q=1
λ(△q ) ≤ K.
2
The next theorem provides a lower bound on α (G∗ ) for the minimizers G∗ of Hω (G) for ω fixed. Theorem 2 Let Ω be a polygonal region on the plane with positive oriented area A for which M0 (Ω) 6= φ. Let G0 ∈ M0 (Ω) and ω such that 0 < ω ≤ α (G0 ), then for G∗ = arg min {Hω (G)} G∈M (Ω)
√ the inequality −2 Aω < α (G∗ ) holds. ˆ q the triangles in the cells of the grids PROOF. Let us denote as △q and △ ˆ a triangle of minimum oriented area of G∗ . G0 and G∗ respectively, and by △ n By choosing ω such that 0 ω. Also the functionals 13
Hω (G) and HAd,ω (G) share similar bounds and properties regarding convergence. Thus, for the sake of brevity, the analogous proofs of theorems 1 and 2 for the functional HAd,ω (G) will not be included.
9
Implementation details
Theorem 3 can be readily seen as a procedure to generate a grid sequence that converges to a harmonic one. The details of our implementation are the following: (1) Convex grid. For a practical implementation the definition of convexity, i.e. α (G∗n ) > 0, should be tested using a tolerance value ǫc and we will say that a grid is ǫc -convex if α (G∗n ) > ǫc . To make it scale invariant we set ǫc = ǫc · α(G). Thus, an ǫc -convex grid satisfies that α (G∗n ) > ǫc · α(G) = ǫc . NOTE: A large value of ǫc is sometimes the reason that the procedure do not converge, but might do for a lower one. (2) Initial value for ω1 . For an initial non convex grid, a reasonable choice for this value is half the average area: ω1 = α(G)/2. If the initial grid is already convex, it is convenient to set ω1 = ǫc . (3) Optimization loop. In order to optimize the grid for each value of ωn , the tolerances for the stopping criteria, tolf for the relative error between successive function values and tolg for the gradient norm, should be set with care, see [17] for details. In the implementation presented below we write the values we proved in our codes. (4) Refinement. Once a convex grid is obtained, the optimization should be executed again with ω = ωmin , where ωmin ≤ ǫc , in order to guarantee that the minimizer is also harmonic. Thus, the algorithm for the generation of harmonic grids using the functional Hω and HAd,ω follows:
Algorithm 1. (Adaptive) Harmonic grid generation. Step 0. Given ǫc , ωmin and an initial grid G1 . Set ǫc = ǫc · α(G), n = 1 if (α (G∗n ) < ǫc ) then (non-convex) Set ω1 = α(G)/2, tolf= 10−4 , tolg= 10−2 else (convex) Set ω1 = ωmin, tolf= (epsmach)1/2 , tolg= 10−5 end 14
Step 1. Solve G∗n = arg min {Hωn (G)} G∈M (Ω)
Step 2.
or
G∗n = arg min {HAd,ωn (G)} G∈M (Ω)
if (α (G∗n ) > ǫc ) then (finish) ‘‘an ǫc -convex grid has been obtained for the current value of ǫc ” else (non-convex) Set ωn+1 = max{ωn /2, ωmin} if(ωn+1 = ωmin ) (finish)‘‘NO ǫc -convex grid has been obtained Try with a lower value for ǫc " else Set tolf= max{(epsmach)1/2 , tolf · 10−1 }, tolg= max{10−5, tolf · 10−1 }, n = n + 1 go back to step 1 end end
It is important to note that if the initial grid for HAd,ω is not convex, then one might begin using Hω and once the convex harmonic grid is obtained, continue minimizing HAd,ω .
10
Numerical tests for Hω
The algorithm of the previous section was applied to four regions defined by quite irregular boundaries, three of them representing the shapes of Havana bay, lake Ucha in Russia and Peru, and the fourth one is just a strongly nonconvex region that we will call the Swan. For the examples presented, a 40 × 40 grid was generated by transfinite interpolation as the initial estimation. The optimization was carried out using a Trust Region Newton Method (TRON) that solves bound constrained nonlinear optimization problems [17]. The final harmonic discrete grids are shown in figure 5. Even though each one of the selected boundaries has a great inherent complexity, the harmonic grids were obtained in at most 5 updates of ω.In order to show this last assertion, we present in figure 6 the grids obtained at the end of the optimization of the harmonic functional at the first four values for the parameter ω in the Swan region.
15
Havana Bay.
Peru.
Ucha lake.
Swan.
Fig. 5. Harmonic grids generated with Hω (G)
As it is well known (see [14], [9]) , the harmonic functional produces coordinate lines moving away from convex portions of the domain. This can be easily controlled using a convex combination of the Harmonic functional with the classical discrete area functional, i.e.
A(G) =
N X
[α(∆k )]2 ,
(14)
k=1
that allows us to control the cell areas without losing the convexity of the optimal grids. Such control is possible due to the continuity of A(G), even though there is no continuation approach for it (Some other useful combinations are addressed in [2]). An example of the effect of this combination of functionals is presented in fig. 7, just taking σ = 0.3 for the convex combination F (G) = σHω (G) + (1 − σ)A(G)
(15)
16
Initial grid
ω=1
ω =0.5
ω =0.25
Fig. 6. Optimal grids for the Region Swan generated by optimizing Hω (G) for the ω values shown. The optimal convex grid for ω =0.125 was already shown in fig. 5
11
Numerical tests for HAd,ω
The algorithm for HAd,ω (G) was used to generate adaptive grids using the following expressions for the monitor functions g(x, y): a) g1 (x, y) = exp(−2(4(2x − 1)2 + 9(2y − 1)2 − 1)2 ). b) g2 (x, y) = sin(2π(x − 2y)). Adaptive grids generated for the already convex grid of fig. 5 obtained with the quasi-harmonic functional Hω (G), for monitor functions g1 and g2 , are presented in figures 8 and 9 respectively.
12
Comparison between Fω and Hω
In order to complete our analysis, the remaining task is to compare the performance of Fω and Hω . The algorithm for the continuation approach for Fω as described in [4] is the following: 17
Fig. 7. Grid for region Swan, using the convex combination of Hω and A(G), with σ =0.3
Algorithm 2. (1) Generate an initial grid G0 (2) Let ω= (3) Solve
−1.05α −α
(G0 ),
α (G0 ) ≤ α(G0 )
(G0 ) + 0.01α(G0 ), otherwise
G∗ = arg min {Fω (G)} G∈M (Ω)
(4) Taking G∗ as a new initial grid, repeat steps 2 and 3 until a convex grid is obtained. To produce the test, we use the regions presented in previous sections, as well as two extra complex regions, Great Britain and Mexico, shown in figure 10. Initial grids were generated by Transfinite Interpolation with different numbers of points by side. The time required to generated a convex grid, usc processor running at ing a personal computer with 512Mb RAM and Intel 1.73GHz, is shown in table 1. These results are quite self explanatory. In general, the algorithm for Fω is very 18
Havana Bay.
Peru.
Ucha lake.
Swan.
Fig. 8. Adaptive grids generated with g1 = exp(−2(4(2x − 1)2 + 9(2y − 1)2 − 1)2 )
HAd,ω (G),
monitor
function
fast, but the step 2 in algorithm 2 is not suitable to guarantee that convexity can be achieved for regions where max{α (G)|G ∈ M(Ω)} is very close to zero. The regions Great Britain and Mexico have this characteristic, causing the algorithm to fail in a number of the test cases. On the other hand, even though the running time was greater, the algorithm for Hω is safer to generate the sought convex grids in all the tests. We are testing a much faster version that will be reported in a future article. One example of harmonic grids with 60 points per side is shown in figure 11.
19
Havana Bay.
Peru.
Ucha lake. Fig. 9. Adaptive g2 = sin(2π(x − 2y))
grids
Swan. generated
with
HAd,ω (G),
Great Britain
monitor
function
Mexico
Fig. 10. Two examples of highly complex regions for the grid generation problem
13
Conclusions and future developments
The cornerstone of an efficient and robust method to generate adaptive harmonic grids was presented. The high reliability of the proposed approach is based on the fact that our functionals HAd,ω and Hω use a very effective general regularization in terms of α, that makes it possible to optimize on the 20
Grid
Hω
Fω
Swan 40x40
0:28.04
0:12.26
Havana 40x40
0:10.02
0:10.83
Ucha 40x40
0:7.47
0:5.62
Swan 60x60
0:46.83
0:30.34
Mexico 60x60
7:50.14
FGCG
Mexico 40x40
1:43.08
FGCG
Havana 60x60
0:32.34
0:28.19
Mexico 30x30
0:33.86
FGCG
Great Britain 40x40
0:22.49
0:13.69
Great Britain 80x80
4:57.46
FGCG
Great Britain 70x70
2:53.29
FGCG
Great Britain 60x60
1:37.50
FGCG
Havana 70x70
0:50.12
0:40.48
Swan 80x80
5:54.35
1:47.00
Ucha 80x80
2:37.08
1:34.50
Table 1 Performance comparison on time between Hω and Fω . Times are measured in minutes. F.G.C.G. stands for “Fail to Get a Convex Grid”.
Great Britain
Mexico
Fig. 11. Examples of harmonic grids generated with Hω for two highly complex regions
whole set of admissible grids for a region, and allows us in consequence to define a safe and robust algorithm. Examples of the numerical solution of partial differential equations using the generated grids are reserved for a future article.
21
Acknowledgements We want to thank to Intercambio Acad´emico of UNAM and the “Macroproyecto: Tecnolog´ıas para la Universidad de la Informaci´on y la Computaci´on”, CIC-UMSNH Grant 9.16 for the financial supporting for this work, and also the Institute of Cybernetics, Mathematics and Physics of Cuba. Many thanks due to the reviewers of this paper for their valuable suggestions and comments.
References [1] C. Armstrong, S. Bridgett, R. Donaghy, T. Li, C. McGleenan, R. McKeag and D. Robinson, Medials for Meshing and More, Proceedings of 4th International Meshing Roundtable, Sandia National Laboratories (1995), p.277-288. [2] P. Barrera-S´ anchez, F.J. Dom´ınguez-Mota and G.F. Gonz´ alez-Flores, Robust Discrete Grid Generation on Plane Irregular Regions, Computational Mathematics and Mathematical Physics, Vol.43, No. 6, (2003) pp. 845-854. [3] P. Barrera-S´ anchez, L. Castellanos and A. P´erez Dom´ınguez, M´etodos variacionales discretos para la generaci´ on de mallas, DGAPA-UNAM, M´exico, (1994). [4] P. Barrera-S´ anchez and J.G. Tinoco Ruiz, Smooth and convex grid generation over general plane regions, Mathematics and computers in simulation 46, (1998) pp. 87-102 . [5] J.U. Brackbill and J.S. Saltzman, Adaptive zoning for singular problems in two dimensions, J. Comput. Phys., 46(3), (1982) p. 342-368. [6] J.E. Castillo, Variational Grid Generation. Ph.D. Thesis, University of New Mexico, (1987). [7] J.E. Castillo, Mathematical Aspects of Numerical Grid Generation. SIAM, (1991). [8] A.A. Charakhch’yan and S.A. Ivanenko, Curvilinear Grids of Convex Quadrilaterals, Comput. Maths. Math. Phys. 28, No.2, (1988) 126-133 . [9] A.A. Charakhch’yan and S.A. Ivanenko, A variational form of the Winslow grid generator, Journal of Computational Physics 135, No.5, (1997) 385-398. [10] V.F. De Almeida, Domain Deformation Mapping: Application to Variational Mesh Generation, SIAM J. Sci. Comput., 20, No. 4, (1999) pp. 1252-1275. [11] A. de la Cruz Uribe, Generaci´ on Num´erica de Mallas Arm´ onicas-Adaptivas y su Aplicaci´ on a la soluci´ on de algunas EDP’s (in spanish) M.Sc. Thesis. Facultad de Ciencias, U.N.A.M, (2006).
22
[12] F.J. Dom´ınguez-Mota, Sobre la generaci´ on variacional discreta de mallas casiortogonales en el plano (in spanish), Ph.D. Thesis., Facultad de Ciencias, U.N.A.M, (2005). [13] S.A. Ivanenko, Adaptive Grids and Grids on Surfaces, Comput. Maths. Math. Phys. 33, No.9, (1993) 1179-1193. [14] S.A. Ivanenko, Harmonic Mappings, Handbook of Grid Generation, CRC Press, pp. 8.1-8.41, (1999). [15] P. Knupp and N. Rovideaux, A Framework for Variational Grid Generation: Conditioning the Jacobian Matrix with Matrix Norms, SIAM J. Sci. Comput., 21, No. 6, (1999) pp. 2029-2047. [16] P. Knupp and S. Steinberg, Fundamentals of grid generation. CRC Press, (1992). [17] C.J. Lin and J. J. Mor´e, Newton’s method for large-scale bound constrained optimization problems, SIAM Journal on Opt. 4, (1999) pp. 1100-1127. [18] P. Sampl, Semi-structured mesh generation based on medial axis, Proceedings of 9th International Meshing Roundtable, Sandia National Laboratories (2000), p.21-32. [19] S. Steinberg and P.J. Roache, Variational Grid Generation, Numerical Methods for P.D.E.s 2, (1992) pp. 71-96. [20] A.M. Winslow, Numerical solution of quasilinear Poisson equation in nonuniform triangle mesh, J. Comput. Phys. 1, (1967) 149. [21] UNAMALLA. An automatic package for numerical grid generation. Web site: http://www.matematicas.unam.mx/unamalla
23
