3
$\begingroup$

I am generating a refined icosahedron as follows:

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

The output is:

enter image description here

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.

$\endgroup$
1
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$
    – bbgodfrey
    Jun 16 at 17:35
4
$\begingroup$

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]
$\endgroup$
3
  • 1
    $\begingroup$ thank you for your detailed answer; I must confess though I'm struggling with understanding it. Could be a version mismatch, but I'm getting an error from MapThread: MapThread::mptd: Object normals at position {2, 3} in MapThread[Association[polygon->#1,edges->Transpose[{#2,RotateLeft[Slot[<<1>>]]}],normal->#3]&, ... has only 0 of required 1 dimensions $\endgroup$
    – Botond
    Jun 16 at 18:36
  • 1
    $\begingroup$ Oh, I forgot to copy and paste the line that computes the normals. It is now edited. If you are having difficulty understand bits of the code, don't hesitate to ask. I can add more verbage. $\endgroup$ Jun 16 at 18:57
  • $\begingroup$ thank you, @Craig Carter. I'm still digesting all the details, but tested it, runs and works and gives the correct result. $\endgroup$
    – Botond
    Jun 16 at 19:43
4
$\begingroup$

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]

enter image description here

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

enter image description here

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

>     <|1 -> 0.251943, 2 -> 0.198251, 3 -> 0.251943|>

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"]

enter image description here

edgeDihedrals @ BoundaryDiscretizeGraphics @ PolyhedronData["Tetrahedron"]

enter image description here

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"]

enter image description here

{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)"
$\endgroup$
1
  • $\begingroup$ thank you, it works! $\endgroup$
    – Botond
    Jun 16 at 21:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.