# How to create a mesh of convex polygons (of maximal size) from a boundary mesh region

I can create a boundary mesh region that consists of a rectangular region with two independently-rotated pentagonal regions excluded using the following code:

d = 10;
x1 = -d/2; y1 = 0; rad1 = 2; \[Theta]1 = 0 Degree;
x2 = d/2; y2 = 0; rad2 = rad1; \[Theta]2 = 20 Degree;
pts1 = CirclePoints[{x1, y1}, {rad1, (Pi/5 - Pi/2) + \[Theta]1}, 5];
pts2 = CirclePoints[{x2, y2}, {rad2, (Pi/5 - Pi/2) + \[Theta]2}, 5];
sx = 10; sy = 8;
boundpts = {{-sx, -sy}, {sx, -sy}, {sx, sy}, {-sx, sy}};
bmr = BoundaryMeshRegion[Flatten[{pts1, pts2, boundpts}, 1],
Line[{1, 2, 3, 4, 5, 1}], Line[{6, 7, 8, 9, 10, 6}],
Line[{11, 12, 13, 14, 11}]]


The graphical output is

I can then create a triangle mesh (with as large cell sizes as possible) with the following code

TriangulateMesh[bmr, MaxCellMeasure -> \[Infinity]]


which results in

My question is whether the triangular mesh can be made into a mesh of the largest possible convex polygons (with greater than 3 sides) instead, so that the number of individual cells is minimized, and, if so, how. I would greatly appreciate any help. Thank you very much!

Addendum: Furthermore, is it possible to specify that the outside boundary edges should not be at all split (or, at most, at one point -- say, the middle -- of each boundary edge)?

• Just to check, do you wish to use only one type of polygon in the solution, ie just triangles as opposed to a combination of triangles and quads? Sep 7 '16 at 9:37
• Do you need the exact solution, i.e. absolute smallest number of polygons? Or just something practical that reduces the polygon count a bit? Sep 7 '16 at 10:25
• @Quantum_Oli No, it can be a mixture of polygons, as long as they are convex. Thank you! Sep 7 '16 at 13:21
• @Szabolcs Just something practical that reduces the number of polygons (even from visual inspection of the example I gave in the question, it's clear that some triangles can be joined together to form convex polygons – I just would like to automate that process). Also, I should mention that I'd like to do this for arbitrary rotations of the pentagonal regions and separation distance between those regions. Thank you!! Sep 7 '16 at 13:23

The generated outside points are called Steiner points. You can prohibit the generation of these with

TriangulateMesh[bmr, MaxCellMeasure -> \[Infinity],
MeshQualityGoal -> 1, "SteinerPoints" -> None]


Now, other mesh elements then triangles are not possible currently. However, as Szabolcs points out in his answer you can merge connected triangles. To do so you can use the "ElementConnectivity" property of an ElementMesh. You can do so in the following way:

Needs["NDSolveFEM"]
mesh = ToElementMesh[bmr, MaxCellMeasure -> \[Infinity],
MeshQualityGoal -> 1, "MeshOrder" -> 1, "SteinerPoints" -> None];
mesh["Wireframe"]


And then extract the connectivity:

mesh["ElementConnectivity"]
{{{0, 3, 2}, {10, 16, 1}, {0, 5, 1}, {0, 6, 5}, {0, 3, 4}, {7, 15,
4}, {0, 12, 6}, {12, 0, 16}, {11, 0, 10}, {2, 0, 9}, {0, 14,
9}, {7, 8, 0}, {15, 0, 14}, {0, 11, 13}, {0, 6, 13}, {8, 2, 0}}}


Which means that the first element has a connection to the boundary (0) and to elements 3 and 2, the second element is connected to elements 10, 16 and 1 and so on. For more details please see the above mentioned documentation that has more examples.

Should you write such an element merging code, I'd be curious to see it.

• Thank you for the great tip about excluding the Steiner points -- that has helped reduce the number of triangles produced by TriangulateMesh, at least! I am still thinking about the element merging problem, but it might not be as pressing an issue now that the extra points/vertices have been eliminated (and the former seems to be rather nontrivial). Sep 7 '16 at 19:29

We can use TriangleLink directly and instruct it not to subdivide any of the edges. I think that TriangulateMesh also uses TriangleLink internally for the 2D case.

The following code relies on the variables defined in the main question.

<<TriangleLink

tri = TriangleCreate[]

TriangleSetPoints[tri, Catenate[{pts1, pts2, boundpts}]]

TriangleSetSegments[tri,
Catenate[Partition[#, 2, 1] & /@ {{1, 2, 3, 4, 5, 1}, {6, 7, 8, 9, 10, 6}, {11, 12, 13, 14, 11}}]]

TriangleSetHoles[tri, Mean /@ N@{pts1, pts2}] (* each of the two points is in one of the holes *)

res = TriangleTriangulate[tri, "pYY"] (* triangulate with appropriate options *)

Graphics[{GraphicsComplex[
TriangleGetPoints[res], {FaceForm[LightGray], EdgeForm[Black],
Polygon@TriangleGetElements[res]}]}]


The options are described in the TriangleTriangulate doc page, under Details. Y means do not divide the boundary edges. YY means do not divide boundary or inner edges. You can use q` to triangles with small angles, but it will result in more triangles being created. Additional documentation on the options is here.

A practical but imperfect way to reduce the polygon count may be to merge neighbouring triangles greedily until it is no longer possible to get a convex merger. Then move on to the next triangle and merge that one with its neighbours. Iterate until no convex merger is possible at all.

• Thank you for the response! This seems very interesting – I will play around with it when I am at my computer. Also, regarding your last paragraph suggesting an iterative combination algorithm, is there an easy way to check for "convexity" in Mathematica? Thanks! Sep 7 '16 at 13:29
• @PhysicsCodingEnthusiast Not sure if there's a fast builtin way. Could make a separate question (if you ask, pls also mention performance). How about verifying that the angle between consecutive line segments (polygon edges) is always < 180 deg? Or rather: you are always rotating in the same direction when moving from segment to segment. Sep 7 '16 at 13:54