Algorithm [preferably fortran] to interpolate data from a 2D unstructured grid to form a cartesian grid

Question

I am modeling fracture propagation using a 2D dynamic unstructured grid. As the fracture propagate over time, the elements move accordingly. For a given time step, I would like to interpolate the data of my unstructured grid zi=f(xi,yi) (where (xi,yi) are the nodes of the unstructured triangles) to obtain the value of the function z on the Cartesian grid (where x and y regularly spaced). My grid size is about 100x100 nodes. There are many ways of doing this interpolation scheme but I do not know enough about the topic to decide which scheme is robust and simple enough. I am using Fortran 90 (do not ask why...). Is there any open-source algorithm available? I do not want to reinvent the wheel.

Thanks a lot!

Answer 1

The Earth System Modeling Framework (ESMF) may be of interest to you:

http://www.earthsystemmodeling.org/

It comes with Fortran, C and more recently Python (limited functionality) APIs and excellent documentation. It is commonly used for data abstraction objects in geophysical models and their coupling, but it also provides routines for regridding between structured and unstructured grids which can be used for offline (stand-alone command line utility) or online (via subroutine calls) interpolation weight generation.

Last time I have looked into it, ESMF provided bilinear, bicubic and quatitity conserving regridding methods.

Answer 2

Even though some references were great (ESMF for instance), I could not find a simple algorithm easily available. So it is sometimes easier to "re-invent" the wheel! In case some of you are interested, here is the simple method I derived:

1. Generate the background Cartesian grid: Node P is at l=j+Nbx(i-1) for i=1:Nbx and j=1:Nby.
2. Set up the node numbering for each background square. There are Nbx*Nby nodes and (Nbx-1)(Nby-1) squares. The squares are labeled as s=j+(i-1)(Nby-1) for i=1:Nbx-1 and j=1:Nby-1. Square S is made of 4 nodes: j+Nby(i-1), j+1+Nby(i-1), j+1+Nby*i, and j+Nby*i
3. For each triangle of the unstrusctured mesh, find the minimun and maximum x and y to describe a box including the triangle. Then, determine the corresponding imin,imax,jmin,jmax.
4. For the given box, go through the Cartesian nodes included in that box (i=imin:imax, j=jmin:jmax) and do the "Point_Inclusion_in_Triangle" Test (using barycentric method) to see if the node is inside the triangle of interest. When doing so, do an additional test to see if the node is on one of the vertices.
5. If the test is positive, the corresponding node of the background Cartesian mesh is eliminated from the list of potential nodes (a node cannot be part of multiple triangles). This makes the process faster we go through the list of triangles.
6. Then, you can do the interpolation at each node of the background mesh using the bilinear interpolaton method in a triangle (using barycentric method).

This method works fine and is pretty fast (for engineering pruposes).