# Repeated interpolation for 3D unstructured grid -> Optimization

Trying to perform mesh-to-mesh 3D interpolation on unstructured grid. Mesh size is about 200k nodes. Painful thing that interpolation shoud be repeated ~1000 times (different thermal map) on the same mesh.

One "operation" performs in less than a minute (w/ InterpolationOrder->1) --> but 1000 steps takes 12 hours. I am wondering are there any options to optimize the procedure for unstructured grid and make the whole calculation faster, cause mesh doesnt change but it takes time to generate Denaunay triangulation at each step.

• Look up NDSolveFEMElementMeshInterpolation. – Michael E2 Sep 25 '15 at 9:17
• If that's no help, then perhaps the question is not clear enough. How is your mesh generated? Can you provide a minimal example? What function are you trying to optimize? Might the optimal values be at the nodes (since you using linear interpolation)? Or did you mean performance-tuning instead of mathematical-optimization? – Michael E2 Sep 25 '15 at 9:23
• There are two list of nodes (for example): inputlist=RandomReal[{-10, 10}, {200000, 3}]  and outputlist=RandomReal[{-10, 10}, {200000, 3}]  (in these example boundaries may be different, but in my case they are equal). Both lists represent unstructured grid. Another list contains data in the inputlist at different timesteps (example): Temp=RandomInteger[{100,200},{200000,1000}]. I made a loop with TempFunc = Interpolation[DataTable, InterpolationOrder -> 1] function, where DataTable=MapThread[List, {inputlist, Temp[[istep]]}]. – Artemiy Sep 28 '15 at 8:18
• Related: mathematica.stackexchange.com/questions/68973/… (but the grid has a semi-regular structure). – Michael E2 Sep 28 '15 at 10:19
• Thanks, Michael E2. NDSolveFEMElementMeshInterpolation helps. Becomes x4 faster. – Artemiy Sep 29 '15 at 15:15

One can use generate an ElementMesh once and for all (15-16 sec.); then use ElementMeshInterpolation on each coordinate to construct interpolations (1.7 sec.).

Needs["NDSolveFEM"]

SeedRandom[0];
inputlist = RandomReal[{-10, 10}, {200000, 3}];
outputlist = RandomReal[{-10, 10}, {200000, 3}];

( mesh = DelaunayMesh[inputlist];
elem = Thread[MeshCells[mesh, 3], Tetrahedron] /.
Tetrahedron -> TetrahedronElement;
(*coords = MeshCoordinates[mesh];
coords == inputlist*)
emesh = ToElementMesh[
"Coordinates" -> inputlist,
"MeshElements" -> {elem}];
) // AbsoluteTiming

ifns = ElementMeshInterpolation[{emesh}, #,
"ExtrapolationHandler" -> {Indeterminate &,
"WarningMessage" -> False}] & /@
Transpose[outputlist]; // AbsoluteTiming
(*
{15.7423, Null}
{1.69745, Null}
*)


One can construct various interfaces to the interpolating functions. The matrix one is faster on lists of points, though only by a little.

Clear[ifn];
ifn[v_?VectorQ] := Through[ifns @@ v];
ifn[m_?MatrixQ] := Through[ifns @@ Transpose[m]];
ifn[x_?NumericQ, y_?NumericQ, z_?NumericQ] := Through[ifns[x, y, z]];


Evaluating the interpolating function is not fast, but I assume that is probably true no matter what. (An unstructured grid is handle by an ElementMesh under the hood, AFAIK.) In any case, constructing emesh only once saves a lot of time. I also am getting a small error with my last iteration. (There wasn't any the first time I tried.)

ifn@inputlist[[;; 10]] - outputlist[[;; 10]] // Abs // Max // AbsoluteTiming
(*  {0.584745, 1.06581*10^-14}  *)


But there doesn't seem to be a problem with the data in the interpolating functions:

Table[
{ifns[[i]]["ElementMesh"]["Coordinates"] == inputlist,
ifns[[i]]["ValuesOnGrid"] == outputlist[[All, i]]}, {i, 3}]
(*  {{True, True}, {True, True}, {True, True}}  *)
`

I suppose it's a weakness in the internal interpolating code.