12
$\begingroup$

(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]]

x y values

Now I invent some data for each coordinate

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

values

Now I can make a mesh with my coordinates

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

mesh

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]

Contour plot

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]

bad mesh

bad contour plot

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]

Good mesh

good contour plot

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

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

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]

enter image description here

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

$\endgroup$
12
  • $\begingroup$ Well, you know, I thought I tried that.... (+1) $\endgroup$ – Michael E2 Dec 2 '20 at 6:13
  • 2
    $\begingroup$ @MichaelE2, maybe I need to document that better. $\endgroup$ – user21 Dec 2 '20 at 6:13
  • $\begingroup$ 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]. $\endgroup$ – Michael E2 Dec 2 '20 at 6:17
  • $\begingroup$ @MichaelE2 returns unevaluated for me, but if possible, I could add that. $\endgroup$ – user21 Dec 2 '20 at 6:19
  • 1
    $\begingroup$ @MichaelE2, good to know that it was not my code the screwed it up... $\endgroup$ – user21 Dec 2 '20 at 15:59
8
$\begingroup$

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] ]

enter image description here

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$ – Hugh Dec 1 '20 at 16:24
  • $\begingroup$ One has to force ToElementMesh to create a simple triangle mesh! I modified my answer. $\endgroup$ – Ulrich Neumann Dec 1 '20 at 18:18
  • 1
    $\begingroup$ 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. $\endgroup$ – user21 Dec 2 '20 at 6:08
7
$\begingroup$

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"]

enter image description here

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]

enter image description here

Check equivalence:

myifn[##] === ifn[##] & @@
 Transpose@RandomPoint[MeshRegion@emesh, 10000]
(*  True  *)
$\endgroup$
1
  • $\begingroup$ If one does not want to mess with the FEM context, your interpolation function approach is the way to go. $\endgroup$ – user21 Dec 2 '20 at 6:17
3
$\begingroup$

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]

Images

$\endgroup$
2
  • $\begingroup$ 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)? $\endgroup$ – Hugh Dec 2 '20 at 11:29
  • 1
    $\begingroup$ @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. $\endgroup$ – Tim Laska Dec 2 '20 at 14:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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