Deleting mesh elements from a mesh

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["NDSolveFEM"]
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.

• Do you need to create such a mesh with your coordinates or are you looking for a way to create an interpolation over this region? – C. E. Jul 20 '17 at 18:16
• I guess I am really looking for an interpolation over this region. When I use my predetermined coordinates, the error message that follows supplies these steps of using the meshes, so I figured that was the best way to do it. If there is something better or easier to implement I would totally do that though. – Aaron Stevens Jul 21 '17 at 4:46
• ok, that's good to know in case anyone has another idea. I fear that it will be difficult to discretize that region using those points, but maybe someone knows how to do it. – C. E. Jul 21 '17 at 16:19
• We will see. Worst case scenario I can just define the mesh using a region with the curve as a boundary. It's just with the specifics of the actual problem I'm trying to solve this will take longer: both for the mesh generation and then solving for my function values at the mesh coordinates to actually perform the interpolation. – Aaron Stevens Jul 22 '17 at 0:58

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["NDSolveFEM"]
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.

• Thanks Alfonso! This seems to do the trick for this example, and it should work for what I'm actually trying to do. This method all makes sense, I'm just not super familiar with working with meshes yet. I understand more now thanks to your answer. One again, thanks a lot! – Aaron Stevens Jul 25 '17 at 17:58
• Something I just noticed with this example is that the new mesh does not include the first point (the point {0,0}), and this happens in my actual problem as well. Is this because we got rid of all mesh elements with that point in it? – Aaron Stevens Jul 27 '17 at 14:13
• Is there any way to get this point back? In my actual problem, points near the origin are very important. – Aaron Stevens Jul 27 '17 at 14:24
• On further inspection, there actually do exist more points to the left of the boundary line from lines in the 1st order mesh that briefly extend into the forbidden region. I think I have enough information to fix this though, or I might find another way around interpolating on these points. Thanks anyway. – Aaron Stevens Jul 27 '17 at 16:16