MIMS Technical Report No.00039 (201208171)

A self-organized mesh generator using pattern formation in a reaction-diﬀusion
systemHirofumi Notsua,∗, Daishin Ueyamab , Masahiro Yamaguchic,da Waseda Institute for Advanced Study, Waseda University, 3-4-1, Okubo, Shinjuku-ku, Tokyo 169-8555, Japanb Meiji Institute for Advanced Study of Mathematical Sciences (MIMS), Meiji University, 1-1-1, Higashimita, Tama-ku, Kawasaki, Kanagawa

214-8571, Japanc Graduate School of Advanced Mathematical Sciences, Meiji University, 1-1-1, Higashimita, Tama-ku, Kawasaki, Kanagawa 214-8571, Japand Research Fellow of the Japan Society for the Promotion of Science

Abstract
A new type of mesh generator is developed by using a self-organized pattern in a reaction-diﬀusion system. The
system is the Gray-Scott model, which creates a spot pattern in a specific parameter region. The spots correspond to
nodes of a mesh. The mesh generator has several advantages: the algorithm is simple and processes to improve themesh, such as smoothing, (locally) addition, and removal of nodes, are automatically performed by the system.
Keywords: Mesh generation, Self-organizing pattern formation, Reaction-diﬀusion system, The Gray-Scott model,
The finite diﬀerence method
2000 MSC: 65L50, 35K57, 65L12

1. Introduction
The reaction-diﬀusion system is one of the well studied mathematical models and includes the so-called selfreplication/self-organization mechanism. According to Turing, the combination of diﬀusion and reaction has thepotential to create a pattern in a self-organized manner [1]. Many researchers believe that the mechanism of patternformation in a reaction-diﬀusion system mimics that in nature [2, 3]. Indeed, Turing patterns have been experimentallyconfirmed in chemical systems [4, 5, 6, 7], apart from being observed on the skin of a marine angelfish [8]. The
development of a mathematical theory for the reaction-diﬀusion system in 50 years is stunning. Moreover, many
applications have been developed in the real world [9, 3, 8].
Here, by using the self-organizing pattern formation mechanism, we propose a new type of application of thereaction-diﬀusion system, i.e., a self-organized mesh generator. The proposed mesh generator is a triangular-meshgenerator suitable for numerical simulations, especially for the finite element method (FEM). The strategy used in themesh generator is quite diﬀerent from previous ones, quadtree/octree [10], Delaunay [11, 12, 13], advancing front [14],bubble [15], and so on. In practical numerical simulations, it is important to obtain a good mesh for a specific domain
shape. It is also necessary to optimize the mesh size according to the domain shape and the properties of the solutions.
Our strategy is based on the self-organized pattern appearing in a reaction-diﬀusion system, the Gray-Scott model(GS model) [16, 17, 18, 19], that produces a spot pattern in a suitable parameter region. The spot pattern is obtained
from appropriate initial data, and subsequently, the spots fill the domain automatically. We first introduce the conceptof the self-organized mesh generator based on the GS model. We then demonstrate the advantages and usefulness ofour strategy.

∗ Corresponding author. Tel./fax: +81 3 52862144

Email address: h.notsu@aoni.waseda.jp (Hirofumi Notsu)

Preprint submitted to Applied Mathematics Letters

August 14, 2012

2. Algorithm for a self-organized mesh generator
2.1. The Gray-Scott model
The GS model represents the two reactions U + 2V → 3V and V → P, where U and V are the reacting chemicals
and P is an inert product. The reaction-diﬀusion system, in normalized units, can be written as

∂u



= ∇ · (Du ∇u) − uv2 + F(1 − u),


 ∂t
(GS)



∂v



= ∇ · (Dv ∇v) + uv2 − (F + k)v,∂t
where u = u(x, t) and v = v(x, t) are the concentrations of the chemicals U and V, respectively, as a function of
position x and time t; Du and Dv are the diﬀusion coeﬃcients of the chemicals, respectively, and F and k are positive
parameters. Details of the GS model can be found in references [16, 17, 18, 19]. In a suitable parameter region, theGS model produces a two-dimensional spot pattern. The features of the obtained spot pattern are as follows. (i) The
domain is filled by spots. (ii) For fixed Du and Dv , the distances between neighbouring spots are almost equal. (iii)
The spot pattern exactly fits the domain shape. It is important to note that these features are automatically realized
when the GS model is used.
We use the spot pattern appearing in the GS model for mesh generation. As expected, the spots correspond to thenodes (nodal points) of a triangular mesh. Although Du and Dv in (GS) are constants in the whole domain, we change
the constants into spatially dependent functions, i.e., Du = Du (x) and Dv = Dv (x), to control the local mesh size. We
note that no spot pattern appears in the region {x; Du (x) = Dv (x) = 0}. Finally, a mesh is generated from a set ofnodes by using the Delaunay criterion [20].
2.2. ProcedureLet Ω be a given domain and Γ be the boundary of Ω. For performing computations in Ω, we need a meshe
Th = {Kl }nl=1, where Kl is the so-called element, i.e., a two-dimensional triangle or a three-dimensional tetrahedron,ne is the total number of elements, and h is the representative mesh size. Let n p be the total number of nodes of Th
and Nh (]Nh = n p ) be the set of the nodes. The procedure for the construction of the mesh generator consists of thefollowing five steps and is explained in Figure 1. The five snapshots correspond to the five steps involved in the meshgeneration for a complex domain near Tokyo Bay.

(a)

(b)

(d)

(c)

(e)

Figure 1: Snapshots of the diﬀerent steps (steps 1–5) in the self-organized mesh generation.

2

Step 1: The domain Ω is given by users. In the grey-scale (bitmap) image shown in Figure 1(a), the non-white regionis the domain. We set a fictitious rectangular domain Ω0 (⊃ Ω) and consider heterogeneous diﬀusion coeﬃcients
Du and Dv : Ω0 → R. The colours in the grey-scale image correspond to the values of Dv , which is zero andpositive in the white and non-white regions, respectively.
Step 2: An artificial spot pattern is set as the initial value in the GS model (see Figure 1(b)); the artificial spots are
put on the interfaces between two diﬀerent colours, including the boundary Γ, and the locations of the artificial
spots are fixed.
Step 3: The system (GS) is numerically solved in Ω0 by the finite diﬀerence method (FDM). Subsequently, the artificial fixed spots on the interfaces produce daughter spots until the domain is filled (Figure 1(c)). The distancesof the spots are decided on the basis of the values of the diﬀusion coeﬃcients.
Step 4: A set of nodes Nh is obtained from the final spot pattern (Figure 1(d)).
Step 5: A mesh Th is generated from Nh by using the Delaunay criterion (Figure 1(e)).
The self-organized mesh generator described above has following advantages: (i) The procedure is quite simple.
(ii) Nodes are produced almost automatically by solving the GS model. It does not need complex techniques and criteria, for example, techniques for smoothing, (local) addition, or removal of nodes; such operations are automaticallyrealized.
2.3. Implementation
We describe the details of the implementation in each step of the procedure. First, we fix F = 0.032 and k = 0.063
and obtain a spot pattern for Du /Dv = 2 in the GS model. Since our implementation is realized through (pixel) imagebased computation, we introduce a discrete domain Ω = {xi, j ; (i, j) ∈ S}, where xi, j ≡ (i δx, j δx) for i, j ∈ Z∪{Z+1/2},
δx ≡ 1 is the computational mesh size for the FDM, and S ⊂ Z2 .
In step 1, we set a discrete rectangular domain Ω0 ≡ {xi, j ; (i, j) ∈ S0 } satisfying Ω ⊂ Ω0 (S ⊂ S0 ) with boundaryΓ0 ≡ {xi, j < Ω0 ; {xi, j±1 , xi±1, j } ∩ Ω0 , ∅}. We assume that white colour is used for Ω0 \ Ω and that the colour inΩ is not white. Let Ω̃0 ≡ Ω0 ∪ Γ0 . In the case of Figure 1, Ω0 = {1, · · · , 1000} × {1, · · · , 800}(= S0 ) ⊂ N2 andΓ0 = ({1, · · · , 1000} × {0, 801}) ∪ ({0, 1001} × {1, · · · , 800}). We define spatially dependent diﬀusion coeﬃcients Du
and Dv : Ω̃0 → R by
(Du , Dv )(x) ≡ (2, 1)d(x),(1)
where d : Ω̃0 → R is a function that takes positive values in Ω and zero in Ω̃0 \ Ω. Then, the relation Du /Dv = 2
holds in Ω and a spot pattern appears only in Ω. The domain Ω can be recognized from the function d. In Figure 1,
d = 1.5×10−1 , 9.6×10−2 , 5.4×10−2 , and 0 for the black, dark-grey, light-grey, and white (Ω̃0 \Ω) regions, respectively.
When d is a constant in whole domain, the rough relation√h ≈ c∗ d(2)
with c∗ = 48 holds by our simulations. We can therefore expect the distance between neighbouring spots in the black,dark-grey, and light-grey regions will be about 18.6, 14.9, and 11.2, respectively.
To solve the GS model numerically, we introduce an approximation (U n , V n )(xi, j ) ∈ R2 of (u, v)(xi, j , tn ) for xi, j ∈Ω̃0 and n = 0, 1, · · ·; here, tn ≡ n δt and δt is a time increment. Let V∗ ≡ 0.6. ‘An artificial spot pattern’ in step 2 isU 0 (xi, j ) = 1 (xi, j ∈ Ω̃0 ), and



 V∗
V 0 (xi, j ) =  0(xi, j ∈ A),(xi, j ∈ Ω̃0 \ A),(3)

where A is a set of lattice points on the interfaces between two diﬀerent√ colours. We note that the distance between
neighbouring spots in A does not have to be consistent exactly with c∗ d in (2). Indeed, in Figure 1 there are three
types of interfaces-between black and dark grey, between dark grey and light grey, and between light grey and white-,and the distance between neighbouring spots in A at these interfaces is almost 18, 14, and 10, respectively.
3

In step 3, we use the following explicit finite diﬀerence scheme to compute the GS model: find {(U n , V n )(xi, j ); xi, j ∈Ω̃0 , n = 0, 1, · · ·} ⊂ R2 such that
 n+1
U − Un(xi, j ) = (Lδx (u)U n )(xi, j ) − U n (V n )2 (xi, j ) + F(1 − U n )(xi, j )(xi, j ∈ Ω0 ),



δt



n+1
n
V −V
(4)



(xi, j ) = (Lδx (v)V n )(xi, j ) + U n (V n )2 (xi, j ) − (F + k)V n (xi, j ) (xi, j ∈ Ω0 \ A),



δt




V n+1 (xi, j ) = V∗(xi, j ∈ A),
for n = 0, 1, · · ·, with the initial condition (3) (artificial spot pattern in step 2), and the zero Neumann boundary
condition for Γ0 ; Lδx (a) : {W : Ω̃0 → R} → {W : Ω0 → R} (a = u, v) is an operator defined by(Lδx (a)W)(xi, j ) ≡2∑[∇δx,k {Da (∇δx,k W)}](xi, j ) (xi, j ∈ Ω0 ),k=1

with (∇δx,k W)(xi, j ) ≡ {W(xi, j + (δx/2)ek ) − W(xi, j − (δx/2)ek )}/δx (k = 1, 2), ek ≡ (δk1 , δk2 ). Further, δ pq (p, q = 1, 2)
is the Kronecker delta and Da (xi, j + (δx/2)ek ) ≡ {Da (xi, j + δx ek ) + Da (xi, j )}/2. We can write
[
]
]
1 [
∇δx,1 {Du (∇δx,1 U n )} (xi, j ) = 2 Du (xi+1/2, j ){U n (xi+1, j ) − U n (xi, j )} − Du (xi−1/2, j ){U n (xi, j ) − U n (xi−1, j )} .
δx
For suﬃciently large n0 , we obtain the numerical stationary solution {(U n0 , V n0 )(xi, j )} xi, j ∈Ω0 . In our computation, we
have used the condition
kV n0 − V n0 −1 kl2 /kV n0 −1 kl2 < ε ≡ 10−5∑for every 5, 000 time steps, where for W : Ω0 → R, k·kl2 is a discrete L2 -norm defined by kWkl2 ≡ {δx2 xi, j ∈Ω0 W(xi, j )2 }1/2 .
We note that a highly accurate solution is not required, i.e., a rough computation is enough if spots are created,
and that the mesh size h can be controlled by varying d. Because of Lδx (a) in scheme (4), δt needs to satisfyδt ≤ δx2 /4 max{kDu k∞ , kDv k∞ } (= δx2 /8kdk∞ when relation (1) is used), where k · k∞ is the maximum norm. In
Figure 1, δt = 2/3.
In step 4, Nh is created from the stationary solution {(U n0 , V n0 )(xi, j )} xi, j ∈Ω0 obtained in step 3, i.e., the local peaksof {V n0 (xi, j )} xi, j ∈Ω0 become nodes. In our computation, for xi, j ∈ Ω0 , if the condition{}
V n0 (xi, j ) > max V n0 (xl,m ); xl,m ∈ Ω0 , xl,m , xi, j , l ∈ [i − ci, j , i + ci, j ], m ∈ [ j − ci, j , j + ci, j ](Ci, j )√is satisfied, we define the point xi, j as a node, where ci, j ≡ c0 d(xi, j ), with c0 ≡ 20. Here, we have set ci, j by using (2).
We then obtain Nh , which is defined byNh ≡ {xi, j ∈ Ω0 ; xi, j satisfies (Ci, j )}.
The set of points A used for the ‘artificial spot pattern” in step 1 will be a subset of Nh . This means that we have theflexibility to set the nodes at appropriate points.
In step 5, the Delaunay criterion gives a triangulation Th from Nh , where the problem of non-convexity of thedomain Ω in the triangulation is solved by using d.
3. Quality of meshes
We examine the quality of meshes generated by the self-organized mesh generator for two sample domains. The
quality is confirmed by the distribution of the inner angles and lengths of meshes.
3.1. Example 1
The first sample domain is Figure 1(a) for the generated mesh shown in Figure 1(e), for which n p = 2, 379 andne = 4, 234. Theoretical results for the FEM [21] imply that regular triangles are preferred. Figure 2 shows graphs ofrelative frequency versus angle (left) and length (right). We can see that there are many almost regular triangles andthat there are three peaks of length around 10, 14, and 18. The three peaks correspond to the three positive values ofd. These results imply that the self-organized mesh generator produces a mesh whose elements are almost regular andwhole local mesh sizes are controlled by d.
4

0.3

GS
Relative frequency

Relative frequency

0.15

0.1

0.05

0

0

π/12

π/6

π/4

GS

0.2

0.1

0

π/3 5π/12 π/2 7π/12 2π/3
Angle

0

5

10
15
20
Length of edges

25

30

Figure 2: Graphs of relative frequency vs. angle (left) and length (right) for the first sample domain near Tokyo Bay.

3.2. Example 2
The other sample domain is shown in Figure 3 (leftmost figure); the boundary is complex, the mesh size is uniform,and Ω0 ≡ (1, 250) × (1, 350). In our mesh generator, d = 0.1 and 0 in the black and white regions, respectively, andδt = 1. For this example, we compare our mesh with that generated by FreeFem++ [22], which is one of the most
famous free FEM softwares. The obtained meshes are shown in the centre (FreeFem++, n p = 231, ne = 389) andrightmost panels (ours, n p = 237, ne = 401) of Figure 3, where we have considered identical boundary nodes. Figure 4
shows plots of relative frequency versus angle (left) and length (right). To see the progress of the self-organized meshgenerator, we show two graphs of n = 5, 000 and 20, 000 in addition to the graphs of the FreeFem++ mesh and ourmesh with n = 55, 000 (final). The red and blue lines indicate a mesh generated by FreeFem++ and our (final) mesh,respectively. The green and magenta lines correspond to n = 5, 000 and 20, 000 in step 3. We can observe that the
distribution of the angles and lengths of our meshes improves as n increases, that the inner angles and lengths of ourmeshes are concentrated near π/3 and 15, respectively, and that the final mesh (n = 55, 000) is better than the meshgenerated by FreeFem++.

Figure 3: Second sample domain (left), a mesh generated by FreeFem++ (center), and our mesh (right).

0.1

0.2

0.1

0.05

0

FreeFem++
GS, n=5,000
GS, n=20,000
GS, n=55,000

0.3

0.15

Relative frequency

Relative frequency

0.4

FreeFem++
GS, n=5,000
GS, n=20,000
GS, n=55,000

0.2

0

π/12

π/6

π/4

π/3
Angle

5π/12

π/2

7π/12

0

2π/3

5

10

15
Length of edges

20

25

Figure 4: Graphs of relative frequency vs. angle (left) and length (right) for the second sample domain and comparison with the graphs for the
FreeFem++ generated mesh.

In order to see the diﬀerence of finite element solutions by the two meshes of the example (i.e., the centre5

(FreeFem++) and rightmost panels (ours) of Figure 3), we solve a Poisson problem with an inhomogeneous Dirichlet boundary condition, where the external force and the boundary value of the problem are given so that the exact
solution is φ(x) ≡ sin(x1 /100) sin(x2 /100). Using the linear (P1) finite element in the FEM, we obtain two finite
element solutions φhFF and φhGS from the mesh generated by FreeFem++ and our mesh, respectively. Let k · kH01
be the norm of the Sobolev space H01 (Ω). The relative errors between the finite element and exact solutions are
kφhFF − φkH01 /kφkH01 = 5.37 × 10−2 and kφhGS − φkH01 /kφkH01 = 4.86 × 10−2 , which imply that φhGS is almost 10 percent
closer to φ than φhFF in the sense of H01 (Ω)-norm.
The computational time for our meshes is much longer than that for the meshes generated by FreeFem++. The
actual computational times for the second sample domain are almost 85 (ours, n = 55, 000) and 0.01 s (FreeFem++)on a notebook computer (Intel Core i7 1.8GHz, 4GB memory), i.e., the computational times diﬀer by a factor of about
8,500. This is the main disadvantage of our mesh generator, although the generated meshes are quite good. Step 3
accounts for the largest part of the computational time of our mesh generation (obtaining the solution of the GS modelby using the FDM). Since the FDM works well in general-purpose computing on graphics processing units (GPGPU)
and shows high performance, the disadvantage relating to the computational time might be considered tolerable.
In general, the computational time of our mesh generation is not so sensitive to the complexity of the domain
shape and mainly depends on the size of Ω0 . Indeed, for a star-shaped domain with the same Ω0 = (0, 250) × (0, 350)
it is almost 90 s (n = 60, 000), and for a twice larger (and similar) domain of Figure 3 (leftmost figure) with Ω0 =(0, 500) × (0, 700) it is almost 455 s (n = 70, 000).
4. Conclusions
We have proposed a self-organized mesh generator based on a spot pattern of the GS model. It shows the featuresof a self-replicating system, i.e., it fills the domain with spots, fits the spot pattern to the domain shape, and maintains
almost equal distances between neighbouring spots. We have discussed a (pixel) image-based implementation, whichworks well and is useful for practical digital images. The domain shape and the mesh size are controlled by d in (1),which gives spatially dependent diﬀusion coeﬃcients Du and Dv in the GS model. The mesh generator has a uniqueadvantage: the fairly simple procedure does not require additional complex processes and criteria, such as smoothing,(local) addition, and removal of nodes. These additional processes are automatically performed by the non-linear
system (GS), i.e., through the self-organization mechanism. Although only a two-dimensional mesh generator has
been realized here, the concept of the generator is expected to be useful for the development of a three-dimensionalmesh generator. There is a possibility that the generator can control the mesh orientation, if Du and Dv in (GS) are
modified as matrix valued diﬀusion coeﬃcients and a modified (Delaunay) triangulation technique is employed.
Acknowledgments
This study was supported by the Meiji University Global COE Program ‘Formation and Development of Mathematical Sciences
Based on Modeling and Analysis’ and the Meiji University Mathematical Modeling and Analysis Projects for Young Researchers(Development of a Mesh Generator Using a Self-Organizing Mechanism of a Reaction-Diﬀusion System). The first author was
supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan, through the Special Coordination
Funds for Promoting Science and Technology (the Tenure Track Program at Waseda Institute for Advanced Study), by the Japan
Society for the Promotion of Science (JSPS) through the Grant-in-Aid for Young Scientists (B) (No. 23740090) and the JapaneseGerman Graduate Externship (Mathematical Fluid Dynamics), and by Waseda University through the Grant for Special Research
Projects (No. 2011A-959).

References[1] A. M. Turing, The chemical basis of morphogenesis, Phil. Trans. R. Soc. London B 237 (1952) 37–72.
[2] H. Meinhardt, Models of Biological Pattern Formation, Academic Press, London, 1982.
[3] J. D. Murray, Mathematical Biology I, II, 3rd ed., Springer, New York, 2002.
[4] K. J. Lee, H. L. Swinney, Lamellar structures and self-replicating spots in a reaction-diﬀusion system, Phys. Rev. E 51 (1995) 1899–1915.
[5] K. J. Lee, W. D. McCormick, Q. Ouyang, H. L. Swinney, Pattern formation by interacting chemical fronts, Science 261 (1993) 192–194.

6

[6] P. De Kepper, V. Castets, E. Dulos, J. Boissonade, Turing-type chemical patterns in the chlorite-iodide-malonic acid reaction, Physica D 49
(1991) 161–169.
[7] P. De Kepper, J.-J. Perraud, B. Rudovics, E. Dulos, Experimental study of stationary Turing patterns and their interaction with traveling waves
in a chemical system, Int. J. Bifurcat. Chaos 4 (1994) 1215–1231.
[8] S. Kondo, R. Asai, A reaction-diﬀusion wave on the skin of the marine angelfish Pomacanthus, Nature 376 (1995) 765–768.
[9] A. Adamatzky, B. de Lacy Costello, T. Asai, Reaction-Diﬀusion Computers, Elsevier, Amsterdam, 2005.
[10] M. A. Yerry, M. S. Shephard, A modified quadtree approach to finite element mesh generation, IEEE Comput. Graph. 3 (1983) 39–46.
[11] P. J. Green, R. Sibson, Computing Dirichlet tessellations in the plane, Comput. J. 21 (1978) 168–173.
[12] S. W. Sloan, A fast algorithm for constructing Delaunay triangulations in the plane, Adv. Eng. Software 9 (1987) 34–55.
[13] T. Taniguchi, Automatic mesh generation for FEM - Application of Delaunay triangulation, Morikita, 1992, (in Japanese).
[14] R. Löhner, Progress in grid generation via the advancing front technique, Eng. Comput. 12 (1996) 186–210.
[15] K. Shimada, D. Gossard, Bubble mesh: automated triangular meshing of non-manifold geometry by sphere packing, Proc. third ACM symp.
(1995) 409–419.
[16] P. Gray, S. K. Scott, Autocatalytic reactions in the isothermal continuous stirred tank reactor, Chem. Eng. Sci. 38 (1983) 29–43.
[17] J. E. Pearson, Complex patterns in a simple system, Science 261 (1993) 189–192.
[18] Y. Nishiura, Far-from-Equilibrium Dynamics, AMS, 2002.
[19] Y. Nishiura, D. Ueyama, A skeleton structure of self-replicating dynamics, Physica D 130 (1999) 73–104.
[20] B. Delaunay, Sur la sphère vide, Bull. Acad. Sci. USSR 6 (1934) 793–800.
[21] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, SIAM, Philadelphia, 2002.
[22] FreeFem++, http://www.freefem.org/.

7