Mathematics and Computers in Simulation 58 (2002) 125–132
A simple efficient algorithm for interpolation between different grids in both 2D and 3D Jichun Li∗ , C.S. Chen Department of Mathematical Sciences, University of Nevada, 4505 Maryland Parkway, Box 454020, Las Vegas, NV 89154-4020, USA Received 1 February 2001; accepted 14 March 2001
Abstract Considering the importance of data transferring between different grids, we present a simple but powerful interpolation scheme using radial basis functions (RBFs) to accomplish such task in both 2D and 3D. Numerical results are presented for such interpolation between truly unstructured grids. © 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Interpolation; Radial basis function and unstructured grid
1. Introduction Due to the continuing innovation of modern computers, numerical simulation has become a common tool in different scientific areas. With so many legacy modeling packages developed in the past several decades, there is a trend to couple different models (originally developed for some specific physics or small individual systems) to model large scale areas or global systems because many physics are naturally interconnected. For instance, the surface runoff carries organic and inorganic contaminants into aquifer recharge zones and rivers. Hence, surface modeling, subsurface modeling and ground-water modeling should be linked together for better understanding of the whole environmental system. Because of the independent development of individual model and its underlying different physics, it is very common that different models have very different grids. For example, coupling ADCIRC [13] (an advanced circulation model for shelves, coasts and estuaries based on unstructured triangular grids) with CE-QUAL-ICM [2] (a widely used multidimensional water quality model by this model employees an unstructured grid system, e.g. unstructured quadrilateral grids) requires transferring the velocity field obtained from ADCIRC on unstructured triangular grids to a velocity field used by CE-QUAL-ICM on quadrilateral grids or other different grids. When coupling component system models such as the ∗
Corresponding author. Tel.: +1-702-895-0357; fax: +1-702-895-4343. E-mail address:
[email protected] (J. Li). 0378-4754/02/$ – see front matter © 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 4 7 5 4 ( 0 1 ) 0 0 3 4 8 - 2
126
J. Li, C.S. Chen / Mathematics and Computers in Simulation 58 (2002) 125–132
Fig. 1. An examplary overlapping mesh for the classical alternating Schwarz method.
NCAR CCM3 and LANL POP ocean model, the situation is that CCM3 uses a low-resolution Gaussian latitude–longitude grid, while POP uses a high-resolution displaced pole grid [17]. Once again, we have to do different grid interpolation for coupling these two models. Hence, accurate data transferring between different grids becomes a critical issue to model coupling. Such data interpolation between different grids also occurs in overlapping domain decomposition methods, where the overlapped area is usually covered by different grids (an exemplary picture is shown in Fig. 1), the implementation of the classical alternating Schwarz method involves such data interpolation between different grids [19]. Similar situation also happens in the mortar element method or multiblock method [1,20], where we have to transform information between the subdomain interfaces when nonmatching grids are used in different blocks. The traditional method for interpolating data between different grids in the finite element method is by using linear or bilinear interpolation scheme [19], which accuracy is low and the implementation becomes very tedious for truly unstructured grids. In this paper, we propose a simple but very efficient algorithm using the state-of-the-art radial basis functions (RBFs) to alleviate the difficulty in coupling different unstructured grid in both 2D and 3D. Over the past 3 decades, RBFs have emerged as an important area in the approximation theory in the mathematical literature. It has been proven to be very effective in interpolating multivariate scattered data in many science and engineering problems. The idea of using RBFs for interpolating scattered multivariate data was first developed by Hardy in 1971 in modeling the earth’s gravitational field [10]. It is a grid free scheme for representing surfaces in an arbitrary number of dimensions. Many successful applications of RBFs were reported in Hardy’s review paper [11] such as geoscience, photogrammetry, topography, estimation, signal processing, hydrology and artificial intelligence, etc. In recent years, RBFs have also been extended to solving partial differential equations. Currently, perhaps the most active area of RBF application is the field of neural network modeling. The significant of Hardy’s work was apparently not receiving much attention in the mathematics community until Franke’s paper in 1982 [5] and Micchelli’s paper in 1986 [15]. Frank’s review paper compared virtually all of the interpolation methods for scattered data sets available at that time. Based on the following criteria: accuracy, visual aspect, sensitivity to parameters, execution time, storage requirements and ease of implementation, Hardy’s multiquadrics (MQ) were ranked the best, followed by Duchon’s thin plate splines (TPS) [3] and several other RBR interpolants. Micchelli gave a proof on the solvability of RBFs and hence established the rigorous theoretical foundation of RBFs. During the past decade, the use of RBFs in multivariate approximation has become a rapidly growing field in
J. Li, C.S. Chen / Mathematics and Computers in Simulation 58 (2002) 125–132
127
mathematical research. We refer reader to a good review paper of the theory of RBFs interpolation by Powell [18]. The paper is organized as follows. In Section 2 we give a detailed description about the algorithm. Extensive numerical results carried out for both 2D and 3D unstructured grids are presented in Section 3. We close with some discussions for further research in Section 4. 2. The algorithm In our interpolation scheme, we use the RBF, which in Rn (n ≥ 2) is defined as a function of the form φ(P − Q), where · denotes the Euclidean distance between the points P and Q, and φ : R+ → R be a continuous function with φ(0) ≥ 0. There are many types of radial basis functions, the widely used ones are TPS, MQ, the Gaussians and compactly supported RBFs. More details about RBFs can be found in [6–8]. For simplicity we only consider MQ, where φ(r) = (r 2 + c2 )β/2 ,
β is an odd integer, constant c > 0.
(1)
The free parameter c is often referred to as the shape parameter, whose choice can greatly affect the accuracy of the approximation. How to choose this shape parameter is still an outstanding research problem. However, there are various empirical results which have shown success in determining this shape parameter. Hardy [10] suggested choosing c to be 0.815 times the average distance from each datum point to its nearest datum point. Golberg et al. [9] applied the method of cross-validation, which is a statistical method, for choosing the optimal shape parameter with excellent results. For details of the cross-validation and other related methods for choosing c, we refer reader to [9] and references cited therein. Other than MQ, there are many RBFs such as r 3 and higher order polyharmonic splines [16] which do no have the free parameter and can also be effective. To consider the interpolation between two different grids, we assume that we are given two different meshes in Rn , n ≥ 2, say, m1 n • Mesh 1 is composed of m1 points {Q[1] j }1 in R , which are connected by the connectivity matrix M1 , [1] m1 m1 and the corresponding values {f (Q[1] j )}1 are given at {Qj }1 ; m2 n • Mesh 2 is composed of m2 points {Q[2] j }1 in R , which are connected by the connectivity matrix M2 .
We want to find highly accurate interpolated values at those points of Mesh 2. Let f (P ) =
m1
aj (rj2 + c2 )β/2 + a0 ,
j =1 1 where the unknown coefficients {aj }m j =1 ∪ {a0 } can be obtained by Collocation method; i.e.
m1
[1] 2 2 β/2 aj (Q[1] + a0 = f (Q[1] i − Qj + c ) i ),
1 ≤ i ≤ m1 ,
(2)
j =1 m1 j =1
aj = 0.
(3)
128
J. Li, C.S. Chen / Mathematics and Computers in Simulation 58 (2002) 125–132
Fig. 2. (a) The given velocity field on an unstructured quadrilateral mesh with 149 vertices and 114 elements; (b) the interpolated velocity field from Fig. 2(a) to an unstructured quadrilateral mesh with 385 vertices and 315 elements.
Micchelli [15] proved that MQ is conditional negative definite of order 0 and hence the coefficient matrix resulting from Eqs. (2) and (3) is negative definite which means they are solvable. In many practical applications, it is observed that the additional term a0 can be set to zero without affecting the accuracy of m2 the method. Once we solve the system of Eqs. (2) and (3), we can evaluate f (P ) at those points {Q[2] j }1 of Mesh 2. For certain analytic functions, Madych [14] showed that the MQ interpolant converges exponentially to f . Note that our interpolation scheme does not use any information about the connectivity, hence it can be used for any unstructured mesh given only the nodal coordinates. Furthermore, it is not difficult to see that the implementation of the RBF interpolation scheme is all the same for any dimensional meshes, the only difference is in calculating the distance r. Note that similar scheme was used for interpolation and approximation of 2D, 3D and 4D scattered data [4,12].
Fig. 3. (a) The interpolated velocity field from Fig. 2(a) to an unstructured quadrilateral mesh with 923 vertices and 843 elements; (b) the interpolated velocity field from Fig. 1(a) to an unstructured triangular mesh with 112 vertices and 142 elements.
J. Li, C.S. Chen / Mathematics and Computers in Simulation 58 (2002) 125–132
129
Fig. 4. (a) The interpolated velocity field from Fig. 2(a) to an unstructured triangular mesh with 310 vertices and 468 elements; (b) the interpolated velocity field from Fig. 1(a) to an unstructured triangular mesh with 868 vertices and 1481 elements.
3. Numerical results Considering the direct application for our project UTPROJ [21], we tested the RBF scheme for interpolating velocity field between different unstructured grids in both 2D and 3D. To investigate the accuracy of the algorithm, we assume that the exact velocity field in 2D is given as u = −π cos(πx) cos(0.5πy),
v = 0.5π sin(πx) sin(0.5πy).
(4)
While in 3D, the velocity field is given as u = −π cos(πx) cos(0.5πy)(1 + z)2 ,
(5)
v = 0.5π sin(πx) sin(0.5πy)(1 + z)2 ,
(6)
w = 2 sin(πx) cos(0.5πy)(1 + z).
(7)
Fig. 5. (a) The given velocity field on an unstructured triangular mesh with 310 vertices and 468 elements; (b) the interpolated velocity field from Fig. 5(a) to an unstructured triangular mesh with 112 vertices and 142 elements.
130
J. Li, C.S. Chen / Mathematics and Computers in Simulation 58 (2002) 125–132
Fig. 6. (a) The interpolated velocity field from Fig. 5(a) to an unstructured triangular mesh with 868 vertices and 1481 elements; (b) the interpolated velocity field from Fig. 5(a) to an unstructured quadrilateral mesh with 149 vertices and 114 elements.
In our implementation, we chose β = 1, i.e. the classical Hardy MQ [10], which were proved to have the best approximation properties [14]. Also we fixed the parameter c = 2 in the definition of MQ (1) for all our tests. The system of Eqs. (2) and (3) is solved by the standard Gaussian elimination method. To show our algorithm’s generality, we presented the interpolated velocity field from an unstructured quadrilateral mesh to different unstructured triangular and quadrilateral meshes in Figs. 2–4. The interpolated velocity field from an unstructured triangular mesh to unstructured triangular and quadrilateral meshes are given in Figs. 5–7. A 3D example is provided in Figs. 8 and 9 for transferring the given velocity field on a cubic mesh to a tetrahedron mesh and another cubic mesh. The results are very impressive. The interpolated velocity field is accurate up to 5 significant digits compared to the exact velocity field. The computer time for each interpolation is all less than 1 s on our Pentium III PC with 700 MHz CPU and 256 K Cache.
Fig. 7. (a) The interpolated velocity field from Fig. 5(a) to an unstructured quadrilateral mesh with 385 vertices and 315 elements; (b) the interpolated velocity field from Fig. 5(a) to an unstructured quadrilateral mesh with 923 vertices and 843 elements.
J. Li, C.S. Chen / Mathematics and Computers in Simulation 58 (2002) 125–132
131
Fig. 8. (a) The given velocity field on a cubic mesh with 343 vertices and 216 elements.
Fig. 9. (a) The interpolated velocity field from Fig. 8(a) to a tetrahedron mesh with 75 vertices and 218 elements; (b) the interpolated velocity field from Fig. 8(a) to acubic mesh with 729 vertices and 512 elements.
4. Conclusions To alleviate the difficulty of coupling different models, a simple and efficient algorithm using radial basis functions for data transferring in both 2D and 3D is presented in this paper. Numerical results are given to demonstrate the accuracy and efficiency of the algorithm. It is not difficult to see that similar interpolations can be carried out for other radial basis functions [6–8].
Acknowledgements The first author would like to thank Professor Mary Wheeler for bringing up the issue about the data transferring between different grids to him.
132
J. Li, C.S. Chen / Mathematics and Computers in Simulation 58 (2002) 125–132
References [1] T. Arbogast, L.C. Cowsar, M.F. Wheeler, I. Yotov, Mixed finite element methods on nonmatching multiblock grids, SIAM J. Numer. Anal. 37 (2000) 1295–1315. [2] C.F. Cerco, T. Cole, User’s Guide to the CE-QUAL-ICM Three-Dimensional Eutrophication Model, Release Version 1.0, Technical Report EL-95-15, US Army Engineer Waterways Experiment Station, Vicksburg, MS, 1995. [3] J. Duchon, Splines minimizing rotation invariant seminorms in Sobolev spaces, in: W. Schemmp, K. Zeller (Eds.), Constructive Theory of Functions of Several Variables, Lecture Notes in Mathematics, Springer, Berlin, 1976, pp. 85–100. [4] T.A. Foley, Interpolation and approximation of 3D and 4D scattered data, Comput. Math. Appl. 13 (1987) 711–740 . [5] R. Franke, Scattered data interpolation: tests of some methods, Math. Comput. 48 (1982) 181–200. [6] M.A. Golberg, C.S. Chen, Discrete Projection Methods for Integral Equations, Computational Mechanics Publications, Southampton, 1997. [7] M.A. Golberg, C.S. Chen, The method of fundamental solutions for potential, Helmholtz and diffusion problems, in: M.A. Golberg (Ed.), Boundary Integral Methods — Numerical and Mathematical Aspects, Computational Mechanics Publications, Southampton, 1998, pp. 103–176. [8] M.A. Golberg, C.S. Chen, A bibliography on radial basis function approximation, Boundary Element Communications 7 (1996) 155–163. [9] M.A. Golberg, C.S. Chen, S.R. Karur, Improved multiquadric approximation for partial differential equations, Eng. Anal. Boundary Elements 18 (1996) 9–17. [10] R.L.Hardy, Multiquadric equations of topography and other irregular surfaces, J. Geophys. Res. 176 (1971) 1905–1915. [11] R.L.Hardy, Theory and applications of the multiquadric-biharmonic method: 20 years of discovery 1968–1988, Comput. Math. Appl. 19 (1990) 163–208. [12] E.J. Kansa, Multiquadrics — a scattered data approximation scheme with applications to computational fluid dynamics I, Comput. Math. Appl. 19 (2000) 127–145. [13] R.A. Luettich, J.J. Westerink, N.W. Scheffner, ADCIRC: An Advanced Three-Dimensional Circulation Model for Shelves, Coasts and Estuaries, Report 1, US Army Corps of Engineers, Washington, DC, 20314-1000, December 1991. [14] W.R. Madych, Miscellaneous error bounds for multiquadric and related interpolants, Comput. Math. Appl. 24 (1992) 121–138. [15] C.A. Micchelli, Interpolation of scattering data: distance matrices and conditionally positive definite functions, Constr. Approx. 2 (1986) 11–22. [16] A.S. Muleshkov, M.A. Golberg, C.S. Chen, Particular solutions of Helmholtz-type operators using higher order polyharmonic splines, Comput. Mech. 23 (1999) 411–419. [17] NCAR 1996 Annual Scientific Report: Computational Science Section, available at http://www.ncar.ucar.edu/archives/ asr/ASR96/scd/css.html. [18] M.J.D. Powell, The theory of radial basis function approximation in 1990, in: W. Light (Ed.), Advances in Numerical Analysis, Vol. II, Oxford Science Publications, Oxford, 1992, pp. 105–210. [19] B. Smith, P. Bjorstad, W. Gropp, Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press, Cambridge, 1996. [20] M.F. Wheeler et al., A parallel multiblock/multidomain approach for reservoir simulation, SPE 51884, in: Proceedings of the 15th SPE Reservoir Simulation Symposium, Houston, TX, 14–17 February 1999. [21] M.F. Wheeler, C. Dawson, J. Li, V. Parr, UTPROJ: The University of Texas Projection Code for Computing Locally Conservative Velocity Fields, ERDC, US Army Engineer Research and Development Center, MSRC (Major Shared Resource Center), Technical Reports 99-04, Vicksburg, MS, March 1999.