// disphenoid67adj.fe

// Adjoint of fundamental region of triply periodic minimal surface
// in disphenoid.  Genus 67

// 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 13, 2002

// parameters found from final TPMS
parameter length1 =  0.405451027240078
parameter length2 =  0.454486867875058
parameter length3 =  0.731395388926525
parameter length4 =  0.545918212617905
parameter length5 =  1.366706447193993

// Constraints for adjoint edges free on planes
constraint 1
formula: z = 0

constraint 7
formula: x - y = 0

// 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
parameter rhs6 = 1
constraint 2 formula -y - z = rhs2
constraint 3 formula -x + z = rhs3
constraint 4 formula y - z = rhs4
constraint 5 formula -x + z = rhs5
constraint 6 formula y - z = rhs6

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

// a: C2 about z axis
swap_colors
-1 0 0 0 
 0 -1 0 0
 0  0 1 0
 0  0 0 1

// b; C2 about (1,-1,0) axis
swap_colors
0 -1 0 0
-1 0 0 0
 0 0 -1 0
 0 0 0 1

// c: mirror about -y - z = rhs2
1 0 0 0 
0 0 -1 -1
0 -1 0 -1
0 0 0 1
 
// d: mirror about  -x + z = rhs3
0 0 1 -1
0 1 0 0
1 0 0 1
0 0 0 1
 
// e: mirror about  y - z = rhs4
1 0 0 0
0 0 1 1
0 1 0 -1
0 0 0 1
 
vertices
1  (length1+length2+length3+length4-length5) \
       (length1+length2+length3+length4-length5) 0 constraints 1,7
2  0 (length1+length2+length3+length4-length5) 0 fixed
3  0 (length4+length3+length2+length1) length5  fixed
4  length1 (length4+length3+length2+length1) (length5-length1) fixed
5  length1 (length1+length3+length4) (length5-length1+length2) fixed
6  (length1+length3) (length1+length3+length4) \
          (length5-length1+length2-length3) fixed
7  (length1+length3) (length1+length3)  \
         (length5+length4+length2-length1-length3) 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 7 fixed
7  7 1 constraint 7

faces
1   1 2 3 4 5 6 7

read

read "adjoint.cmd"

// adjointing command
adj := { unset vertex constraint 1; unset vertex constraint 7; 
         unset edge constraint 1; unset edge constraint 7;
         unfix vertices; unfix edges;
         adjoint; 
}

// evolution command
gg := { refine edge where valence == 1;
          g 20; r; g 20; u; V; g 20;
          r; 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;
          refine edge where fixed and valence == 1; 
          u; V; u; V; 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].y - vertex[2].z;
  rhs3 := -vertex[3].x + vertex[3].z;
  rhs4 := vertex[4].y - vertex[4].z;
  rhs5 := -vertex[5].x + vertex[5].z;
  rhs6 := vertex[6].y - vertex[6].z;
  maxrhs := rhs2;
  if rhs3 > maxrhs then maxrhs := rhs3;
  if rhs4 > maxrhs then maxrhs := rhs4;
  if rhs5 > maxrhs then maxrhs := rhs5;
  if rhs6 > maxrhs then maxrhs := rhs6;
  set vertex x x/maxrhs;
  set vertex y y/maxrhs;
  set vertex z z/maxrhs;

  // set constraints
  rhs2 := -vertex[2].y - 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].y - vertex[4].z;
  foreach edge ee where original == 4 do 
  { set ee constraint 4;  set ee.vertex constraint 4; };
  rhs5 := -vertex[5].x + vertex[5].z;
  foreach edge ee where original == 5 do 
  { set ee constraint 5; set ee.vertex constraint 5; };
  rhs6 := vertex[6].y - vertex[6].z;
  foreach edge ee where original == 6 do 
  { set ee constraint 6;  set ee.vertex constraint 6; };

  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 == 7 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 rhs6 > maxrhs then maxrhs := rhs6;

  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;
  if ( maxrhs - rhs6 < delta ) then rhs6 := maxrhs
    else rhs6 += 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==4,length)/sqrt(2);
                  printf "length3 = %18.15f\n", 
                    sum(edge where original==5,length)/sqrt(2);
                  printf "length4 = %18.15f\n", 
                    sum(edge where original==6,length)/sqrt(2);
                  printf "length5 = %18.15f\n", 
                    sum(edge where original==2,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;
}

// 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 backcolor yellow }

gogo := { gg; adj; frame; kill; 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.
*/