Creating MeshRegions of various 2D lattices

Question

This is a soft question about the best approach to design a function that creates MeshRegion objects of various 2D lattices, e.g. kagome or hexagonal, all through one consistent interface.

A reasonable approach is to provide a unit cell as a Polygon, along with translation vectors. This is what many solutions on this site use, e.g.

• Drawing a trihexagonal tiling
• Making a 2D hexagonal mesh

The problem is that MeshRegion expects a set of points, then a set of polygons that refer to those points by index. Neighbouring cells of the lattice will share some of these points. Thus it is not simple to compute the translation of points: after translating the coordinates, we must check if some of the translated ones coincide with ones we already had. Doing this is also likely to result in a "messy" indexing of the points. Ideally, they should be indexed in a consistent order.

Here's an example to show what I mean, using hexTile:

DiscretizeGraphics[
 Graphics@Flatten@hexTile[3, 5],
 MeshCellLabel -> {0 -> "Index"}
 ]

Other than the messy indexing, I am not happy to use DiscretizeGraphics. It is an extremely general function which must handle lots of different inputs. Such functions in Mathematica are not always very robust. Will it always correctly identify points which are the same? What if the coordinates are floating point? Even ConvexHullMesh won't always detect collinear points, not even with exact coordinates.

How would you approach this problem? There are several potential solutions. I am looking for "clean" ones which perform well and are not too complicated to implement.

Answer 1

I modify hexTile slightly for convenience by skipping the wrapup with Polygon and by Flattening a level.

hexTile1[n_, m_] := 
  With[{hex = 
     Table[{Cos[2. Pi k/6] + #, Sin[2. Pi k/6] + #2}, {k, 6}] &}, 
   Flatten[Table[hex[3 i + 3 ((-1)^j + 1)/4, Sqrt[3]/2 j], {i, n}, {j, m}], 1]];

Let's assume that we have a ragged array polygonpts whose entries are lists containing the vertex positions of each polygon. First, we compute a clustering tolerance ϵ lower than any edge length in the tesselation. Second, we cluster duplicates with the help of Nearest in order to produce a reduced point set newpts. Afterwards, we use the NearestFunction of newpts as a lookup function for the vertex indices.

polygonpts = hexTile1[3, 5];
ϵ = 0.5 Min[Sqrt[Total[(# - RotateLeft[#])^2, {2}]] & /@ polygonpts];
pts = Flatten[polygonpts, 1];
plist = Sort[DeleteDuplicates[
   Compile[{{idx, _Integer, 1}}, 
    Min[idx],
    RuntimeAttributes -> {Listable},
    Parallelization -> True
    ][Nearest[pts -> Automatic, pts, {∞, ϵ}]]
   ]];
newpts = pts[[plist]];
nf = Nearest[newpts -> Automatic];
polygons = Internal`PartitionRagged[Flatten[nf[pts]], Length /@ polygonpts];

MeshRegion[newpts, Polygon[polygons], MeshCellLabel -> {0 -> "Index"}]

It appears that this is the same mesh as the one returned by DiscretizeGraphics.

Sorting the reduced point set first leads to somewhat nicer enumeration of points:

newpts = Sort[pts[[plist]]];
nf = Nearest[newpts -> Automatic];
polygons = Internal`PartitionRagged[Flatten[nf[pts]], Length /@ polygonpts];

MeshRegion[newpts, Polygon[polygons], MeshCellLabel -> {0 -> "Index"}]

Answer 2

A slight modification of hexMesh[] from this answer can generate hexagonal MeshRegion[] objects:

(* https://mathematica.stackexchange.com/questions/975 *)
multisegment[lst_List, scts : {__Integer?Positive}, offset : {__Integer?Positive}] :=
             Module[{n = Length[lst], k, offs},
                    k = Ceiling[n/Mean[offset]];
                    offs = Prepend[Accumulate[PadRight[offset, k, offset]], 0];
                    Take[lst, #] & /@ 
                    TakeWhile[Transpose[{offs + 1, offs + PadRight[scts, k + 1, scts]}],
                              Apply[And, Thread[# <= n]] &]] /;
             Length[scts] == Length[offset];

multisegment[lst_List, scts : {__Integer?Positive}] := 
             multisegment[lst, scts, scts] /; Mod[Length[lst], Total[scts]] == 0

hexMesh[n_Integer, m_Integer] := 
   MeshRegion[Flatten[Delete[NestList[TranslationTransform[{0, Sqrt[3]}], 
                                      FoldList[Plus, {-1, Sqrt[3]}/2,
                                               Table[AngleVector[-π Sin[k π/2]/3],
                                                     {k, 4 n + 1}]], m],
                             {{1, -1}, {-1, 1}}], 1], 
              Polygon[Flatten[{multisegment[#1, {4, 2}, {3, 1}], 
                               Reverse[multisegment[Rest[#2], {2, 4}, {1, 3}], 2]} & @@@ 
                              Partition[Join[{PadRight[Range[4 n + 1], 4 n + 2]}, 
                                             Partition[Range[4 n + 2, m (4 n + 2) - 1],
                                                       4 n + 2],
                                             {PadLeft[m (4 n + 2) - 1 + Range[4 n + 1],
                                                      4 n + 2]}], 2, 1],
                              {{1, 3}, {2, 4}}]]]

hexMesh[3, 4]

Answer 3

The method for generating the 籠目 lattice in this answer can straightforwardly be converted to something that generates MeshRegion[] objects (multisegment[] is defined in my other answer):

kagomeMesh[m_Integer?Positive, n_Integer?Positive] := 
      MeshRegion[Flatten[Table[If[EvenQ[j] && EvenQ[k], Nothing,
                                  {j + (k - 3)/2, (k - 1) Sqrt[3]/2}],
                               {k, 2 n + 1}, {j, 2 m + 1}], 1], 
                 Polygon[Flatten[Apply[{Append[Most[#], First[#2]], 
                                        Flatten[Riffle[#2, {Rest[#], Reverse[Most[#3]]}]], 
                                        Prepend[Rest[#3], Last[#2]]} &, 
                                       Transpose[Partition[MapIndexed[
                         With[{l = Mod[First[#2], 2]}, Partition[#, l + 2, l + 1]] &,
                         Most[multisegment[Range[(n + 1) (3 m + 2)], {2 m + 1, m + 1}]]],
                         3, 2], {1, 3, 2}], {2}], 2]]]

kagomeMesh[5, 3]