Interpolation on an unstructured mesh

Question

(Edit: This question is about making a mesh on a 2D surface to work with a 3D surface look at ElementMeshInterpolation on a BoundaryMesh )

I have experimental data in the form of { {x1, y1, z1}, {x2, y2, z2}...} I wish to interpolate it so my thoughts turned to using the finite element function ElementMeshInterpolation.

Here is a minimum working example of how some data might look. First I create my x and y values.

pts = Table[t { Cos[t], Sin[t]}, {t, 0, 100, 0.1}];
    Graphics[Point[pts]]

Now I invent some data for each coordinate

values = {#[[1]], #[[2]], #[[1]] #[[2]]^2} & /@ pts;
Graphics3D[Point[values], BoxRatios -> {1, 1, 0.5}]

Now I can make a mesh with my coordinates

Needs["NDSolve`FEM`"];
dm = DelaunayMesh[pts];
mesh = ToElementMesh[dm, "MeshOrder" -> 1];
mesh["Wireframe"]

Now the problem comes because when I do

int = ElementMeshInterpolation[{mesh}, values[[All, 3]]];

I get an error because the Delaunay mesh has generate more points than in my data.

{Length[pts], Length[mesh["Coordinates"]]}

(*  {1001, 1744}  *)

This is obvious on reflection. How do I create a mesh with only the points in my data? I can't use ToElementMesh because I don't know the triangulation of the points. Any suggestions? Thanks.

EDIT

Users Ulrich Neumann and Michael E2 have both come up with ways to keep the number of points in the mesh the same as the number of values. This is done by setting MeshQualityGoal -> 0. However, user21, who knows the code, thinks this is dubious and works only by luck. user21 has come up with a very simple solution which I illustrate here because I also want to demonstrate a further problem that I hit when I tried with my actual data.

Needs["NDSolve`FEM`"];
pts = Table[t { Cos[t], Sin[t]}, {t, 0, 100, 0.1}];
values = {#[[1]], #[[2]], #[[1]] #[[2]]^2} & /@ pts;
mesh = ToElementMesh[pts];
int = ElementMeshInterpolation[{mesh}, values[[All, 3]]];
ContourPlot[int[x, y], {x, y} ∈ mesh, 
 ContourShading -> False, Contours -> 25]

Now let me show a problem that happens when the x and y values are very different in scale. I rescale the data so that the x values are much smaller than before.

pts2 = pts /. {x_, y_} -> {x/1000., y};
mesh2 = ToElementMesh[pts2];
Show[mesh2["Wireframe"], AspectRatio -> 1]
int2 = ElementMeshInterpolation[{mesh2}, values[[All, 3]]];
ContourPlot[int2[x, y], {x, y} ∈ mesh2, 
 ContourShading -> False, Contours -> 25, PlotRange -> All]

As you can see the mesh is badly formed and this puts errors into the contour plot. The solution is to rescale the data as follows.

{x1, x2} = MinMax[pts2[[All, 1]]];
{y1, y2} = MinMax[pts2[[All, 2]]];
pts3 = {Rescale[#[[1]], {x1, x2}, {-1, 1}], 
     Rescale[#[[2]], {y1, y2}, {-1, 1}]} & /@ pts2;
mesh3 = ToElementMesh[pts3];
Show[mesh3["Wireframe"], AspectRatio -> 1]
int3 = ElementMeshInterpolation[{mesh3}, values[[All, 3]]];
ContourPlot[int3[x, y], {x, y} \[Element] mesh3, 
 ContourShading -> False, Contours -> 25, PlotRange -> All]

Everything is fine again except that all the data is scaled to {-1,1} in both directions. Thanks for the help.

Answer 1

Here is how to do it, just let ToElementMesh create the mesh:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[pts];
values = {#[[1]], #[[2]], #[[1]] #[[2]]^2} & /@ pts;
int = ElementMeshInterpolation[{mesh}, values[[All, 3]]];
{Length[pts], Length[mesh["Coordinates"]]}

{1001, 1001}

Plot3D[int[x, y], {x, y} \[Element] mesh]

Generallly speaking, given a set of points ToBoundayMesh will return a convex hull and ToElementMesh will return a Delaunay triangulation.

Answer 2

Try

mesh = ToElementMesh[DelaunayMesh@pts, MeshQualityGoal -> 0 ,"MeshOrder" -> 1 ];

This is a simple triangle mesh (no additional points!)

Length[pts]==Length[mesh["Coordinates"]]
(*True*) 

values = #[[1]] #[[2]]^2 & /@ mesh["Coordinates"];
fFE = ElementMeshInterpolation[{mesh}, values];
Plot3D[fFE[x, y], Element[{x, y}, mesh] ]

Answer 3

I think Interpolation does what you want under the hood (basically what @MarcoB said):

ifn = Interpolation[Transpose@{pts, values[[All, 3]]}, 
   InterpolationOrder -> 1];

emesh = ifn@"ElementMesh";
emesh["Wireframe"]

Note that it controlled the construction of the mesh in the way you wanted.

Update. You can construct the mesh this way:

mymesh = ToElementMesh[ConvexHullMesh@pts, MeshQualityGoal -> 0, 
   MaxCellMeasure -> Infinity, "IncludePoints" -> pts, 
   "MeshOrder" -> 1];
Normal@mymesh["Wireframe"] === Normal@emesh["Wireframe"]
(* True  *)

The coordinates are ordered differently from pts, so we need to permute the values accordingly to construct the ElementMeshInterpolation:

myvalues = 
  values[[All, 3]][[
    Ordering@ pts]][[
    InversePermutation@ Ordering@ mymesh@"Coordinates"]];
myifn = ElementMeshInterpolation[{mymesh}, myvalues];

Plot3D[myifn[x, y], {x, y} ∈ mymesh]

Check equivalence:

myifn[##] === ifn[##] & @@
 Transpose@RandomPoint[MeshRegion@emesh, 10000]
(*  True  *)

Answer 4

This answer refers to the latest edit in the OP about dramatically changing the scale of one dimension. The current mesher is geared towards making isotropic triangles. Therefore, it prefers a square domain.

Perhaps a simpler way to tackle the problem is to scale the coordinates of your high aspect ratio domain to be square and construct the mesh. Then, explicitly re-mesh by scaling the coordinates back to their high aspect ratio domain and keeping the connectivity the same.

Here's an example:

pts2 = pts /. {x_, y_} -> {x/1000., y};
(*Scale the coordinates so that the domain is square*)
pts3 = pts2 /. {x_, y_} :> {1000 x, y};
mesh2 = ToElementMesh[pts3];
Show[mesh2["Wireframe"]]
(*Rescale mesh coordinates back to original scale*)
pts4 = mesh2["Coordinates"] /. {x_, y_} :> { x/1000, y};
(*Re-mesh using rescaled coordinates*)
mesh3 = ToElementMesh["Coordinates" -> pts4, 
   "MeshElements" -> mesh2["MeshElements"]];
Show[mesh3["Wireframe"]]
int3 = ElementMeshInterpolation[{mesh3}, values[[All, 3]]];
ContourPlot[int3[x, y], {x, y} \[Element] mesh3, 
 ContourShading -> False, Contours -> 25, PlotRange -> All]