Mesh with Prescribed Vertices

Question

Mathematica has documentation for triangulating regions bounded by curves using e.g. TriangulateMesh. This is very useful for the work that I (aspire to) do drawing pictures of complicated surfaces with integral representations. I triangulate a region, solve my problem numerically on the vertices, and then I can make a graph over the triangulated region out of triangles which will closely resemble a solution, if one exists.

But when I want to start this numerical scheme, I will want to use a value that I know as initial data. Unfortunately, Mathematica will not give me a triangulation which is guaranteed to have any specific points as vertices, which is annoying.

Can I use TriangulateMesh to incorporate a given point? For example, suppose I want to triangulate something simple, like a disk. How can I guarantee that the origin appears as a vertex?

There is a related question here, about adding specific edges.. I think this is more complicated than I need, since I am not (at the moment, at least) I am only concerned with adding points, and letting Mathematica draw any edges it wants.

Update, some months later: I am revisiting this question because I have once again encountered this problem, and this time made some progress myself, but the progress has in turn lead to more difficulties.

Suppose I have a MeshRegion called m consisting of vertices, edges, and faces in a planar region. This is the kind of thing naturally made by TriangulateMesh, for example. Here is a way to add a specific point to the region.

First, retrieve the list vertices = Apply[Complex, MeshCoordinates[m], 1];

Define NewVertices=Append[vertices,NewPoint]. Then, find the cell which contains the vertex. In my relatively simple case, I just use Show to see where the point is and find the cell by zooming in. Then make a new MeshRegion by writing:

NewMesh = MeshRegion[ReIm[NewVertices],Join[MeshCells[m, 1], {Line[{453, 645}], Line[{459, 645}], Line[{505, 645}]}]]

Numbers are chosen from my specific example but don't mean anything here. This works great, if you only want a mesh skeleton of dimension $1$. Problem is, if I try to go further and add to this code stuff like ,Join[MeshCells[m,2],{Polygon[{a,b,c}]}, then Mathematica complains and will not do this. In fact, it even complains if I don't try to Join anything at all, and just try to incorporate any of the old faces. In particular, its complaint is nonopt, it expects options where my 2-cells are supposed to go, even though the documentation seems to suggest I can have cells in any dimension.

What is the cause of the error, and how do I fill in my skeleton? I'd like to keep using TriangulateMesh if possible.

Answer 1

I am not sure how to do this with TriangulateMesh but you can do this with the finite element mesher ToElementMesh (See the options for more details on "IncludePoints":

Needs["NDSolve`FEM`"]
includeCoord = N[{{0, 0}, {1/Sqrt[2], 1/2}}];
mesh = ToElementMesh[Disk[], "IncludePoints" -> includeCoord];
Position[mesh["Coordinates"], #] & /@ includeCoord
{{{49}}, {{50}}}

This will generate a second order mesh - MeshRegion and friends understand that but can not really do anything with that. So to generate a first order FEM mesh and convert that to MeshRegion you can use:

Needs["NDSolve`FEM`"]
includeCoord = N[{{0, 0}, {1/Sqrt[2], 1/2}}];
mesh = MeshRegion[
   ToElementMesh[Disk[], "IncludePoints" -> includeCoord, 
    "MeshOrder" -> 1]];
Position[MeshCoordinates[mesh], #] & /@ includeCoord

Answer 2

Actually, I have answered my own question. The problem was a mis-reading of the documentation. To adjoin the cells, do what I do in the question, but when Joining, join them all at once. Join[MeshCells[m,2],{Polygon[{a,b,c}]} belongs in the same set of braces as Join[MeshCells[m, 1], {Line[{453, 645}]. To be specific, the formatting of MeshRegion is NOT MeshRegion[{List of points},{List of edges},{List of faces}, ....]. Rather, it's MeshRegion[{List of points},{List of all cells of any dimension}].

Here is the code which corrects a region. I've taken the unit disk, and seed point the origin, but you can adapt it to any region you like:

`In[1]:= R = DiscretizeRegion[Disk[], MeshCellLabel -> {0 -> "Index"}]  
SeedPoint = {0, 0};
v = Length[MeshCells[R, 0]];`

Next we detect which cell contains the point you want to add, and call that cell j

`For[i = 1, i <= Length[MeshCells[R, 2]], i++,
  If[TestPoint \[Element] ConvexHullRegion[{MeshCoordinates[R] 
     [[MeshCells[R, 2][[i, 1, 1]]]], 
       MeshCoordinates[R][[MeshCells[R, 2][[i, 1, 2]]]], 
       MeshCoordinates[R][[MeshCells[R, 2][[i, 1, 3]]]]}], 
       j = i;
    Print[j]
    ];
   ];`

This might not be optimal, but it seems that Mathematica does not recognize Polygon[{x,y,z}] as a region, it is using the addresses x,y,z in the list of vertices to know how to draw the mesh, but not thinking of it as a subset of the plane. So this code works around by using ConvexHullRegion and then takes the convex hull of the three points that define each triangle. Once we know which triangle it's in, we define a new MeshRegion consisting of the old vertices and the seed point for its list of vertices, the edges of the old triangulation, all of the old faces except the one flagged j, and then we add in an edge from SeedPoint to each of the vertices of j and a face consisting of SeedPoint and two out of the three vertices of j. This is done in the following block:

`a = MeshCells[R, 2][[j, 1]][[1]];
 b = MeshCells[R, 2][[j, 1]][[2]];
 c = MeshCells[R, 2][[j, 1]][[3]];
 CorrectedVertices = Append[MeshCoordinates[R], SeedPoint];
 CorrectedMesh = MeshRegion[CorrectedVertices,
   Join[
   MeshCells[R, 1],
   Complement[MeshCells[R, 2], {MeshCells[R, 2][[j]]}],
   {Line[{v + 1, a}], Line[{v + 1, b}], Line[{v + 1, c}]},
   {Polygon[{v + 1, a, b}], Polygon[{v + 1, a, c}], Polygon[{v + 1, b, c}]}
   ]
 ]`

This code adds only a single vertex. For a fixed finite number of additional vertices, just repeat the process with each subsequent vertex. Lastly I apologize for any issues with formatting, this is my first time submitting code to this site.

Answer 3

We can identify the polygon (p) that is nearest the added vertex (v) using the function Region`Mesh`MeshNearestCellIndex and replace p with triangles formed by v and successive pairs of vertices of p.

ClearAll[addVertex]
addVertex[mesh_, pt : {__?NumericQ}, threshold_: 10^-3] :=
 Module[{mc = Append[pt]@MeshCoordinates[mesh],
   nci = Region`Mesh`MeshNearestCellIndex[mesh]@pt},
   If[RegionDistance[RegionBoundary@MeshPrimitives[mesh, nci]]@pt < threshold,
   mesh,
   MeshRegion[mc,
    {DeleteCases[MeshCells[mesh, nci]]@MeshCells[mesh, 2],
     Polygon[Append[Length@mc]@#] & /@ Partition[First@MeshCells[mesh, nci], 2, 1, 1]}]]]

For the case of multiple new vertices, Fold the function addVertex over the new vertices:

addVertex[mesh_, pts : {{__?NumericQ} ..}] := Fold[addVertex, mesh, pts]

Examples:

region = DiscretizeRegion[Disk[], MaxCellMeasure -> .2, ImageSize -> 500];

newpoints = {{0, 0}, {.3, .5}};

Row[{Show[region, Graphics[{Red, PointSize @ Medium, Point@newpoints}]],
   Show[addVertex[region, newpoints], 
     Graphics[{Red, PointSize@Medium, Point@newpoints}], ImageSize -> 500],
   Show[HighlightMesh[addVertex[region, newpoints], 
     Style[#, EdgeForm[White], Opacity[.5], RandomColor[]] & /@ 
      Complement[Flatten@MeshCells[addVertex[region, newpoints], 2], 
        Flatten@MeshCells[region, 2]], ImageSize -> 500], 
   MeshRegion[MeshCoordinates@region, MeshCells[region, All], 
    PlotTheme -> "Lines", 
    MeshCellStyle -> {{1, All} -> Directive[Opacity[1], Black, Thick]}], 
   Graphics[{Red, AbsolutePointSize@5, Point@newpoints}]]}]

Use

SeedRandom[7]
points = RandomReal[{0, 10}, {12, 2}];
region = VoronoiMesh[points, ImageSize -> 500];
newpoints = RandomPoint[vor, 7];

to get

Note: If the new vertex v is close to the boundary of the polygon p nearest v (that is, its distance from the boundary of p is less than optional third argument threshold), then p is not touched.