Message Boards Message Boards

A tetrahedral chain challenge

Stan Wagon, Macalester College
Posted 12 years ago


In 1958 it was proved that a chain of regular tetrahedra, meeting face to face, cannot form a perfect loop: the last tetrahedron cannot equal the first, exactly. I found a loop that comes close. Picture here:



This chain of 48 tetrahedra has a gap, but we can introduce some small error at each step so as to give the illusion of this being a perfect ring. My attempt is below, and that has about 1% error in each tetrahedron. The challenge is to reduce this error to 0.5% or less so as to create a better illusion of regularity. Here is what it looks like.



So the challenge is to reduce the error.  Ask me for the code or more details.
POSTED BY: Stan Wagon
9 Replies
Sort By:
That looks like a nice challenge! But far from trivial! Does it need to be a 'loop' or can it be any shape?
POSTED BY: Sander Huisman
Stan, this sounds very exciting. Would you mind providing your original code so people have a starting point? If you have trouble posting it yourself - just email it to me and I will post it.
POSTED BY: Vitaliy Kaurov
I will get the code up eventually.

We certainly want to avoid self-intersecting tetrahedra here. So the idea is to have any sort of a loop -- spiral, knot, etc -- that leaves a gap that is small enough that it can be absorbed by changing some of the tests by 1% or less. I will have a 3D printer model of my first try soon, but I want to work a bit now on getting the error under 1%. I feel that the human eye might well be able to sense a 1% error, tho maybe in the context of the loop it will be hard to sense. I will see when I get my model. It seems pretty clear that by going up from 48 tetrahedra one can get a smaller gap. Gatting a small gap with fewer tetrahedra would be interesting, but for the "hoax" it is the "average gap" that matters I suppose.

Vitaly: I thought I was "following" this thread, but I got no email saying posts are up.

Edit: Working with Michael Elgersma, we now have error of around 0.25%. That should be well good enough for a model where each tetrahedron looks regular to our eye. Of course, there is still the chance of reducing the error, or finding a completely different sequence of non-intersecting tetrahedra.

Edit: Working hard here. Confirmed that the error can be as low as 0.25%. Here is a picture of a model from Ken Sloan's lab.
POSTED BY: Stan Wagon
Posting detailed code here. The initial image, with the gap: I zeroed out the "bestVec", which was used in an earlier attempt to optimize
  faces = {{2, 3, 4}, Reverse@ {1, 3, 4}, {1, 2, 4}, {1, 2, 3}};
  faces = {{2, 3, 4}, Reverse@ {1, 3, 4}, {1, 2, 4},
     Reverse@{1, 2, 3}}; (* important *)
  facesRev = Reverse /@ faces;
  startV = N@{{0, 0,
       Sqrt[2/3] - 1/(2 Sqrt[6])}, {-(1/(2 Sqrt[3])), -(1/2), -(1/(
        2 Sqrt[6]))}, {-(1/(2 Sqrt[3])), 1/2, -(1/(2 Sqrt[6]))}, {1/
       Sqrt[3], 0, -(1/(2 Sqrt[6]))}};
  
 
 bestVec =
   0 Flatten[{{0.0030786234640517588`, -0.0012713757639448969`,
       0.005921062251657448`}, {-0.0017401845032882991`, \
 -0.0020681936130036126`,
       0.005393869601475211`}, {0.0020298099367479444`,
       0.0005102453293635576`,
       0.0029724999287053682`}, {-0.00023558451188598613`,
       0.0009452979460132543`,
       0.005798188590794252`}, {-0.000027554327839240944`,
       0.0020459864295050674`,
       0.0034435421203985097`}, {0.0025981484730485446`,
       0.0007355749150532`,
       0.0012539109619916077`}, {-0.006626782180945919`, \
 -0.0005807683214184328`, 0.0038058021075369495`}}];
 
 
 choices = {1, 2, 3, 3, 2, 1, 1, 2, 3, 3, 2, 1, 1, 2, 3, 3, 2, 1, 1, 2,
     3, 3, 2, 1, 1, 2, 3, 3, 2, 1, 1, 2, 3, 3, 2, 1, 1, 2, 3, 3, 2, 1,
    1, 2, 3, 3, 2,
    1}; (* found by luck and trial and error in August 2013  *)
 
 
 perms = Permutations[Range[4]];
 dist[tet1_, tet2_] :=
   Min[Table[Max[Norm /@ ((tet1 - tet2[[pp]]))], {pp, perms}]];
 
 
 nextTetPushed[currentV_, freeVert_, moveVertInd_, pVec_] := (
   vertToMove = DeleteCases[Range[4], freeVert][[moveVertInd]];
   vertToMovVal = currentV[[vertToMove]];
   reflectface = DeleteCases[Range[4], vertToMove];
   centroid = Mean[currentV[[reflectface]]];
   newpos = vertToMovVal + 2 (centroid - vertToMovVal);
   newVerts = ReplacePart[currentV, vertToMove -> newpos];
   {# + pVec & /@ newVerts, vertToMove})
 
 
 count = 0; Clear[func, st];
 func[{a1_?NumericQ, a2_, a3_, a4_, a5_, a6_, a7_, a8_, a9_, a10_,
     a11_, a12_, a13_, a14_, a15_, a16_, a17_, a18_, a19_, a20_,
     a21_}] :=  (st[0] = {startV, 1};
    change =
     Flatten[Table[{{a1, a2, a3}, {a4, a5, a6}, {a7, a8, a9}, {a10,
         a11, a12}, {a13, a14, a15}, {a16, a17, a18}, {a19, a20,
         a21}}, {10}], 1];
    
    Do[st[i] =
      nextTetPushed[st[i - 1][[1]], st[i - 1][[2]],  choices[[i]],
       change[[i]]], {i, Length[choices]}];
    L = Length[choices];
    data = st /@ Range[0, L];
    store = {{a1, a2, a3}, {a4, a5, a6}, {a7, a8, a9}, {a10, a11,
       a12}, {a13, a14, a15}, {a16, a17, a18}, {a19, a20, a21}};
    bb = dist[data[[1, 1]], data[[-1, 1]]]; count++;
    If[Mod[count, 500] == 0, Print[{bb, Norm[Flatten[store]]}]]; bb);
 
 
 funcPercent[{a1_?NumericQ, a2_, a3_, a4_, a5_, a6_, a7_, a8_, a9_,
     a10_, a11_, a12_, a13_, a14_, a15_, a16_, a17_, a18_, a19_, a20_,
     a21_}] :=  (st[0] = {startV, 1};
    change =
     Flatten[Table[{{a1, a2, a3}, {a4, a5, a6}, {a7, a8, a9}, {a10,
         a11, a12}, {a13, a14, a15}, {a16, a17, a18}, {a19, a20,
         a21}}, {10}], 1];
    
    Do[st[i] =
      nextTetPushed[st[i - 1][[1]], st[i - 1][[2]],  choices[[i]],
       change[[i]]], {i, Length[choices]}];
    L = Length[choices];
    data = st /@ Range[0, L];
    store = {{a1, a2, a3}, {a4, a5, a6}, {a7, a8, a9}, {a10, a11,
       a12}, {a13, a14, a15}, {a16, a17, a18}, {a19, a20, a21}};
    bb = dist[data[[1, 1]], data[[-1, 1]]]; count++;
    dataAll = Flatten[First /@ data, 1];
    dubs =
     Table[ Count[dataAll,
       p_ /;  Norm[p - dataAll[[i]]] < .35 /1000]  , {i,
       Length[dataAll]}];
    dups = dataAll[[Position[dubs, 2] // Flatten]];
    close[p_] := Select[dataAll, Norm[# - p] < .35 /1000 &];
    close[p_] := False;
    dataClose = close /@ dataAll;
    dataAllFixed = Table[ Mean[dataClose[[i]]], {i, Length[dataAll]}];
    dataTetFixed = Partition[dataAllFixed, 4];
    ans = {Max[#], Min[#]} & @
      Flatten[Table[
        Norm[#[[1]] - #[[2]]] & /@ Subsets[dataTetFixed[[i]], {2}], {i,
          48}]];
    ans = ans[[1]] - ans[[2]];
   If[ans < best, best = ans; storedBest = store];
    If[Mod[count, 200] == 0,
    Print[{count, ans, bb, Norm[Flatten[store]], best}]]; ans + 10 bb);


funcPercent[Flatten[bestVec]]; (* important *)
dataAll = First /@ data;
aa = {Max[#], Min[#]} & @
   Flatten[Table[
     Norm[#[[1]] - #[[2]]] & /@ Subsets[dataAll[[i]], {2}], {i,
      48}]];
aa[[1]] - aa[[2]];
dataPols = First /@ data;
{vp, vv} = {{-0.5736193192422712`, -1.2442823717517009`,
    0.40843880317618597`}, {-0.6180067894290332`,
    0.8698216744526492`, -0.3761742628418529`}};
cols = Lighter /@
   Prepend[Join @@ Table[{Blue, Green, Green, Blue, Red, Red}, {50}],
    Red];
gg = Graphics3D[
  Table[{ EdgeForm[Black], FaceForm[{Opacity[.9], cols[[i]]}],
    FaceForm[Lighter@Red, Green],
    Polygon[(dataPols)[[i]][[#]]] & /@
     If[EvenQ[i], facesRev, faces]}, {i, Length@data - 1}],
  BoxRatios -> Automatic, ViewAngle -> .47 2, ViewPoint -> vp,
  ViewVertical -> vv, PlotRegion -> {{-.15, 1.15}, {-.25, 1.45}},
  PlotRange -> {All, {{-1, 1}, {-1, 1}, {-1, 1}}}[[1]],
  Boxed -> False, Axes -> ! True, ImageSize -> 700,
  Lighting -> {Automatic, "Neutral"}[[2]]]



There are some orientation issues.... but they were fixed above. Inside is green. Outside is red. So orientation is all good. Next it makes sense to remove the interior faces.

 polys = Cases[gg, Polygon[_], \[Infinity]];
 Length[polys];
 polys1 = polys /. Polygon[z_] :> Polygon[Sort[z]];
 Last /@ (tt = Tally[polys]);
 (polys1 = First /@ Select[tt, #[[2]] >= 1 &]) // Length
 g1 = Graphics3D[{Opacity[1], FaceForm[Lighter@Red, Green],
    Table[Polygon[polys1[[i, 1]]], {i, Length[polys1]}]},
   Boxed -> False, Axes -> ! True, AxesLabel -> {x, y, z},
   ImageSize -> 600, Lighting -> "Neutral"]
polys = Cases[g1, Polygon[_], \[Infinity]];


Moving the error around to make the last = first. Joint work with Michael Elgersma.This defines the rim edges.
special[{L1_, L2_}] := ((( L1[[1]] < .4 - .3  L1[[ 2]]) &&
       ( L2[[1]] < .4 - .3  L2[[ 2]])) ||
     (((
         L1[[1]] > .4 - .3  L1[[ 2]]) &&
        ( L2[[1]] > .4 - .3 L2[[ 2]])))) && (Norm[
       L1[[{2, 3}]] - cen[[{2, 3}]]] > 1.7 &&
    
     Norm[L2[[{2, 3}]] - cen[[{2, 3}]]] > 1.7 ) ;

This shows the rim edges in yellow. Sort of like railway tracks.
 cen = {.7, 2, 0};
 qq = Flatten[
    Cases[g1,
     Polygon[z_] :> {Line[{z[[1]], z[[2]]}], Line[{z[[2]], z[[3]]}],
       Line[{z[[3]], z[[1]]}]}, \[Infinity]], 1];
 Show[g1, Graphics3D[{Thickness[.04],
    Table[{www =
       Which[special@q, Directive[{Yellow, Thickness[.02]}] ,
        special@q, Directive[{Yellow, Thickness[.02]}] ,
       True, Directive[{Black, Thickness[.004]}]],
     Line@q}, {q, First /@ qq}]}](*, plane;*) ,
PlotRange -> {{-2.2, 2.2}, {-1, 5}, {-2.5, 2.5} 1.1},
ImageSize -> 500, Boxed -> False, Axes -> False]




Now the idea is to start the process from the beginning with a parameter, \, that is used to length the rim edges, while all the other edges remain at length 1. This is a little trickier than I thought. I use SOLVE at each point. Here we change the start to allow \ for the one rim edge.
 startVerts[\[Epsilon]_] := {startV [[1]], {x, y, z} /. Quiet@Solve[{
         ({x, y, z} - startV[[1]]).({x, y, z} - startV[[1]]) == (
            1 + \[Epsilon])^2,
         ({x, y, z} - startV[[3]]).({x, y, z} - startV[[3]]) == 1,
         ({x, y, z} - startV[[4]]).({x, y, z} - startV[[4]])  ==
          1}, {x, y, z}][[1]], startV[[3]], startV[[4]]};
 
 startVerts[.01 0]
 Out[] = {{0., 0., 0.612372}, {-0.288675, -0.5, -0.204124}, {-0.288675, 0.5, -0.204124}, {0.57735, 0., -0.204124}}

 nextTetElgersma[currentV_, freeVert_, moveVertInd_, \[Epsilon]_] := (
    vertToMove = DeleteCases[Range[4], freeVert][[moveVertInd]];
    vertToMovVal = currentV[[vertToMove]];
    reflectface = DeleteCases[Range[4], vertToMove];
    centroid = Mean[currentV[[reflectface]]];
    newpos = vertToMovVal + 2 (centroid - vertToMovVal);
    oldVerts = Delete[currentV, vertToMove];
    edges = {#, newpos} & /@ oldVerts;
    pos = Position[edges, _?special];
   newVerts =  ReplacePart[currentV, vertToMove ->
      If[pos == {},
       Select[{x, y, z} /. Quiet@Solve[{
             ({x, y, z} - edges[[1, 1]]).({x, y, z} -
                 edges[[1, 1]]) == 1,
             ({x, y, z} - edges[[2, 1]]).({x, y, z} -
                 edges[[2, 1]]) == 1,
             ({x, y, z} - edges[[3, 1]]).({x, y, z} -
                 edges[[3, 1]]) == 1}, {x, y, z}],
         Norm[# - vertToMovVal] > .5 &][[1]],
      
       ss = edges[[pos[[1, 1]]]];
       others = Complement[Range[3], {pos[[1, 1]]}]; specCount++;
       Select[{x, y, z} /. Quiet@Solve[{
             ({x, y, z} - ss[[1]]).({x, y, z} - ss[[1]]) == (
                1 + \[Epsilon])^2,
             ({x, y, z} - edges[[others[[1]], 1]]).({x, y, z} -
                 edges[[others[[1]], 1]]) == 1,
             ({x, y, z} - edges[[others[[2]], 1]]).({x, y, z} -
                 edges[[others[[2]], 1]]) == 1}, {x, y, z}],
         Norm[# - vertToMovVal] > .5 &][[1]]  ]];
   {newVerts, vertToMove});

The following was found using FindMinimum and also some trial and error at the start. The value of 0.0027 means that each edge is between 1. and 1.0027 in length. So the error in terms of difference from a regular tetrahedron of side length 1.00137 is about 0.13%, or 1/7 of one percent. This is likely within the tolerances of a 3D printing machine. Still, it would be nice to make it even smaller.
bestepsilonvalue = 0.002743438209350663;

The last tetrahedron is now indistinguishable from the first. 
funcElgersma[.002743438209350663]
Out[] = 2.06186*10^-14

We fix the last points to match exactly, though at 10^-14 difference this hardly matters.
data[[-2]]
data[[1]]
data[[-2]] = {{data[[-2, 1, 1]], data[[1, 1, 2]], data[[1, 1, 3]], data[[1, 1, 4]]}, data[[1, 2]]}

Here is a picture of the fake tetrahedral torus.
 gg = Graphics3D[
   Table[{ EdgeForm[Black],
     FaceForm[{Opacity[1], GrayLevel[RandomReal[{.3, .9}]];
       cols[[i]]}], Polygon[(First /@ data)[[i]][[#]]] & /@ faces}, {i,
      Length@data - 1}], BoxRatios -> Automatic, ViewAngle -> .47 2,
   ViewPoint -> vp, ViewVertical -> vv,
   PlotRegion -> {{-.15, 1.15}, {-.25, 1.45}},
   PlotRange -> {All, {{-1, 1}, {-1, 1}, {-1, 1}}}[[1]],
   Boxed -> False, Axes -> ! True, ImageSize -> 500,
  Lighting -> {Automatic, "Neutral"}[[2]]]
POSTED BY: Stan Wagon
Stan thank you for sending the code. These are some very nice developments. The 3D printed image looks great! 
POSTED BY: Vitaliy Kaurov
We now have much much improved results, including one example where the error is about the diameter of a proton!

Contact me for details or images.

Stan Wagon
wagon@macalester.edu
POSTED BY: Stan Wagon
Interesting problem.  A practical case is twinning in face centered cubic crystals, where multiple twinning can produce what appear to be 5-fold symmetric facets by reflections through 111 planes.  At the 1990 Materials Research Society meeting Linus Pauling invoked this process to refute evidence of quasicrystals (he was wrong).

Can you post the citation for the 1958 proof?  I found one in the Mathematical Gazette vol 56 no 397 pp194-197 (1978).  
POSTED BY: Fred Doty
The first proof of nonexistence is:
S. ?wierczkowski, On a free group of rotations of the Euclidean space, Indag. Math. 20 (1958) 376\378.
Problem first posed:  H. Steinhaus, Problem 175, Coll. Math., 4 (1957) 243.
I now have a demo on the web site allowing some exploring, and showcasing some of our nicer loops. 

One does wonder whether these loops have physical significance. 

We continue to work hard and are continually finding new good stuff. 

Send me private email if you want to be on mailing list for paper.

Stan 
POSTED BY: Stan Wagon
New record. 148 tetrahedra in a loop such that the gap at the end has error 7 . 10^-14. So instead of the diameter of a proton, we (Mike Elgersman and I) are now down to the radius of a proton. Our previous record was 174 tetrahedra with a gap of about 1.4 10^-13.

Stan Wagon
POSTED BY: Stan Wagon
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract