// hybrid-1.fe
// Adjoint of Schoen's hybrid P, F-RD 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.
*/

parameter asize = 0.44790391042324   // shape parameter, for period killing


// Constraints for use after adjoint transformation
constraint 3 formula: z = 1
constraint 4 formula: x + y = 0
constraint 5 formula: x - y = 0 
constraint 6 formula: y - z = 0

view_transform_generators 6
 0 1 0 0    1 0 0 0   0 0 1 0   0 0 0 1  // a: x-y swap
 0 -1 0 0  -1 0 0 0   0 0 1 0   0 0 0 1  // b: x+y swap
 1 0 0 0    0 0 1 0   0 1 0 0   0 0 0 1  // c: y-z swap
-1 0 0 2    0 1 0 0   0 0 1 0   0 0 0 1  // d: x=1 mirror
 1 0 0 0    0 -1 0 2  0 0 1 0   0 0 0 1  // e: y=1 mirror
 1 0 0 0    0 1 0 0   0 0 -1 2  0 0 0 1  // f: z=1 mirror
  
vertices
1   0 0 0 fixed
2   0 0 asize fixed
3   -.5 -.5 asize fixed
4   -.5 -.5 1 fixed
5    0  -1 1 fixed

edges
1  1 2 fixed
2  2 3 fixed
3  3 4 fixed
4  4 5 fixed
5  5 1 fixed

faces
1   -5 -4 -3 -2 -1

read
hessian_normal

// good evolution, getting lots of facets near vertex 2 cusp.
gg := { refine edge where valence == 1; g 5; r; g 10; u; V;
          g5; hessian;hessian; 
          g 5; hessian; hessian;
          r; g 5; u; V; u; g 5; hessian; hessian;
          g5; hessian;hessian; 
          r; g 5; u; V; u; g 5; hessian; hessian;
          g 5; V; u; V; hessian; hessian;
          r; refine edge where original==5;
          g 5; u; V; u; g 5; hessian; hessian;
          g 5; V; u; V; hessian; hessian;
        }

// Some distances in the adjoint
calc := {
          edge2dz := sum(edge ee where original==2,
           sum(ee.facet ff, (ff.x*ee.y-ff.y*ee.x)/sqrt(ff.x^2+ff.y^2+ff.z^2)));

          printf " edge2dz %g   \n",edge2dz;
        }


read "adjoint.cmd"

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

// Applying constraints after adjointing
frame := {
    aa := max(vertex,x-y);
    bb := min(vertex,x+y);
    set vertex x  x-(aa+bb)/2;
    set vertex y  y-(bb-aa)/2;
    dd := vertex[5].y - vertex[5].z;
    set vertex z  z+dd;
    mag := max(vertex,z);
    set vertex x x/mag;
    set vertex y y/mag;
    set vertex z z/mag;
    foreach edge ee where original==1 do
    { set ee.vertex constraint 3; set ee constraint 3; };
    foreach edge ee where original==2 do
    { set ee constraint 4; set ee.vertex constraint 4; };
    foreach edge ee where original==3 do
    { set ee constraint 3; set ee.vertex constraint 3; };
    foreach edge ee where original==4 do
    { set ee constraint 5; set ee.vertex constraint 5; };
    foreach edge ee where original==5 do
    { set ee constraint 6; set ee.vertex constraint 6; };
}


// To get true asize after evolving after adjointing
true_asize := { printf "True asize: %20.15f\n",
   sum(edge where original == 1,length)/sum(edge where original==2,length)*sqrt(.5); }

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

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

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

// Show unit cell
showcube := { transform_expr "defabcabc"; 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.
*/