13
$\begingroup$

I would like to know a good method for interpolating data between two unstructured grids, where one grid is a coarser version of the other.

Efficiency is very important to me since I'm solving a transient PDE problem where I need to transfer data between the grids at every time step of the solution.

I thought about using kd-tree for searching the nearest node of a given point, then I would use the shape functions of that element (FEM simulation) to interpolate the data. Is this a good solution? Are there better ones?

Do you also know any robust and reliable library in C/C++ for this task?

*I know there is a similar question, but it asks for the most accurate method on a structured grid.

Improve this question
$\endgroup$
1

4 Answers 4

Reset to default
7
$\begingroup$

Unstructured grids have their place.

You may want to look at the Earth System Modeling Framework (ESMF). They have some code for re-gridding -- specifically for this purpose -- and they've done some nifty stuff with parallel code, too. The whole system is designed to couple models, so there may be other useful stuff there as well.

Some other notes:

"no way to do this efficiently for any significant number of points"

well, efficient is a relative thing -- once you've got the grid in a tree structure, you can search it in O(logn), which can be pretty darn fast, though not O(1), as searching a regular grid is.

Also, it sounds like while the interpolation needs to be done at every time step, if the grids aren't adapting, then the mapping from one grid to another remains constant. So you can compute that mapping (i.e. which element in each grid corresponds to which element in the other) in whatever way is convenient, store it, and then you never need to compute it agin (until the grids change).

That leaves you with the interpolation code -- where you will want to balance accuracy with performance -- simple linear interpolation across a triangle is fast, and may be good enough.

"I thought about using kd-tree for searching the nearest node of a given point, then I would use the shape functions of that element"

remember that the nearest node doesn't get you the element -- so you'll want to do a bit more to find the element you want. One option would be to use an rtree instead, which stores/searches by bounding box -- you'll get more than one element with each search, but you can then check which of those is correct directly.

Improve this answer
$\endgroup$
2
  • $\begingroup$ This looks nice. I don't need to adapt the meshes, so the mapping from one grid to the other will be made only once. Thanks for the tip about the r-tree data structure. $\endgroup$
    – Bernardo M. R.
    Jun 28, 2012 at 18:38
  • 1
    $\begingroup$ I still stand by my view that this can't be done efficiently :-) There are actually two problems: (i) Typically you have $O(N)$ points, so an $O(\log N)$ lookup for each point still leads to a superlinear overall behavior. (ii) The most expensive part is actually the interpolation on the cell you have just found because you have to transform back into the coordinate system of that cell, evaluate shape functions, multiply with coefficients, etc. That's a lot of work if you need to do it for large numbers of interpolation points. $\endgroup$
    – Wolfgang Bangerth
    Jun 28, 2012 at 22:25
6
$\begingroup$

If I understand you correctly, you want to fill in the values of the finer grid by interpolating on the coarser grid. One way to do linear interpolation on an unstructured grid is with Delaunay triangulations (this is how Matlab's griddata and TriScatteredInterp commands are implemented). After constructing a triangulation of your grid points, interpolation boils down to locating the triangle containing the target point, calculating its barycentric coordinates, and using the function values at the vertices to compute the interpolated value. CGAL can construct n-dimensional triangulations (for medium n), and also has a built-in 2d interpolation module.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Yes. But I also want to "inject" the values from the fine grid into the coarse grid as well, that's why I said transfer. $\endgroup$
    – Bernardo M. R.
    Jun 28, 2012 at 11:26
3
$\begingroup$

This is what I'm doing at the moment, except I'm transferring function values at quadrature points, not nodes. I'm implementing the technique explained in the chosen answer to my question here: Finding which triangles points are in .

Basically, say you have 2 grids $A$ and $B$, and want to transfer information from grid $A$ to grid $B$

  1. Based on node or quadrature points in grid $B$, construct a list of points $p_i$ to be evaluated in grid $A$,
  2. order the evaluation points $p_i$ based on a Hilbert curve,
  3. construct the constrained Delaunay extension of grid $A$ (all triangles in $A$, plus some more exterior triangles to make the shape convex)
  4. starting at a triangle in the extended version of $A$, randomly walk from triangle to triangle "in the general direction" of point $p_1$ until you get there. Then keep walking to $p2$, $p3$, etc. until you have found which triangles in $A$ contain all the evaluation points.

If you have $N$ evaluation points and $M$ triangles in the extended version of $A$, and if you ignore the Hilbert curve ordering, and you walk deterministically, worst case the time could be $O(N \sqrt{M})$. Including the Hilbert curve ordering and walking with randomness, I'm seeing much more efficient results in practice, more like $O(max(N,M))$.

Improve this answer
$\endgroup$
2
$\begingroup$

This is the kind of job for which you really want to avoid unstructured meshes since there is no way to do this efficiently for any significant number of points. You should consider using meshes that are at least somehow related to each. For example, if they are both obtained from hierarchical refinement of a coarse mesh, then you can relatively easily and efficiently find out where the interpolation points of one mesh are located on the other mesh.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I think that might be the best option (hierarchy of grids). If this is the case, do you know any good data structure, or specific method to use? $\endgroup$
    – Bernardo M. R.
    Jun 28, 2012 at 11:29
  • $\begingroup$ Yes, hierarchic meshes are all stored as quad/oct-trees (if they start on a single coarse cell) or forests of such trees (if the coarse mesh has more than once cell). $\endgroup$
    – Wolfgang Bangerth
    Jun 28, 2012 at 22:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.