4
$\begingroup$

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 doesn`t change but it takes time to generate Denaunay triangulation at each step.

Improve this question
$\endgroup$
5
  • $\begingroup$ Look up NDSolve`FEM`ElementMeshInterpolation. $\endgroup$
    – Michael E2
    Commented Sep 25, 2015 at 9:17
  • 1
    $\begingroup$ 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? $\endgroup$
    – Michael E2
    Commented Sep 25, 2015 at 9:23
  • $\begingroup$ 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]]}]. $\endgroup$
    – Artemiy
    Commented Sep 28, 2015 at 8:18
  • $\begingroup$ Related: mathematica.stackexchange.com/questions/68973/… (but the grid has a semi-regular structure). $\endgroup$
    – Michael E2
    Commented Sep 28, 2015 at 10:19
  • $\begingroup$ Thanks, Michael E2. NDSolveFEMElementMeshInterpolation helps. Becomes x4 faster. $\endgroup$
    – Artemiy
    Commented Sep 29, 2015 at 15:15

1 Answer 1

Reset to default
5
$\begingroup$

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["NDSolve`FEM`"]

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.

Improve this answer
$\endgroup$

Your Answer

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

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