# How to build a neighbor table for the Hexagonal lattice

I want to consider a collection of vertices arranged in a finite Hexagonal lattice, with say $$n_{r}$$ rows and $$n_{x}$$ vertices per row, for a total of $$N = n_{x}n_{r}$$ vertices. The goal is to construct a neighbor table, which is a matrix $$A$$ with dimensionality $$N\times N$$. Each element $$A_{ij}$$ is one if there is a bond connecting vertices $$i$$ and $$j$$, and zero otherwise. This can also be seen as an adjacency matrix for a particular graph.

I can't think of a way to build this matrix that isn't frustrating to code as a function of $$n_{r}$$ and $$n_{x}$$, but I'm aware there are some lattice functions built into Mathematica. Maybe those can make the process more smooth if anyone has suggestions.

Also, to make things a little harder, I also want the option of including periodic boundary conditions. This just means if you walk off the finite lattice in a particular direction (say, off the left-hand side) you reappear on the opposite side of the lattice (in this case the right hand side). This means there are new connections that would look long-ranged in any planar representation of the graph, or you can think of it as putting the graph on a torus. Help with that case would be particularly appreciated.

I will take this opportunity to showcase the abilities of IGraph/M for lattice generation and mesh / graph / matrix conversions. IGraph/M thrives on user feedback, so if you find it useful, please take some time to write a few comments about your experience. It will help me improve the package.

### Non-periodic case

You can directly generate a (non-periodic) lattice with IGraph/M.

<< IGraphM

mesh = IGLatticeMesh["Hexagonal", Polygon@CirclePoints[3, 6],
MeshCellLabel -> {2 -> "Index"}]


The second argument of IGLatticeMesh may be a region. This region will be filled with cells. In this case, I chose a big hexagon to be filled with small hexagonal cells.

am = IGMeshCellAdjacencyMatrix[mesh, 2]


"2" means 2-dimensional cells, i.e. little hexagons. "1" would mean edges and "0" points.

MatrixPlot[am]


If you need the graph,

graph = IGMeshCellAdjacencyGraph[mesh, 2,
VertexCoordinates -> Automatic]


Notice that this is actually a triangular connectivity, which could also be generated directly (in some shapes) with IGTriangularLattice. Demo:

{IGTriangularLattice[4], IGTriangularLattice[{3, 5}]}


We could have use IGLatticeMesh too:

IGLatticeMesh["Triangular", {3, 3}]


Let's get the point-to-point connectivity now (instead of the cell-to-cell one):

IGMeshCellAdjacencyGraph[%, 0]


### Periodic case

Now let us do the periodic case.

We start with a hex lattice arranged in an $$n\times m$$ grid.

{n, m} = {5, 6};
mesh = IGLatticeMesh["Hexagonal", {n, m}, MeshCellLabel -> {2 -> "Index"}]


Convert it to a graph. This time I will not preserve the vertex coordinates so that we can get a clearer layout after we make the lattice periodic.

graph = IGMeshCellAdjacencyGraph[mesh, 2, VertexLabels -> "Name"];
graph = VertexReplace[graph, {2, i_} :> i]


I have also converted the vertex names, which were of the form {2, index} (2 indicating two-dimensional mesh cells) to simply index.

We add the extra edges needed for periodic boundary conditions.

extraEdges = DeleteDuplicates@Flatten@Table[
{If[Mod[i, m] == 0, {i <-> i - m + 1, i <-> Mod[i - 2 m + 1, m n, 1]}, {}],
If[i <= m, {i <-> i + m n - m, i <-> Mod[i + m n - m + 1, m n, 1]}, {}]},
{i, m n}
]



Then we can get (or plot) the graph's adjacency matrix.

IGAdjacencyMatrixPlot[pgraph]


am = AdjacencyMatrix[pgraph]


Extra visualization: here's the graph in 3D with {m,n} = {10,20}:

(* remember to re-create graph and extraEdges after setting {m,n} *)
pgraph = Graph3D[EdgeAdd[graph, extraEdges], VertexLabels -> None]


### An alternative solution for the periodic case

The adjacency relations of hexagonal cells form a triangular lattice. There is a function in IGraph/M for directly generating a triangular lattice graph, and it has an option to make it periodic:

IGTriangularLattice[{5, 10}]


IGTriangularLattice[{5, 10}, "Periodic" -> True]


Then you can just get the adjacency matrix again.

Note that the {m,n} syntax in IGLatticeMesh and IGTriangularLattice do not have the exact same meaning—pay attention to the difference if you mix these approaches! The vertex labelling will also be different. Presumably, at some point you will want to use the visualization of the hex lattice mesh to plot your results. Thus it is useful to be able to map back to mesh cell indices.

### Update: How to do this for the connectivity of vertices in a hexagonal graph?

OP is asking how to do this if the vertices of the graph are the vertices (not faces) of the hexagonal mesh.

The simplest way is to use the same method as above, but start with the dual lattice of the hexagonal one, i.e. a triangular lattice.

IGLatticeMesh["Triangular", {4, 5}]


IGMeshCellAdjacencyGraph[triMesh, 2, VertexCoordinates -> Automatic]


We can also do it directly with the vertices of a hexagonal lattice, but it is a bit more trouble because of those two hanging out points you can see in the graph above.

Let us start by creating the graph directly from a hexagonal mesh.

{n, m} = {4, 5};
graph = IGMeshGraph[
IGLatticeMesh["Hexagonal", {n, m}],
VertexShapeFunction -> "Name",
PerformanceGoal -> "Quality"
]


Now we need to add periodicity. This time, I am not going to add extra edges to connect the left and right, top and bottom of the lattice. If we simply repeat this partial lattice in both directions to see which node would need to be connected to which other one, we will immediately see that it is not enough to add connections. It would also be necessary to add two new vertices (red dots in the illustration below).

We are going to merge corresponding vertices at bottom and top, left and right of the lattice. The formulas for correspondences are easy to figure out by making drawings like the one above. For convenience, we will use VertexReplace instead of VertexContract.

bottom = Range[m + 1, 2 n (m + 1), m + 1];
repl1 = Thread[bottom + m -> bottom]
(* {11 -> 6, 17 -> 12, 23 -> 18, 29 -> 24, 35 -> 30, 41 -> 36, 47 -> 42, 53 -> 48} *)

left = Range[1, 2 m];
repl2 = Thread[left + 2 n (m + 1) -> left]
(* {49 -> 1, 50 -> 2, 51 -> 3, 52 -> 4, 53 -> 5, 54 -> 6, 55 -> 7, 56 -> 8, 57 -> 9, 58 -> 10} *)


If you look carefully at the replacement lists, you will notice that we are not done yet. I kept the output for this specific size of lattice so you can see that vertex 53 is replaced by 48 in the top -> bottom replacement and the same vertex 53 is replaced by 5 in the right -> left replacement. This creates an inconsistency. To get the correct result, we also need to merge 5 and 48 in a third step.

repl3 = {2 n (m + 1) -> m}
(* {48 -> 5} *)


The replacement lists must be applied successively and in the correct order, rather than concurrently, because of repeated treatment of the same vertices. We use Fold for this.

pgraph = SimpleGraph@Fold[VertexReplace, graph, {repl1, repl2, repl3}]


In version 11.3, the vertex coordinates are lost in this process. Let us re-add them so we can see the result better, and we can verify that it is correct.

coord = AssociationThread[VertexList[graph], GraphEmbedding[graph]];

pgraph = Graph[pgraph,
VertexCoordinates -> Normal@KeyTake[coord, VertexList[pgraph]],
VertexShapeFunction -> "Name", PerformanceGoal -> "Quality"
]


Notice that with this layout, 5 and 46 are the two vertices that would have been missing if we naively repeat the lattice in every direction and try to add edges (instead of contracting vertices).

I was still not completely confident about the result. As you can see from the necessity for repl3, it is easy to make mistakes. Thus, let us make further checks. We expect the result to be vertex transitive. That means that for any two vertices, the graph has a symmetry which transforms them into each other. Loosely speaking, all vertices look the same, they cannot be distinguished based on their position in the graph (at least not without a reference point).

IGraph/M has a function for this.

IGVertexTransitiveQ[pgraph]
(* True *)


Are all edges interchangeable as well? That is not the case. Clearly, we have three categories of edges, running in three different directions in the geometrically laid out lattice.

To show this, let us make a function that categorizes edges based whether they may be transformed into each other by any graph automorphisms.

edgeCategory[graph_] := With[{lg = LineGraph[graph]},
IGPartitionsToMembership[lg]@
GroupOrbits@PermutationGroup@IGBlissAutomorphismGroup[lg]
]


This function returns a category number for each edge, in the same order as EdgeList.

We can use these numbers for colouring:

Graph[pgraph, EdgeStyle -> Thick] //
IGEdgeMap[ColorData[100], EdgeStyle -> edgeCategory]


Again, everything looks good. Every vertex is incident to three edges of distinct categories, and there are precisely three categories.

pgraph has the symmetries we expect for an infinite hexagonal lattice.

Just for fun, here's a force directed layout for a $$12\times 16$$ periodic lattice.

• This is incredible, thank you so much. I will get a lot out of this package I'm sure - and don't worry I will leave comments when I'm done! Commented Oct 28, 2018 at 20:24
• Hi Szabolcs - sorry to bother you for more details, but I noticed your approach to the periodic case assumes that the fundamental units of the problem are the hexagons. That is, when you draw the extra edges for the periodic case, you connect the faces to those on the opposite edge of the lattice. I actually want to do this with the vertices, so that I have one node per vertex, not face, and walking off the edge of the lattice in one direction connects me to a vertex on the opposite side. Is there an easy way to add edges for that purpose? I think there must be but I'm not seeing it. Commented Oct 31, 2018 at 7:22
• @mflynn I updated the post. Commented Oct 31, 2018 at 9:46
• @Szabolcs, this seems like an excellent package for helping to solve my still-question: mathematica.stackexchange.com/q/159311/7936 Commented Oct 31, 2018 at 16:01
• @Szabolcs: Does Mathematica 12 contains an in build function with which I can calculate the points of a triangular lattice?
– lio
Commented May 14, 2019 at 16:32

## Edges

This problem can be handled elegantly and efficiently using ListCorrelate.

ntab[r_, x_, pad_: "Cyclic"] := (
2^Partition[Range[r x - 1, 0, -1], x]
// ListCorrelate[1 - IdentityMatrix[3], #, 2, pad] &
// IntegerDigits[Join @@ #, 2, r x] &
)


With and without wrapping:

ntab[5, 6] // MatrixPlot

ntab[5, 6, 0] // MatrixPlot


## Vertices

Regarding your comment to Szabolcs the same methods can be applied to a vertex graph.

Now much faster and using far less memory via SparseArray rather than powers-of-two.

ntabV[n_, m_] :=
Module[{r = m + 1, x = 2 (n + 1), a, k},
a = Partition[Hold /@ Range[r x], r];
k[1] = {{1, 1}, {1, 0}};
k[2] = {{0, 1}, {1, 1}};
ListCorrelate[k[#], a[[# ;; ;; 2]], 2 (-1)^#] & /@ {2, 1}
// Thread[{Level[Riffle @@ #, {-1}]}] &
// SparseArray[Automatic, {r x, r x}, 0,
{1, {Range[0, 3 r x, 3], #}, ConstantArray[1, 3 r x]}] &
]



Borrowing Szabolcs's illustration we can label vertices like this:

From the adjacency table output of ntabV we can create a matching explicit list:

ntabV[4, 5]["AdjacencyLists"];
Thread[Range[0, 59] -> (% - 1)]

{0 -> {6, 11, 54}, 1 -> {6, 7, 55}, 2 -> {7, 8, 56}, 3 -> {8, 9, 57},
4 -> {9, 10, 58}, 5 -> {10, 11, 59}, 6 -> {0, 1, 12}, 7 -> {1, 2, 13},
8 -> {2, 3, 14}, 9 -> {3, 4, 15}, 10 -> {4, 5, 16}, 11 -> {0, 5, 17},
12 -> {6, 18, 23}, 13 -> {7, 18, 19}, 14 -> {8, 19, 20}, 15 -> {9, 20, 21},
16 -> {10, 21, 22}, 17 -> {11, 22, 23}, 18 -> {12, 13, 24}, 19 -> {13, 14, 25},
20 -> {14, 15, 26}, 21 -> {15, 16, 27}, 22 -> {16, 17, 28}, 23 -> {12, 17, 29},
24 -> {18, 30, 35}, 25 -> {19, 30, 31}, 26 -> {20, 31, 32}, 27 -> {21, 32, 33},
28 -> {22, 33, 34}, 29 -> {23, 34, 35}, 30 -> {24, 25, 36}, 31 -> {25, 26, 37},
32 -> {26, 27, 38}, 33 -> {27, 28, 39}, 34 -> {28, 29, 40}, 35 -> {24, 29, 41},
36 -> {30, 42, 47}, 37 -> {31, 42, 43}, 38 -> {32, 43, 44}, 39 -> {33, 44, 45},
40 -> {34, 45, 46}, 41 -> {35, 46, 47}, 42 -> {36, 37, 48}, 43 -> {37, 38, 49},
44 -> {38, 39, 50}, 45 -> {39, 40, 51}, 46 -> {40, 41, 52}, 47 -> {36, 41, 53},
48 -> {42, 54, 59}, 49 -> {43, 54, 55}, 50 -> {44, 55, 56}, 51 -> {45, 56, 57},
52 -> {46, 57, 58}, 53 -> {47, 58, 59}, 54 -> {0, 48, 49}, 55 -> {1, 49, 50},
56 -> {2, 50, 51}, 57 -> {3, 51, 52}, 58 -> {4, 52, 53}, 59 -> {5, 48, 53}}


Observe that the seams wrap, e.g.

 2  -> { 7,  8, 56}
24 -> {18, 30, 35}
35 -> {24, 29, 41}
59 -> { 5, 48, 53}


# Explanation

Szabolcs implied that I need to explain this better. Let's start with a drawing of the lattice:

hex = {Polygon@CirclePoints[#, {1.1, 90 °}, 6], Yellow, Text[i++, #]} &;
i = 1;
Array[hex@{2 #2 + #, -Sqrt[3] #} &, {5, 6}] // Graphics


We can see that this is a skewed rectangular matrix, equivalent to:

MatrixForm[m = Partition[HoldForm /@ Range[5*6], 6]]


$$\left( \begin{array}{cccccc} 1 & 2 & 3 & 4 & 5 & 6 \\ 7 & 8 & 9 & 10 & 11 & 12 \\ 13 & 14 & 15 & 16 & 17 & 18 \\ 19 & 20 & 21 & 22 & 23 & 24 \\ 25 & 26 & 27 & 28 & 29 & 30 \\ \end{array} \right)$$

In that rectangular matrix the neighbors of 8 are {1, 2, 3, 7, 9, 13, 14, 15} but in this skewed version 1 and 15 are too far. We eliminate these, as well as the center, using zeros in the "mask" that is our convolution kernel:

k = {{0,  1,  1},            (* visually skewed to make the hexagon apparent *)
{1,  0,  1},
{1,  1,  0}};


We apply this kernel using ListCorrelate:

nt = ListCorrelate[k, m, 2];


The neighbors for 8:

nt[[2, 2]]

2 + 3 + 7 + 9 + 13 + 14


Note also that by default ListCorrelate wraps around:

nt[[1, 1]]  (* neighbors of 1 *)

2 + 6 + 7 + 12 + 25 + 26


(ListCorrelate has a fourth parameter that controls padding; if 0 is specified it effectively does not wrap.)

In this example I used HoldForm so that the integers would not sum. In the complete function ntab I used powers of two to allow them to sum, then get the binary output you wanted using IntegerDigits, e.g.

{2, 6, 7, 12, 25, 26};
2^(% - 1)
Total[%]
IntegerDigits[%, 2, 5*6] // Reverse
Position[%, 1] // Flatten

{2, 32, 64, 2048, 16777216, 33554432}

50333794

{0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, \
0, 0, 0}

{2, 6, 7, 12, 25, 26}


The vertex graph case is similar but I needed two different kernels for the "up" and "down" triads.

• Can you add some explanations to the code, please? Commented Oct 29, 2018 at 7:40
• Okay, now I see how it works. But it is not trivial, with all the bit-position-encodings of the cells. Why not make the answer a bit more educational? Commented Oct 29, 2018 at 7:47
• @Szabolcs Please see my update and tell me if it is sufficient. Commented Oct 29, 2018 at 14:55
• I forgot to respond, but yes, of course it was. Commented Oct 31, 2018 at 17:46
• @Szabolcs Now that I have tuned ntabV I think you might find its performance interesting, e.g. I can run ntabV[300, 300] in under a second, whereas IGLatticeMesh["Hexagonal", {300, 300}]` alone takes almost eleven. Commented Nov 1, 2018 at 0:03