// schoen12adj.fe
// Adjoint surface of Schoen's p. 12 surface

// Programmer: Ken Brakke, brakke@susqu.edu, http://www.susqu.edu

/* Commands:
   gogo - typical evolution
   showcube - show one unit cell
   showcubelet - show 1/8 of unit cell, as on web page
   transforms off - show just single fundamental region
   setcolor - to color one side yellow, as in my web page.

   To turn off showing all the edges in the graphics display,
      hit the "e" key in the graphics window.
*/

// Adjoint frame edge lengths
parameter asize = 1
parameter bsize = 2
parameter dsize = 2.778

// Constraints for use after adjointing
parameter minx = 0
parameter maxx = 0
parameter miny = 0
parameter maxy = 0
parameter minz = 0
parameter maxz = 0
constraint 1 formula: x = minx
constraint 2 formula: x = maxx
constraint 3 formula: y = miny
constraint 4 formula: y = maxy
constraint 5 formula: z = minz
constraint 6 formula: z = maxz

view_transform_generators 6 
-1 0 0 2*minx  0 1 0 0         0 0 1 0        0 0 0 1  // a: mirror in x=minx
-1 0 0 2*maxx  0 1 0 0         0 0 1 0        0 0 0 1  // b; mirror in x=maxx
 1 0 0 0       0 -1 0 2*miny   0 0 1 0        0 0 0 1  // c: mirror in y=minx
 1 0 0 0       0 -1 0 2*maxy   0 0 1 0        0 0 0 1  // d: mirror in y=maxx
 1 0 0 0       0 1 0 0         0 0 -1 2*minz  0 0 0 1  // e: mirror in z=minx
 1 0 0 0       0 1 0 0         0 0 -1 2*maxz  0 0 0 1  // f: mirror in z=maxx
vertices
1  -dsize 0 0 fixed
2  asize 0 0  fixed
3  asize -dsize 0 fixed
4  0 -dsize 0 fixed
5  0 -dsize bsize fixed
6  0 asize bsize fixed
7  -dsize asize bsize fixed
8  -dsize 0 bsize fixed

edges
1  1 2 fixed
2  2 3 fixed
3  3 4 fixed
4  4 5 fixed
5  5 6 fixed
6  6 7 fixed
7  7 8 fixed
8  8 1 fixed

faces
1    -8 -7 -6 -5 -4 -3 -2 -1


read
hessian_normal

// calculate gap
calc := { print sum(edge ee where original==1, 
           sum(ee.facet ff, ff.z/sqrt(ff.x^2+ff.y^2+ff.z^2)*ee.x));
          print sum(edge ee where original==8, 
           sum(ee.facet ff, ff.x/sqrt(ff.x^2+ff.y^2+ff.z^2)*ee.z))
        }

// For changing parameters
new_a := asize;
new_d := dsize;
do_a := { set vertex y y*new_a/asize where y > 0;
          set vertex x x*new_a/asize where x > 0;
          asize := new_a;
        }
do_d := { set vertex y y*new_d/dsize where y < 0;
          set vertex x x*new_d/dsize where x < 0;
          dsize := new_d;
        }

// initial evolution
gg := { r; g 5; refine edge where valence == 1;
          r; g 12; u; V; g 12; u; V; g 12;
          r; g 12; u; V; g 12; u; V; g 12;
          r; g 12; u; V; g 12; u; V; g 12;
  }
 
// with less refinement
gg1 := { r; g 5; refine edge where valence == 1;
          r; g 12; u; V; g 12; u; V; g 12;
}

// color the edges
trim := { foreach edge ee where original >= 1 and original <= 8 do
            set ee.facet color ee.original
        }
 
// A little exploration of parameter space
explore := {
  while ( asize > .4 ) do
  {
    new_d := 1; do_d;
    while ( dsize > .4 ) do
    { 
      g 5; u; V; g 12; u; V; g 20; hessian;  hessian;
      printf "a_size %g  dsize %g \n",asize,dsize;
      calc;
      new_d := dsize - 0.1; do_d;
    };
    new_a := asize - 0.1;
    do_a;
  }
}


//  Shrinking unit cell down to cube?  can't get there.
shrink := { 
   set vertex z  minz+(z-minz)*.99;
   maxz := minz + (maxz-minz)*.99;
   hessian; hessian;
}

read "adjoint.cmd"

// Call this to do adjoint transformation!
adj := { 
         adjoint;
       }


// Applying box constraints after adjointing
frame := { 
    minx := min(vertex,x);
    maxx := max(vertex,x);
    miny := min(vertex,y);
    maxy := max(vertex,y);
    minz := min(vertex,z);
    maxz := max(vertex,z);
    foreach edge ee where original == 1 do
     { set ee constraint 2; set ee.vertex constraint 2; };
    foreach edge ee where original == 2 do
     { set ee constraint 3; set ee.vertex constraint 3; };
    foreach edge ee where original == 3 do
     { set ee constraint 1; set ee.vertex constraint 1; };
    foreach edge ee where original == 4 do
     { set ee constraint 6; set ee.vertex constraint 6; };
    foreach edge ee where original == 5 do
     { set ee constraint 4; set ee.vertex constraint 4; };
    foreach edge ee where original == 6 do
     { set ee constraint 1; set ee.vertex constraint 1; };
    foreach edge ee where original == 7 do
     { set ee constraint 3; set ee.vertex constraint 3; };
    foreach edge ee where original == 8 do
     { set ee constraint 5; set ee.vertex constraint 5; };
    unfix vertices;
    unfix edges;
}




// Set inside color to yellow
setcolor := { set facet backcolor yellow }

// Typical evolution 
gogo := { gg; adj; frame; show_trans "R"; V; V; hessian; hessian; }

//  Evolution with less refinement
gogo1 := { gg1; adj; frame; show_trans "R"; V; V; hessian; hessian; }

// Show 1/8 unit cell
showcubelet := { transform_expr "ca"; show_trans "R"; }

// Show unit cell
showcube := { transform_expr "fca"; show_trans "R"; }


/* Commands:
   gogo - typical evolution
   showcube - show one unit cell
   showcubelet - show 1/8 of unit cell, as on web page
   transforms off - show just single fundamental region
   setcolor - to color one side yellow, as in my web page.

   To turn off showing all the edges in the graphics display,
      hit the "e" key in the graphics window.
*/