Deleting mesh elements from a mesh

Question

I am trying to make an interpolation function over my data. The points I am selecting make an unstructured grid, so I am following through the steps in the documentation "ref/message/Interpolation/udeg" from the message given after trying to interpolate over an unstructured grid. The following is not what I am looking at specifically, but I believe it captures the essence of my problem.

Let's say the points we want to define the interpolation function over is set up like this:

points = Flatten[Table[{x, z}, {x, 0, 10}, {z, 0, x^2, x^2/9}], 1];

so that we essentially have a region of points "cut off" by the function z=x^2

Following the steps in the documentation, I create a mesh:

Needs["NDSolve`FEM`"]
testMesh = ToElementMesh[points]

And to visualize the mesh with the boundary function I use the HighlightMesh function to get the following plot:

Now, I want this mesh to be of a higher order, so still following the documentation I increase the order:

testMesh2 = MeshOrderAlteration[testMesh, 2]

And then visualize again with the HighlightMesh function:

The process of making the new mesh puts midpoints between each of the original points.

The problem is that I don't want any points to the left of the boundary function.From the documentation it seems like I cannot just throw those coordinates out, since the entire second order mesh is needed for the interpolation to be of second order. When I try removing the triangle elements corresponding to these points and creating a new 1st order mesh, I get another error saying the set of mesh elements did not contain the 1st coordinate, so it seems like these extra elements are what are keeping that lower left coordinate in the mesh.

How can I create a mesh with my coordinates while excluding any mesh elements to going beyond that boundary function? I know I could just define the boundary lines and create a mesh in the region, but the initial points are picked for a specific reason, so I can't just initially have points chosen for me.

Of course, this problem could have some other solution that I am not considering here. If so, bringing that to my attention would be awesome.

Answer 1

Here's a way of constructing a mesh with the properties you want: start by constructing your initial mesh in the way you did:

points = Flatten[Table[{x, z}, {x, 0, 10}, {z, 0, x^2, x^2/9}], 1];
Needs["NDSolve`FEM`"]
testMesh = ToElementMesh[points]

As we can see in your plot, to get rid of the region you are not interested in we need to drop some 2-dimensional cells from the mesh. The cells we need to drop are those which contain the "first point" of the mesh so let's find them:

In[]:= labels1d = First@First@testMesh[[2]];
First /@ Position[labels1d, 1];
cellstodrop = labels1d[[#]] & /@ %

Out[]= {{1, 2, 3}, {1, 3, 4}, {4, 5, 1}, {1, 5, 6}, {6, 7, 1}, {7, 8, 1}, {9, 10, 1}, {10, 11, 1}, {8, 9, 1}, {21, 1, 11}, {31, 41, 1}, {41, 51, 1}, {1, 21, 31}, {1, 71, 81}, {91, 1, 81}, {1, 91, 101}, {71, 1, 61}, {61, 1, 51}}

These are the cells which we need to drop from the full list of 2-dim cells. So the new list of 2-dim cells becomes:

 cellstokeep = Complement[First@First@testMesh[[2]], cellstodrop]

We check that the mesh with the new 2-dimensional cells is the sought one:

 mesh1 = HighlightMesh[MeshRegion[testMesh[[1]], Polygon@cellstokeep, Frame -> True, AspectRatio -> 1/2], 0] 

We can now use this mesh to carry out your interpolation procedure:

 newmesh = ElementMesh@mesh1;
 testMesh2 = MeshOrderAlteration[newmesh, 2];
 HighlightMesh[MeshRegion[testMesh2, Frame -> True, AspectRatio -> 1/2], 0]

You mentioned that you attempted to remove the "triangles" (= 2-dim cells) of the region you don't need but it seems that somehow you also removed points. In this procedure we only removed 2-dim cells, but no point of the mesh got removed.