# Interpolation on an unstructured mesh

(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 = {#[], #[], #[] #[]^2} & /@ pts;
Graphics3D[Point[values], BoxRatios -> {1, 1, 0.5}] Now I can make a mesh with my coordinates

Needs["NDSolveFEM"];
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["NDSolveFEM"];
pts = Table[t { Cos[t], Sin[t]}, {t, 0, 100, 0.1}];
values = {#[], #[], #[] #[]^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[#[], {x1, x2}, {-1, 1}],
Rescale[#[], {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.

• Is Interpolation[{{#1, #2}, #3}& @@@ values] not an option because of the reduced interpolation order? – MarcoB Dec 1 '20 at 15:02
• Automatic rescaling is probably difficult, I'll keep it in mind though. – user21 Dec 2 '20 at 14:18
• The scaling problem is the underlying issue here: mathematica.stackexchange.com/questions/68973/… – Michael E2 Dec 2 '20 at 14:18

## 4 Answers

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

Needs["NDSolveFEM"]
mesh = ToElementMesh[pts];
values = {#[], #[], #[] #[]^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.

• Well, you know, I thought I tried that.... (+1) – Michael E2 Dec 2 '20 at 6:13
• @MichaelE2, maybe I need to document that better. – user21 Dec 2 '20 at 6:13
• I think maybe what I did was ToElementMesh["Coordinates" -> pts]. My computer crashed and I lost the code. I must not have thought of simply ToElementMesh[pts]. – Michael E2 Dec 2 '20 at 6:17
• @MichaelE2 returns unevaluated for me, but if possible, I could add that. – user21 Dec 2 '20 at 6:19
• @MichaelE2, good to know that it was not my code the screwed it up... – user21 Dec 2 '20 at 15:59

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 = #[] #[]^2 & /@ mesh["Coordinates"];
fFE = ElementMeshInterpolation[{mesh}, values];
Plot3D[fFE[x, y], Element[{x, y}, mesh] ] • Thanks but this is no good. I only have measured data at my coordinates. I don't have data at the mesh coordinates. The mesh has more coordinates than I have because Delaunay has filled in an absent area. – Hugh Dec 1 '20 at 16:24
• One has to force ToElementMesh to create a simple triangle mesh! I modified my answer. – Ulrich Neumann Dec 1 '20 at 18:18
• I don't think setting MeshQualityGoal->0 is a good idea. This has the potential to return a bad solution. On other words I think it's probably a coincidence that it works at all. – user21 Dec 2 '20 at 6:08

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  *)

• If one does not want to mess with the FEM context, your interpolation function approach is the way to go. – user21 Dec 2 '20 at 6:17

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] • Thank you. What a good approach. Would you call the fact that Plot3D and related functions don't do this automatically when given data that is unstructured a bug (or almost a bug)? – Hugh Dec 2 '20 at 11:29
• @HughYou are welcome! I think it could be difficult to implement automatic scaling in a general way. Depending on the model geometry and the physics you're trying to solve, I could see automatic scaling leading to unintended results. If you have a model with multiple material domains that require different mesh refinements, you may want to apply a different and potentially non-uniform scaling to each material domain for visualization. I may be exaggerating the complexity of the problem. – Tim Laska Dec 2 '20 at 14:12