// disphenoid43adj.fe // Adjoint of fundamental region of triply periodic minimal surface // in disphenoid. Genus 43. // Programmer: Ken Brakke, brakke@susqu.edu, http://www.susqu.edu /* Commands: gogo - typical evolution showcube - display unit cell (not exactly cubic, but same idea) showcubelet - display 1/8 of cubic unit cell showrhombic - show rhombic dodecahedron unit cell transforms off - show just single fundamental region setcolor - to color one side yellow, as in my web page. draw_edges - creates some edges of disphenoid; useful to show rhombic dodecahedron outline. To turn off showing all the edges in the graphics display, hit the "e" key in the graphics window. */ // Follows data from Alan Schoen, December 12, 2002 // parameters found from final TPMS parameter length1 = 0.474779088657332 parameter length2 = 1.042450061923738 parameter length3 = 0.240025282395047 parameter length4 = 1.485105573958310 // Constraints for adjoint edges free on planes constraint 1 formula: z = 0 constraint 6 formula: x = y // Constraints for edges after conjugating, formulated so // right hand sides sides should be positive and equal. parameter rhs2 = 1 parameter rhs3 = 1 parameter rhs4 = 1 parameter rhs5 = 1 constraint 2 formula x - z = rhs2 constraint 3 formula -x - z = rhs3 constraint 4 formula x - z = rhs4 constraint 5 formula y + z = rhs5 constraint 11 formula: x + y = 0 constraint 12 formula: z = 0 constraint 13 formula: x - y = 0 // Transform generators for disphenoid view_transform_generators 5 // C2 about z axis swap_colors -1 0 0 0 0 -1 0 0 0 0 1 0 0 0 0 1 // C2 about (1,-1,0) axis swap_colors 0 -1 0 0 -1 0 0 0 0 0 -1 0 0 0 0 1 // mirror about -x - z = rhs2 0 0 -1 -1 0 1 0 0 -1 0 0 -1 0 0 0 1 // mirror about x - z = rhs3 0 0 1 1 0 1 0 0 1 0 0 -1 0 0 0 1 // mirror about y + z = rhs4 1 0 0 0 0 0 -1 1 0 -1 0 1 0 0 0 1 vertices 1 0.5 0.5 0 constraints 1,6 2 (length3-length1+length2) length4 0 fixed 3 (length3-length1) length4 length2 fixed 4 length3 length4 (length1+length2) fixed 5 0 length4 (length1+length2+length3) fixed 6 0 0 (length1+length2+length3-length4) fixed edges 1 1 2 constraint 1 2 2 3 fixed 3 3 4 fixed 4 4 5 fixed 5 5 6 fixed 6 6 1 constraint 6 faces 1 1 2 3 4 5 6 read read "adjoint.cmd" // adjointing command adj := { unset vertex constraint 1; unset vertex constraint 6; unset edge constraint 1; unset edge constraint 6; unfix vertices; unfix edges; adjoint; } // evolution command gg := { refine edge where valence == 1; g 20; r; g 20; u; V; g 20; r; refine edge where original == 5; g 20; u; V; u; V; g 20; conj_grad; g 20; u; V; u; V; g 20; r; g 20; u; V; g 20; } // to reposition and normalize surface after conjugating frame := { // translate xbase := vertex[1].x; ybase := vertex[1].y; zbase := vertex[1].z; set vertex x x-xbase; set vertex y y-ybase; set vertex z z-zbase; // scale so rhs's all turn out 1 at the end rhs2 := vertex[2].x - vertex[2].z; rhs3 := -vertex[3].x - vertex[3].z; rhs4 := vertex[4].x - vertex[4].z; rhs5 := vertex[5].y + vertex[5].z; maxrhs := rhs2; if rhs3 > maxrhs then maxrhs := rhs3; if rhs4 > maxrhs then maxrhs := rhs4; if rhs5 > maxrhs then maxrhs := rhs5; set vertex x x/maxrhs; set vertex y y/maxrhs; set vertex z z/maxrhs; // set constraints rhs2 := vertex[2].x - vertex[2].z; foreach edge ee where original == 2 do { set ee constraint 2; set ee.vertex constraint 2;}; rhs3 := -vertex[3].x - vertex[3].z; foreach edge ee where original == 3 do { set ee constraint 3; set ee.vertex constraint 3; }; rhs4 := vertex[4].x - vertex[4].z; foreach edge ee where original == 4 do { set ee constraint 4; set ee.vertex constraint 4; }; rhs5 := vertex[5].y + vertex[5].z; foreach edge ee where original == 5 do { set ee constraint 5; set ee.vertex constraint 5; }; foreach edge ee where original == 1 do { set ee.vertex constraint 11; set ee.vertex constraint 13; set ee constraint 11; set ee constraint 13; }; foreach edge ee where original == 6 do { set ee.vertex constraint 11; set ee.vertex constraint 12; set ee constraint 11; set ee constraint 12; }; fix vertices where on_constraint 11; // avoids some spurious neg eigenvalues fix edges where on_constraint 11; } // For equalizing right sides after adjointing. delta := 0.01 kill := { // move lesser rhs's towards biggest by .01 maxrhs := rhs2; if rhs3 > maxrhs then maxrhs := rhs3; if rhs4 > maxrhs then maxrhs := rhs4; if rhs5 > maxrhs then maxrhs := rhs5; if ( maxrhs - rhs2 < delta ) then rhs2 := maxrhs else rhs2 += delta; if ( maxrhs - rhs3 < delta ) then rhs3 := maxrhs else rhs3 += delta; if ( maxrhs - rhs4 < delta ) then rhs4 := maxrhs else rhs4 += delta; if ( maxrhs - rhs5 < delta ) then rhs5 := maxrhs else rhs5 += delta; recalc; } true_lengths := { printf "length1 = %18.15f\n", sum(edge where original==3,length)/sqrt(2); printf "length2 = %18.15f\n", sum(edge where original==2,length)/sqrt(2); printf "length3 = %18.15f\n", sum(edge where original==4,length)/sqrt(2); printf "length4 = %18.15f\n", sum(edge where original==5,length)/sqrt(2); } // Draw disphenoid edges draw_edges := { vert1 := new_vertex(0,0,-rhs2); vert2 := new_vertex(0,2*rhs2,-rhs2); vert3 := new_vertex(rhs2,rhs2,0); vert4 := new_vertex(-rhs2,rhs2,0); newe1 := new_edge(vert1,vert2); set edge[newe1] bare; set edge[newe1] no_refine; newe2 := new_edge(vert2,vert3); set edge[newe2] bare; set edge[newe2] no_refine; newe3 := new_edge(vert2,vert4); set edge[newe3] bare; set edge[newe3] no_refine; } length1min := .80; length1max := 1.2; length2min := .8; length2max := 1.2; length3min := .8; length3max := 1.2; dlength := 0.05; goflag := 1; tester := { read "params.txt"; vertex[2].x := length3 + length2 - length1; vertex[2].y := length4; vertex[3].x := length3-length1; vertex[3].y := length4; vertex[3].z := length2; vertex[4].x := length3; vertex[4].y := length4; vertex[4].z := length1+length2; vertex[5].y := length4; vertex[5].z := length1 + length2 + length3; vertex[6].z := length1 + length2 + length3 - length4; gg; adj; frame; diff := abs(rhs2-rhs3) + abs(rhs2-rhs4) + abs(rhs3-rhs4) + abs(rhs5-rhs2) + abs(rhs5 - rhs3) + abs(rhs5-rhs4); printf "diff %8.6f at %f %f %f\n",diff,length1,length2,length3 >> "43.out"; if ( length1 < length1max ) then length1 += dlength else if ( length2 < length2max ) then { length1 := length1min; length2 += dlength; } else if ( length3 < length3max ) then { length1 := length1min; length2 := length2min; length3 += dlength; } else { goflag := 0; return;}; printf "length1 := %f; length2 := %f; length3 := %f;\n", length1,length2,length3 >>> "params.txt"; } //tester //if goflag then load "disphenoid43adj"; // Handy transformations of TMPS showcubelet := { transform_expr "cecb"; show_trans "R"; } showcube := { transform_expr "cdecdeb"; show_trans "R"; } showrhombic := { transform_expr "cdecdeab"; show edge where valence == 1 or ( bare and id != newe1); show_trans "R"; } setcolor := { set facet frontcolor yellow } gogo := { gg; adj; frame; kill; u;V; hessian; hessian; show_trans "R"; } /* Commands: gogo - typical evolution showcube - display unit cell (not exactly cubic, but same idea) showcubelet - display 1/8 of cubic unit cell showrhombic - show rhombic dodecahedron unit cell transforms off - show just single fundamental region setcolor - to color one side yellow, as in my web page. draw_edges - creates some edges of disphenoid; useful to show rhombic dodecahedron outline. To turn off showing all the edges in the graphics display, hit the "e" key in the graphics window. */