Getting dihedrals of a polyhedron object

Question

I am generating a refined icosahedron as follows:

Needs["PolyhedronOperations`"]
Geodesate[PolyhedronData["Icosahedron"], 3]

The output is:

I can also exract the face coordinates as

faces = Cases[Normal@Geodesate[PolyhedronData["Icosahedron"], 3], _Polygon, Infinity];
f = faces // N;

My question is, how to get all the dihedral angle values between all neighboring faces? I'm using Mathematica 11.3.0.0 and neither the functions DihedralAngle nor PolyhedronAngle seem to be implemented.

Answer 1

I am afraid that one has to resort to extracting the geometry. Here is an example:

Needs["PolyhedronOperations`"]
hedron = Geodesate[PolyhedronData["Icosahedron"], 3]

Extract points.

pts = Cases[hedron, GraphicsComplex[pts_, other__] :> pts][[1]]

Extract indices of polygons

polygons = 
 Cases[hedron, GraphicsComplex[pts_, Polygon[other_]] :> other][[1]]

Write a function to give the normal to the polygon from its points and indices.

normal[pts_, indices_] := 
 With[{a = pts[[indices[[2]]]] - pts[[indices[[1]]]], 
   b = pts[[indices[[3]]]] - pts[[indices[[2]]]]},
  Simplify[Normalize[Cross[a, b]]]
  ]

Compute the normals for each facet (I made this numerical, otherwise the expressions become very long)

normals = normal[N@pts, #] & /@ polygons

Create a handy set of associations for each polygon

polygonAssocs = 
 MapThread[<|"polygon" -> #1, 
    "edges" -> Transpose[{#2, RotateLeft[#2]}], 
    "normal" -> #3 |> &, {Range[Length[polygons]], polygons, normals}]

e.g.

RandomChoice[polygonAssocs]

Make a list of all edges on the polygons, keep only the ordered ones.

allEdges = 
 Select[Flatten[polygonAssocs[[All, "edges"]] , 1], OrderedQ]

A function to see if a polygon contains an edge, either ordered or not. Arrange the pair so that the first polygon is the one containing the oriented edge.

containsEdge[pgonAssoc_, edge : {e1_, e2_}] :=
   With[{polygons = 
    Select[pgonAssoc, MemberQ[#edges, {e1, e2} | {e2, e1}] &]},
  Which[
   MemberQ[polygons[[1, "edges"]], edge], polygons,
   True, Reverse[polygons]
   ]
  ]

For each edge, find the neighboring polygons, and compute the dot-product of their normals.

edgePolygons =
 MapThread[
  With[
    {edge = #1, polygonPair = containsEdge[polygonAssocs, #2]},
    <|"edge" -> edge, "polygons" -> polygonPair, 
     "dihedral" -> 
      ArcCos[
       First[polygonPair]["normal"] . Last[polygonPair]["normal"]]
     |>] &, {Range[Length[allEdges]], allEdges}]

e.g.,

RandomChoice[edgePolygons]

Answer 2

Update: A more direct approach using Region`Mesh`MeshCellNormals and "EdgeFaceConnectivity" to get a Dataset with "edge" (edge index) "faces" (indices of faces connected to "edge"), "facenormals" (normal vectors for "faces"), and "dihedral" (dihedral angle of "faces") as columns:

ClearAll[edgeDihedrals]
edgeDihedrals[bmr_] := Module[{fn = Region`Mesh`MeshCellNormals[bmr, 2]}, 
 Dataset @ MapIndexed[Association @ {"edge" -> #2[[1]], "faces" -> #, 
   "facenormals" -> AssociationThread[#, fn[[#]]], 
   "dihedral" -> ArcCos[Dot @@ fn[[#]]]} &] @ bmr["EdgeFaceConnectivity"]]

Examples:

bdg = BoundaryDiscretizeGraphics @  Geodesate[PolyhedronData["Icosahedron"], 3];

edgeDihedrals[bdg]

edgeDihedrals[bdg][Range[10], {"edge", "faces", "dihedral"}]

edgeDihedrals[bdg][{1, 2, 3}, {"edge", "dihedral"}]

HighlightMesh[bdg, 
  {Style[{1, 1}, Thick, Red], 
   ## & @@ Thread[{2, Normal @ edgeDihedrals[bdg][1, "faces"]}], 
   Style[{1, 15}, Thick, Green], 
   ## & @@ Thread[{2, Normal @ edgeDihedrals[bdg][15, "faces"]}]}, 
 PlotTheme -> "FaceNormals"]

edgeDihedrals @ BoundaryDiscretizeGraphics @ PolyhedronData["Tetrahedron"]

Original answer:

First, use BoundaryDiscretizeGraphics to get a BoundaryMeshRegion object:

bdg = BoundaryDiscretizeGraphics @  Geodesate[PolyhedronData["Icosahedron"], 3];

We can identify faces connected thru an edge using the properties "FaceVertexConnectivity" and "FaceFaceConnectivity". Then, we can use the function Region`Mesh`MeshCellNormals to get normals for each pair of neighboring faces and use Dot + ArcCos to get the dihedral angles.

ClearAll[neighboringFaces, dihedralAngle]
neighboringFaces[bmr_] := Module[{faces = bmr["FaceVertexConnectivity"]}, 
   Association @ MapIndexed[Function[{x, ind}, ind[[1]] -> 
     DeleteCases[x, _?(Length[Intersection[faces[[ind[[1]]]], faces[[#]]]] != 2 &)]], 
    bmr["FaceFaceConnectivity"]]];


dihedralAngle[bmr_][i_, j_] /; MemberQ[neighboringFaces[bmr][i], j] :=
    ArcCos[Dot @@ Region`Mesh`MeshCellNormals[bmr, {{2, i}, {2, j}}]]

Examples:

neighboringFaces[bdg][1]

 {7, 10, 37}

neighboringFaces[bdg][15]

 {3, 18, 100}

HighlightMesh[bdg, {Style[{2, 1}, Red], ## & @@ Thread[{2, neighboringFaces[bdg][1]}], 
  Style[{2, 15}, Green], ## & @@ Thread[{2, neighboringFaces[bdg][15]}]}, 
 PlotTheme -> "FaceNormals"]

{1, #} -> dihedralAngle[bdg][1, #] & /@ neighboringFaces[bdg][1]

 {{1, 7} -> 0.198251, {1, 10} -> 0.251943, {1, 37} -> 0.251943}

{15, #} -> dihedralAngle[bdg][15, #] & /@ neighboringFaces[bdg][15]

 {{15, 3} -> 0.251943, {15, 18} -> 0.198251, {15, 100} -> 0.251943}

$Version

"11.3.0 for Microsoft Windows (64-bit) (March 7, 2018)"