Some Unistable Polyhedra

(*program za subdivizijo vej na d delov, dani so cikli*)
EdgesOnFace[lice_]:=Module[{i},(*iz cikla dobimo veje*)
Append[Table[{lice[[i]],lice[[i+1]]},{i,1,Length[lice]-1}],{Last[lice],First[lice]}]];

EdgesOnFaces[cic_]:=Block[{edgs={},n=Length[cic],ed,i},(*vse veje iz ciklov*)
		Do[ed=EdgesOnFace[cic[[i]]];edgs=Union[edgs,ed],
			{i,1,n}];
		edgs=Union[Map[Sort,edgs,1]]];

f1[g_][{i_,j_}]:=Norm[g[i]-g[j]];
f2[g_][{i_,j_}]:=Norm[g[i]-g[j]](g[i]+g[j])/2;
volume1[telo_]:=Apply[Plus,Map[f1[telo[[2,#]]&],EdgesOnFaces[telo[[1]]]]];
moment1[telo_]:=Apply[Plus,Map[f2[telo[[2,#]]&],EdgesOnFaces[telo[[1]]]]];

volume0[telo_]:=Length[telo[[2]]];

moment0[telo_]:=Apply[Plus,telo[[2]]];

volume3[solid_]:=Module[{f,v,vol=0,faces,vertices,face,n},faces=solid[[1]];f=Length[faces];vertices=solid[[2]];
Do[n=Length[faces[[k]]];face=faces[[k]];Do[vol=vol+Det[{vertices[[face[[1]]]],vertices[[face[[j]]]],vertices[[face[[j+1]]]]}] ,{j,2,n-1}],
			{k,1,f}];vol/6];

moment3[solid_]:=Module[{f,v,vol=0,volz=0,faces,zt,face,vertices},faces=solid[[1]];f=Length[faces];vertices=solid[[2]];
Do[n=Length[faces[[k]]];face=faces[[k]];zt=aver[Map[solid[[2,#]]&,face]][[1]];vol=0;Do[vol=vol+1/6Det[{vertices[[face[[1]]]],vertices[[face[[j]]]],vertices[[face[[j+1]]]]}] ({vertices[[face[[1]]]]+vertices[[face[[j]]]]+vertices[[face[[j+1]]]]}/4),{j,2,n-1}];volz=volz+vol,
			{k,1,f}];volz[[1]]];

volume2[solid_]:=Module[{f,v,vol=0,faces,vertices,face,n},faces=solid[[1]];f=Length[faces];vertices=solid[[2]];
Do[n=Length[faces[[k]]];face=faces[[k]];Do[vol=vol+1/2 Norm[Cross[vertices[[face[[j]]]]-vertices[[face[[1]]]],vertices[[face[[j+1]]]]-vertices[[face[[1]]]]] ] ,{j,2,n-1}],
			{k,1,f}];vol];

moment2[solid_]:=Module[{f,v,vol=0,volz=0,faces,zt,face,vertices},faces=solid[[1]];f=Length[faces];vertices=solid[[2]];
Do[n=Length[faces[[k]]];face=faces[[k]];zt=aver[Map[solid[[2,#]]&,face]][[1]];vol=0;Do[vol=vol+1/6 Norm[Cross[vertices[[face[[j]]]]-vertices[[face[[1]]]],vertices[[face[[j+1]]]]-vertices[[face[[1]]]]] ]({vertices[[face[[1]]]]+vertices[[face[[j]]]]+vertices[[face[[j+1]]]]}),{j,2,n-1}];volz=volz+vol,
			{k,1,f}];volz[[1]]];

aver[list_]:=Apply[Plus,list]/Length[list];

volumes={volume0,volume1,volume2,volume3};
moments={moment0,moment1,moment2,moment3};

solid[m_,ep_,n1_,k1_]:=Module[{k,fi,delta,r1,r2,s11,s21,zac,cikli,kor1,del=(Pi-2(Pi/m+ep))/(m-2)},fi=Join[{0},Table[Pi/m+ep+i del,{i,0,m-2}],{Pi}];delta=Join[{0,Pi/m+ep},Table[del,{m-2}],{Pi/m+ep}];r1[1]=3;
r1[n_]:=r1[n-1]/Cos[delta[[n]]]+0.0//N;r2=Table[r1[n],{n,1,m+1}];s11=Graphics[{{Thickness[0.01],Line[Table[{Re[r2[[n]]Exp[I(-Pi/2+fi[[n]])]],Im[r2[[n]]Exp[I(-Pi/2+fi[[n]])]]},{n,1,m+1}]]},{Thickness[0.001],Table[Line[{{0,0},{Re[r2[[n]]Exp[I(-Pi/2+fi[[n]])]],Im[r2[[n]]Exp[I(-Pi/2+fi[[n]])]]}}],{n,1,m+1}]}}];s21=Graphics[{Thickness[0.01],Line[Table[{Re[r2[[n]]Exp[I(-Pi/2-fi[[n]])]],Im[r2[[n]]Exp[I(-Pi/2-fi[[n]])]]},{n,1,m+1}]]}];
k=(n1-k1)/r2[[m+1]];
zac=Table[{n,If[n==2m-1,1,n+1],If[n==2m-1,1,n+1]+2m-1,n+2m-1},{n,1,2m-1}];
cikli=Join[zac,{Reverse[Range[2m-1]]},{Range[2m,4m-2]}];
kor1=Join[Table[{Re[r2[[n]]Exp[I(-Pi/2+fi[[n]])]],Im[r2[[n]]Exp[I(-Pi/2+fi[[n]])]],k  Im[r2[[n]]Exp[I(-Pi/2+fi[[n]])]]-n1},{n,2,m+1}],Reverse[Table[{Re[r2[[n]]Exp[I(-Pi/2-fi[[n]])]],Im[r2[[n]]Exp[I(-Pi/2-fi[[n]])]],k  Im[r2[[n]]Exp[I(-Pi/2-fi[[n]])]]-n1},{n,2,m}]],Table[{Re[r2[[n]]Exp[I(-Pi/2+fi[[n]])]],Im[r2[[n]]Exp[I(-Pi/2+fi[[n]])]],-k  Im[r2[[n]]Exp[I(-Pi/2+fi[[n]])]]+n1},{n,2,m+1}],Reverse[Table[{Re[r2[[n]]Exp[I(-Pi/2-fi[[n]])]],Im[r2[[n]]Exp[I(-Pi/2-fi[[n]])]],-k  Im[r2[[n]]Exp[I(-Pi/2-fi[[n]])]]+n1},{n,2,m}]]];{{s11,s21},{cikli,kor1},r2[[m+1]]}];
(* code for nets*)
AdjacencyPairsHI[lis_] := Cases[Union @@ Table[{a, b}, 
       {a, 1, Length[lis] - 1}, {b, a + 1, Length[lis]}], 
     {x_, y_} /; AdjacentHI[lis][x, y]];
 AdjacentHI[lis_][a_, b_] := Length[Intersection[lis[[a]], lis[[b]]]] == 2;
TreeHI[pove_][n_][m_] := Module[{tr = {pove[[n]]}}, 
     While[ !MaksiTreeHIQ[tr, m], Do[If[DodajHIQ[tr, pove[[k]]], 
         tr = Append[tr, pove[[k]]]; Break[]], {k, 1, Length[pove]}]]; tr];
MaksiTreeHIQ[tr_, n_] := Module[{flag = True}, 
     Do[If[NiNode[i][tr], flag = False; Break[]], {i, 1, n}]; flag];
NiNode[a_][seznam_] := Module[{flag = True}, 
     Do[If[a == seznam[[i]][[1]] || a == seznam[[i]][[2]], 
        flag = False; Break[]], {i, 1, Length[seznam]}]; flag];
DodajHIQ[tr_, {a_, b_}] := JeNode[a][tr] && NiNode[b][tr] || 
     JeNode[b][tr] && NiNode[a][tr];
 JeNode[a_][seznam_] := MemberQ[seznam, {a, x_}] || MemberQ[seznam, {y_, a}];
UrediTree2HI[drevo_] := Module[{drev = drevo, len = Length[drevo]}, 
     Do[Do[If[drev[[k]][[2]] == drev[[i]][[2]], 
         drev = ReplacePart[drev, {drev[[k]][[2]], drev[[k]][[1]]}, k]], 
        {i, 1, k - 1}], {k, 2, len}]; drev];
ChangeOrientation[pari_][lica_] := Module[{nnov, nprej, n, 
      len = Length[lica], ai, zi, indlicap, indlican, licep, licen, presek, 
      p1, p2, q1, q2, a1, ainov, r1, r2, aiprej, aikompl, aizaceten, qq1, 
      uni1, uni2, newfaces = lica, avindex, edges, q1n, q2n}, 
     Do[indlicap = pari[[m]][[1]]; indlican = pari[[m]][[2]]; 
        nnov = Length[lica[[indlican]]]; nprej = Length[lica[[indlicap]]]; 
        licep = newfaces[[indlicap]]; licen = newfaces[[indlican]]; 
        presek = Intersection[licep, licen]; {{p1}} = Position[licep, 
          presek[[1]]]; {{p2}} = Position[licep, presek[[2]]]; 
        {{q1}} = Position[licen, presek[[1]]]; 
        {{q2}} = Position[licen, presek[[2]]]; q1n = q1; q2n = q2; 
        If[PrejQHI[nprej][p1, p2] && PrejQHI[nnov][q1, q2] || 
          PrejQHI[nprej][p2, p1] && PrejQHI[nnov][q2, q1], 
         newfaces = ReplacePart[newfaces, Reverse[licen], indlican]; 
          q1n = nnov - q1; q2n = nnov - q2]; newfaces = ReplacePart[newfaces, 
          If[PrejQHI[nnov][q1n, q2n], RotateLeft[newfaces[[indlican]], 
            q1n - 1], RotateLeft[newfaces[[indlican]], q2n - 1]], indlican], 
       {m, 1, len - 1}]; newfaces];
PrejQHI[n_][p_, q_] := q == p + 1 || p == n && q == 1;
NetOfConvex3D[solid_, facetree_, ulomek_] := 
    Module[{linijaset, linija, poin1, poin2, d1, d2, len, p1, p2, indn, indp, 
      indlica, novatocka, len1, maksi, dih, afaces, vertis}, 
     afaces = solid[[1]]; vertis = solid[[2]];
      linijaset = Table[{}, {Length[afaces]}];
      normals = Table[{}, {Length[afaces]}];
      maksi = Table[0, {Length[afaces]}]; dih = DihedralD[solid]; 
      indlica = facetree[[1,1]]; maksi = ReplacePart[maksi, 
        Max[afaces[[indlica]]], indlica];  
      normal = ort[Cross[vertis[[afaces[[indlica]][[2]]]] - 
          vertis[[afaces[[indlica]][[1]]]], 
         vertis[[afaces[[indlica]][[3]]]] - vertis[[afaces[[indlica]][[
            1]]]]]]; normals = ReplacePart[normals, normal, indlica]; 
      len = Length[afaces[[indlica]]]; linija = (vertis[[#1]] & ) /@ 
        afaces[[indlica]];  linijaset = ReplacePart[linijaset, linija, 
        indlica]; len1 = Length[facetree]; linija = {};
      Do[indn = facetree[[i,2]]; indp = facetree[[i,1]]; 
        maksi = ReplacePart[maksi, Max[afaces[[indn]]], indn]; 
        {{p1}} = Position[afaces[[indp]], afaces[[indn,1]]]; 
        {{p2}} = Position[afaces[[indp]], afaces[[indn,2]]]; linija = {}; 
        poin1 = linijaset[[indp,p1]]; poin2 = linijaset[[indp,p2]]; 
        linija = Append[linija, poin1]; linija = Append[linija, poin2]; 
        len = Length[afaces[[indn]]]; edgeOrt = ort[poin2 - poin1]; 
        normalp = normals[[indp]]; edgs = EdgesOnFaces[afaces]; 
        {{indedge}} = Position[edgs, Sort[{afaces[[indn,1]], 
            afaces[[indn,2]]}]]; normaln = 
         Rot[edgeOrt, -((Pi - dih[[indedge]])*ulomek)][normalp]; 
        normals = ReplacePart[normals, normaln, indn]; 
        Do[d1 = DistanceHI[vertis[[afaces[[indn,1]]]], 
            vertis[[afaces[[indn,k]]]]]; d2 = DistanceHI[
            vertis[[afaces[[indn,2]]]], vertis[[afaces[[indn,k]]]]]; 
          novatocka = TwoCirclesIntersectionNewHI3D[poin1, poin2, {d1, d2}, 
            normaln]; linija = Append[linija, novatocka], {k, 3, len}]; 
        linijaset = ReplacePart[linijaset, linija, indn], {i, 1, len1}]; 
      {linijaset, maksi}];
DihedralD[telo_] := Module[{c1, solid, f, edgs, v, e, i, pari,koor,normals,rezul}, 
     solid = telo; koor = solid[[2]]; c1 = solid[[1]]; f = Length[c1]; 
      normals = Table[ort[Cross[koor[[c1[[i]][[2]]]] - koor[[c1[[i]][[1]]]], 
          koor[[c1[[i]][[3]]]] - koor[[c1[[i]][[1]]]]]], {i, 1, f}];
      edgs = EdgesOnFaces[c1]; v = Length[Union[Flatten[c1]]];
      e = Length[edgs]; pari = Table[FacesOfEdge[edgs[[i]], c1], {i, 1, e}];
     rezul= Table[N[Pi - ArcCos[normals[[pari[[i]][[1]]]] . 
           normals[[pari[[i]][[2]]]]]], {i, 1, e}];rezul];
ort[x_] := x/Sqrt[x[[1]]^2 + x[[2]]^2 + x[[3]]^2];
EdgesOnFaces[cic_] := Block[{edgs = {}, n = Length[cic], ed, i}, 
     Do[ed = EdgesOnFace[cic[[i]]]; edgs = Union[edgs, ed], {i, 1, n}]; 
      edgs = Union[Map[Sort, edgs, 1]]]
EdgesOnFace[lice_] := Module[{i,rez},rez= Table[{lice[[i]], lice[[i + 1]]}, 
       {i, 1, Length[lice] - 1}];AppendTo[rez, {Last[lice], First[lice]}]];
FacesOfEdge[{a_, b_}, c1_] := Module[{f, v, e}, f = Length[c1]; 
      v = Length[Union[Flatten[c1]]]; e = f + v - 2; 
      Select[Range[f], Length[Intersection[c1[[#1]], {a, b}]] == 2 & ]];
Rot[om_, al_][x_] := x . om*om + (x - x . om*om)*Cos[al] + 
     Cross[om, x]*Sin[al];
DistanceHI[{a_, b_, c_}, {x_, y_, z_}] := 
    N[Sqrt[(a - x)^2 + (b - y)^2 + (c - z)^2]];
TwoCirclesIntersectionNewHI3D[{x1_, y1_, z1_}, {x2_, y2_, z2_}, {r1_, r2_}, 
     n_] := Module[{c, a, b, e}, 
     c = Sqrt[(x2 - x1)^2 + (y2 - y1)^2 + (z2 - z1)^2]; a = r2; b = r1; 
      e = ({x2, y2, z2} - {x1, y1, z1})/c; 
      ConstructTriangleHI3D[{x1, y1, z1}, e, n][a, b, c]];
ConstructTriangleHI3D[{x1_, y1_, z1_}, e_, n_][a_, b_, c_] := 
    Module[{p, hc, cosa, AB, x2, y2, x3, y3, z2, z3}, 
     p = (a + b + c)/2; hc = 2*(Sqrt[p*(p - a)*(p - b)*(p - c)]/c); 
      cosa = (b^2 + c^2 - a^2)/(2*b*c); AB = c*e; 
      {x2, y2, z2} = {x1, y1, y1} + AB; {x3, y3, z3} = 
       {x1, y1, z1} + (b*cosa)*e + Cross[n, e]*hc; {x3, y3, z3}];
ShowMrezaP[p1_,drevo_,ulomek_]:=Block[{tr2,lica3,stran,net},tr2=UrediTree2HI[drevo];stran=p1[[1]];
lica3=ChangeOrientation[tr2][stran];
net=NetOfConvex3D[{lica3,p1[[2]]},tr2,ulomek];Map[Polygon,net[[1]]]];
tree[m_]:=Join[Table[{i,i+1},{i,1,m-1}],{{1,2m-1}},Table[{i,i-1},{i,2m-1,m+2,-1}],{{m-1,2m},{m-1,2m+1}}];
mreza2D[solid_,facetree_]:=Block[{mr,tr2,lica3,stran,net,p1,p2,p3,bv1,bv2,pf},tr2=UrediTree2HI[facetree];stran=solid[[1]];
lica3=ChangeOrientation[tr2][stran];
mr=NetOfConvex3D[{lica3,solid[[2]]},tr2,0];p1=mr[[1,1,1]];p2=mr[[1,1,2]];p3=mr[[1,1,3]];bv1=Normalize[p2-p1];bv2=Normalize[Cross[Cross[p2-p1,p3-p1],p2-p1]];pf[r_]:={(r-p1).bv1,(r-p1).bv2};
Map[pf,mr[[1]],{2}]];
barve=Table[Hue[Random []],{20}];

indexb1[m_]:=Join[Range[m-1],Reverse[Range[m-1]],{m, m+1,m+1}];
ClosedLine[a_]:=Line[Append[a,First[a]]];

Manipulate[
solidpicture=solid[m,ep,n1,k];
With[{kk=(n1-k)/solidpicture[[3]],k1=1,
solid1=solidpicture[[2]],
picture=solidpicture[[1]],
tree1=tree[m],
indexb=indexb1[m]},
With[{
voly=(moments[[dim+1]][solid1])[[2]],
vol=volumes[[dim+1]][solid1],nor=Normalize[Cross[solid1[[2,solid1[[1,2m,2]]]]-solid1[[2,solid1[[1,2m,1]]]],solid1[[2,solid1[[1,2m,3]]]]-solid1[[2,solid1[[1,2m,1]]]]]],
net1=mreza2D[solid1,tree1]},
With[{
vy=voly/vol,
net1A=Table[{barve[[indexb[[i]]]],Polygon[net1[[i]]]},{i,1,Length[net1]}],net1B={GrayLevel[0.0],Thickness[0.001],Map[ClosedLine,net1]}},
Text@Column[{
ExpressionCell[Switch[x1,
2,Show[{picture,Graphics[{Red,PointSize[0.02],Point[{0,vy}]}]},Axes->True,ImageSize->{400,400}],
1,Graphics3D[{Opacity[op],PointSize[0.02],Point[{0,0,0}],{Red,Point[{0,vy,0}]},If[sh1,{{Blue,Thickness[0.003],Line[{{0,0,-n1},{0,n1/kk,0}}]},{Thickness[0.003],Green,Line[{{0,vy,0},{0,vy,0}+14 nor}]}},{}],{Map[Polygon,Map[solid1[[2,#]]&,solid1[[1]]]]}},PlotRange->{{-10,10},{-10,15},zoom{-20,20}},Boxed->False,SphericalRegion->True,Lighting->"Neutral",ImageSize->{400,400},ViewPoint->{10,3,3},ViewVertical->{0,1,0}],3,Graphics[{net1A,net1B},ImageSize->{400,400}]],Deployed->True],Row[{Style["y",Italic]," coordinate of the mass center = ",Style[vy,If[vy≥0,Blue,Red]]}],
Row[{"number of faces = ",2m+1}]}, 
Alignment-> Center]]]],
Row[{"number of faces is 2",Style["m",Italic],"+1"}],
{{m,9,Style["m",Italic]},Range[5,11]},
"",
"angle correction ϵ",
{{ep,0,""},0,0.3,0.05, ImageSize-> Tiny,Setter},
"",
"treat polyhedron as (set of)",
{{dim,3,""},{0->"vertices",1->"edges",2->"faces",3->"solid"}},
{{n1,12,Style["a",Italic]},10,50,0.05,Appearance->"Labeled", ImageSize-> Tiny},
{{k,4,Style["b",Italic]},0.0001,10,0.01,Appearance->"Labeled", ImageSize-> Tiny},
"",
Row[{"show segment of the line ",Style["z",Italic]," = ",Style["h y",Italic]," - ",Style["a",Italic]}],"and perpendicular from mass center","to one prism's base",
{{sh1,True,""},{True,False},Enabled->x1==1},
"",
"opacity",
{{op,0.7,""},0.5,1,Appearance->"Labeled", ImageSize-> Tiny,Enabled->x1==1},
"",


"","show",
Control[{{x1,1,""},{1->"polyhedron",2->"cross-section",3->"2D net"},ControlType->SetterBar}],
"","",
Style["unistable polyhedral…",Bold,10],
Button["… solid with 19 faces",m=9;dim=3;n1=30;k=1.35;ep=0.;zoom=2.5,ImageSize->200],Button["… surface with 15 faces",m=7;dim=2;n1=29;k=1.83;ep=0.2;zoom=2.5,ImageSize->200],
Button["… skeleton with 13 faces",m=6;dim=1;n1=36.35;k=3.1;ep=0.;zoom=2.5,ImageSize->200],
{solidpicture,solid[9,0,12,5],ControlType->None},"",
{zoom,1,3, ImageSize-> Tiny,Enabled->x1==1},
SaveDefinitions->True,SynchronousInitialization->False,
AutorunSequencing-> {1,2,4}]

A polyhedral solid with uniform density is unistable (or monostable) if it is stable on exactly one face. Beside these, this Demonstration also treats polyhedra as surfaces, skeletons, and as set of vertices.

In 1968, Guy constructed a unistable polyhedron with 19 faces. It is a kind of truncated prism. A half of the cross-section of it consists of m=9 right-angled triangles (see cross-section), each of the angles at the coordinate center equals π/m. In our construction we allow the first and the last angle to be a little larger, namely π/m+ϵ (ϵ is called angle correction), and the others to be a bit smaller. If by truncation we succeed that mass center falls below the coordinate origin, only the bottom face (see cross-section) and possibly two truncated bases are stable. But we could make the last two unstable by taking the length 2b of the shortest lateral edge sufficiently large. In this case the perpendicular from mass center (green line) to the coresponding prism's new basis must not intersect it. So the blue line on the basis and green line must meet outside of the solid.

For example, to find a unistable polyhedral surface with 15 faces, choose m=7, "faces", a=29, b=1.83. To get a unistable polyhedron, its center of mass (red point) must be below the origin (black point). The meaning of the parameters a is: z=±a  are z coordinates of intersections of the polyhedron with z axis. So the slope of the blue line is h=(a-b)/r, where r is the distance of the origin to the shortest lateral edge of the polyhedron.

A face of a polyhedron is stable if, and only if, the orthogonal  projection of the center of mass onto the face lies inside the face or on the edge [3, p. 273].

References

[1] J. Bryant and C. Sangwin, How Round Is Your Circle?: Where Engineering and Mathematics Meet, Princeton University Press, 2008 pp. 274–276.

[2] R. K. Guy, A Unistable Polyhedron, University of Calgary Dept. of Mathematics, 1968.

[3] J. H. Conway, M. Goldberg, and R. K. Guy, Problem 66-12, SIAM Review 11 (1969), pp. 78–82.

Contributed by: Izidor Hafner