Generating icosahedra with any T number

Results

12.1.0 for Linux x86 (64-bit) (March 14, 2020)

To begin at the end, this figure shows a $(m,n)=(5,2)$ boundary mesh on a unit sphere. The boundary mesh was generated using functions provided below. The icosahedron inside the sphere is there to orient the viewer. The icosahedron is slightly smaller than the sphere, so its vertices are not on the sphere.

The Code

The code provided in this answer is a re-write of an earlier answer, which is now deleted. All of the code required to generate the boundary mesh and the figure above is provided in this section. There are two supporting functions. An example of their usage and an explanation of the code are given below.

Generator Function

This function returns an association of functions. Each function generates a different set of points, as described below. One of the functions is used to generate some the mesh points shown in the figure above.

ClearAll[createGeneratorFunction] createGeneratorFunction[{m_Integer, n_Integer}] := Block[{ origin, reg, uDir, vDir, allPts, boundary, outer, vertices = ConstantArray[False, (m^2 + 3 m + 3 n + n^2)/2 + 2 m *n + 1]}, allPts = origin + Flatten[Join[ Table[(k - irow)*uDir + irow*vDir, {irow, 0, n}, {k, 0, m + irow}], Table[(k - n) uDir + (m + n - irow) vDir, {irow, 0, m - 1}, {k, 0, n + irow}]], 1]; boundary = With[{org = {0, 0}, aDir = {1, 0}, bDir = {Cos[\[Pi]/3], Sin[\[Pi]/3]}, line = Line[{origin, m*uDir + n*vDir, (m + n)*vDir - n*uDir, origin}]}, (Element[#, line] /. {origin -> org, uDir -> aDir, vDir -> bDir}) & /@ allPts ]; outer = Not /@ With[{org = {0, 0}, aDir = {1, 0}, bDir = {Cos[\[Pi]/3], Sin[\[Pi]/3]}}, reg = Polygon[{org, org + m*aDir + n*bDir, org + (n + m)*bDir - n*aDir}]; (# /. {origin -> org, uDir -> aDir, vDir -> bDir}) \[Element] reg & /@ allPts]; vertices[[{1, m*(n + 1) + (n + 2) (n + 1)/2, m*(n + 1) + (n + 2) (n + 1)/2 + 1}]] = True; Association[ "AllPoints" -> Function[{origin, uDir, vDir}, Evaluate[allPts]], "BoundaryPoints" -> Function[{origin, uDir, vDir}, Evaluate[Pick[allPts, boundary]]], "ExteriorPoints" -> Function[{origin, uDir, vDir}, Evaluate[Pick[allPts, outer]]], "InteriorPoints" -> Function[{origin, uDir, vDir}, Evaluate[Pick[allPts, Not /@ Thread[Or[outer, boundary]]]]], "IntervertexPoints" -> Function[{origin, uDir, vDir}, Evaluate[Pick[allPts, Thread[And[boundary, Not /@ vertices]]]]], "Vertices" -> Function[{origin, uDir, vDir}, Evaluate[Pick[allPts, vertices]]] ] ] 

Geodesate Function

This function returns the coordinates of the mesh points on the unit sphere.

ClearAll[geodesate] geodesate[{m_Integer, n_Integer}, polyhedron_String] := Block[ {polygons = Cases[PolyhedronData[polyhedron, "Polygons"], _Polygon, ∞], vertices = PolyhedronData[polyhedron, "Vertices"], regNear = RegionNearest[Sphere[PolyhedronData[polyhedron, "Midcenter"]]], betweenPts = Subdivide[First[#], Last[#], GCD[m, n]][[2 ;; -2]] & , edgeEndpts = PolyhedronData[polyhedron, "Vertices"][[#]] & /@ PolyhedronData[polyhedron, "Edges"], gen = createGeneratorFunction[{m, n}], lsf = LinearSolve[{{m, 0, 0, n, 0, 0}, {0, m, 0, 0, n, 0}, {0, 0, m, 0, 0, n}, {-n, 0, 0, m + n, 0, 0}, {0, -n, 0, 0, m + n, 0}, {0, 0, -n, 0, 0, m + n}}], coords, pts, uDir, vDir}, coords = {{PolyhedronData[polyhedron, "Vertices"]}, betweenPts /@ edgeEndpts, Table[pts = First@poly; {uDir, vDir} = ArrayReshape[ lsf[Flatten[{pts[[2]] - pts[[1]], pts[[3]] - pts[[1]]}]], {2, 3}]; gen["InteriorPoints"][pts[[1]], uDir, vDir], {poly, polygons}]}; regNear /@ Flatten[coords, 2] ] 

The list of point coordinates begins with the coordinates of the 12 vertices of the underlying icosahedron. Next are the points that fall on icosahedron edges between the vertices. Such intervertex points occur only when $m$ and $n$ are relatively prime. The number of intervertex points on each edge is GCD[m,n]-1.

Usage

This code generates a (5,2) boundary mesh based on the icosahedron. This code produces the graphical output shown above.

$Version coords = geodesate[{5, 2}, "Icosahedron"]; dmesh = DelaunayMesh[coords]; bmesh = BoundaryMesh[dmesh, MeshCellStyle -> {{2, All} -> Directive[Opacity[.3], Lighter@Green], {1, All} -> Red, {0, All} -> Black}]; Show[{bmesh, PolyhedronData["Icosahedron"], Graphics3D[{Red, PointSize[0.001], Point@coords}]}, Axes -> False, AxesLabel -> {"X", "Y", "Z"}, Boxed -> False, BoxRatios -> {1, 1, 1}, ImageSize -> Medium, ViewPoint -> 2 {Sin[\[Theta]] Cos[\[Phi]], Sin[\[Theta]] Sin[\[Phi]], Cos[\[Theta]]} /. {\[Theta] -> 7 \[Pi]/36, \[Phi] -> 6 \[Pi]/36}] 

triHex Function

This function and the code following it produce the triaxial hexagon shown in the figure used below, in the Explanation section.

ClearAll[triHex] triHex[{m_, n_}] := Block[ {uDir = {1, 0}, vDir = {Cos[\[Pi]/3], Sin[\[Pi]/3]}, allPts, vertices, boundary, exterior}, allPts = With[{dir = vDir - uDir}, Flatten[Join[ Table[irow*uDir + k*dir, {irow, 0, m}, {k, 0, n + irow}], Table[ m*uDir + irow*vDir + k*dir, {irow, 1, n}, {k, 0, m + n - irow}]], 1]]; vertices = {{0, 0}, m*uDir + n*vDir, (m + n) vDir - n*uDir}; boundary = Select[allPts, Element[#, Line[{vertices[[1]], vertices[[2]], vertices[[3]], vertices[[1]]}]] &]; exterior = Select[allPts, Not[Element[#, Polygon[vertices]]] &]; Association[{"AllPoints" -> allPts, "BoundaryPoints" -> boundary, "ExteriorPoints" -> exterior, "InteriorPoints" -> Complement[allPts, exterior, boundary], "IntervertexPoints" -> Complement[boundary, vertices], "Vertices" -> vertices}] ] thx = triHex[{5, 2}] // N; grfx = Graphics[{PointSize[0.05], Red, Point@thx["Vertices"], PointSize[0.025], Green, Point@thx["InteriorPoints"], PointSize[0.04], Blue, Point@thx["IntervertexPoints"], EdgeForm[{Thick, Magenta}], Opacity[0], Polygon[thx["Vertices"]]}]; Needs["NDSolve`FEM`"]; emesh = ToElementMesh[thx["AllPoints"]] // Quiet; Show[{emesh["Wireframe"], grfx}, ImageSize -> Small]; 

Explanation

An irregular hexagon with three symmetry axes is well-suited to a quick understanding of the code provided above.

The interior of the hexagon has been meshed with equilateral triangles. The triaxial hexagon can be generated by starting at one of the red vertices and going $m=5$ steps in the u-direction, then turning 60-degrees (CCW) and going $n=2$ steps in the v-direction. The next $m$ steps are in the $\vec{v}-\vec{u}$ direction followed by $n$ steps in the $-\vec{u}$ direction. And so on.

In this 2D example, all of the grid points can be easily generated from unit vectors $\vec{u}=(1,0)$ and $\vec{v}=(\cos\frac{\pi}{3},\sin\frac{\pi}{3})$, which is exactly what the triHex function does. This function is only used to produce the triaxial hexagon figure. It is not used to produce the spherical mesh. However, the triHex function is also a prototype for the function that is used to generate the points in the interior of the large magenta equilateral triangle circumscribed by the triaxial hexagon.

The magenta triangle is similar to the triangular faces of the icosahedron. Further, we can easily find the points that are in the interior of the magenta triangle, especially in 2D. The interior points, shown in green in the figure above, can be mapped to the iscosahedron and then to the sphere.

The beauty of it is, the same expression for the interior points in terms of the 2D vectors $\vec{u}$ and $\vec{v}$ can be used with 3D vectors.

triHex returns an association of sets of 2D points that contain only the boundary points, only the interior points, only the intervertex points, etc. These sets of points are replaced by a sets of expressions that can be used in 3D to generate the corresponding points on faces of the undelying icosahedron.

Generator Function

We want a function that will take two appropriate vectors $\vec{u}$ and $\vec{v}$ as arguments and return a list of the coordinates of the interior points. The number of points and the expression for each of the points in terms of $\vec{u}$ and $\vec{v}$ depend only upon $(m,n)$. So, given $(m,n)$ the function createGeneratorFunction returns a function, which maps $\vec{u}$ and $\vec{v}$ to the points. createGeneratorFunction is similar to the triHex function. Analogously to triHex, createGeneratorFunction generates expressions for each point in the triaxial hexagon. Then it evaluates the expressions with 2D basis vectors to determine which expressions yield points in the interior of the magenta triangle. It then returns those expressions as functions of an origin point and two basis vectors. createGeneratorFunction returns an association that can be used to select the generator function appropriate the set of points to be generated.

The same generator function can be used with different basis vectors for each face of the icosahedron. The basis vectors are given as the (2D or 3D) coordinates of the triangle corresponding to a face of the icosahedron. Here is a code snippet that shows the usage of createGeneratorFunction to generate an $m=12$, $n=8$ hexagon with color-coded point sets.

With[{gen = createGeneratorFunction[{12, 8}], face = {{0, 0}, {1, 0}, {Cos[\[Pi]/3], Sin[\[Pi]/3]}}}, Graphics[{PointSize[0.01], Point[gen["ExteriorPoints"] @@ face], PointSize[0.04], Red, Point[gen["Vertices"] @@ face], PointSize[0.04], Blue, Point[gen["IntervertexPoints"] @@ face], PointSize[0.025], Green, Point[gen["InteriorPoints"] @@ face] }, ImageSize -> Small]] 

Basis vector calculation

The basis vectors for each face of the icosahedron can be calculated from the three vertices of the face. Taking one point as the origin, $P_1$, the second point, $P_2$, is found from $m\vec{u}+n\vec{v}$. The third point, $P_3$, is at $(m+n)\vec{v}-n\vec{u}$ from the first point, as can be seen from the triaxial hexagon diagram.

In other words, there is a linear transformation for the vectors $\vec{P_1P_2}$ and $\vec{P_1P_3}$ from the vectors $\vec{u}$ and $\vec{v}$. The LinearSolve function is used obtain $\vec{u}$ and $\vec{v}$ from the vertices $P_1, P_2, P_3$.

geodesate

The geodesate function takes $(m,n)$ and the string "Icosahedron" as arguments and returns the coordinates of all the points in the boundary mesh of the sphere. The geodesate function uses exact numbers in its calculations.

From the values of $(m,n)$ geodesate first obtains a generator function (from createGeneratorFunction) and a linear solver function from LinearSolve. geodesate first takes the vertex coordinates of the icosahedron. Next, geodesate subdivides the icosahedron edges to generate any intervertex points.

geodesate then loops over the faces of the icosahedron, calculates the basis vectors for each face, and applies the generator function to obtain the interior points for each face.

This process does not produce any duplicate points.

Creating the mesh

The code in the usage section above should be fairly self-explanatory. Once the coordinates are found, DelaunayMesh is used to produce a 3D tetrahedral mesh of the interior of the sphere. Then BoundaryMesh is applied to the Delaunay mesh to obtain the mesh on just the spherical surface.

Remarks

In the above $(m,n)=(5,2)$ example, the number of coordinates is Length[coords] (* 392 *). The number of points in the boundary mesh, however, is MeshCellCount[bmesh,1] (* 1170 *), almost 3 times as many.

geodesate appears to work with $n=0$.