// triplane5adj.fe
// Adjoint of 5nd surface in triplane series.
//     5 crossings of z axis edge in fundamental tetrahedron
//     Five periods to kill.

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

/* Commands:
   gogo - typical evolution; do before the following commands.
   showcube - display cubic unit cell
   showrhombic - show rhombic dodecahedron unit cell
   rhombic_edges - create edges outlining rhombus
   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 = 1.000105537718361 // for period killing;
parameter bsize = 0.999593653771789 // for period killing;
parameter csize = 1.004876615471779 // for period killing;
parameter dsize = .931036630870728 // for period killing; 
parameter esize = .425512594794573 // for period killing; 


constraint 1 // mirror plane in adjoint
formula: x = 0 

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

view_transform_generators 7
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=0 mirror
0,0,1,0    0,1,0,0   1,0,0,0   0,0,0,1  // c: x,z swap
swap_colors 1 0 0 2 , 0  1 0 0, 0 0  1 0, 0 0 0 1 // d: x translation
swap_colors 1 0 0 0 , 0  1 0 2, 0 0  1 0, 0 0 0 1 // e: y translation
swap_colors 1 0 0 0 , 0  1 0 0, 0 0  1 2, 0 0 0 1 // f: z translation
swap_colors 1 0 0 0,  0 -1 0 0, 0 0 -1 2, 0 0 0 1 // g: C2 rotation

vertices
1   0 0 0 fixed
2   -1 -1 0 fixed
3   -1-asize -1+asize 0 fixed
4   -1-asize-bsize  -1+asize-bsize  0 fixed
5   -1-asize-bsize-csize  -1+asize-bsize+csize  0 fixed
6   -1-asize-bsize-csize-dsize  -1+asize-bsize+csize-dsize  0 fixed
7   -1-asize-bsize-csize-dsize-esize  -1+asize-bsize+csize-dsize+esize  0 fixed
8   0 -1+asize-bsize+csize-dsize+esize -1-asize-bsize-csize-dsize-esize 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 constraint 1

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

read
hessian_normal

// good evolution
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; g 5; u; V; u; g 5; hessian; hessian;
          g 5; V; u; V; hessian; hessian;
        }

// Some distances in the adjoint
calc := {
          edge1dy := sum(edge ee where original==1,
           sum(ee.facet ff, (ff.z*ee.x-ff.x*ee.z)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          printf " edge1dy: %g   \n",edge1dy;

          edge2dy := sum(edge ee where original==2,
           sum(ee.facet ff, (ff.z*ee.x-ff.x*ee.z)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          printf " edge2dy: %g   \n",edge2dy;

          edge3dy := sum(edge ee where original==3,
           sum(ee.facet ff, (ff.z*ee.x-ff.x*ee.z)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          printf " edge3dy: %g   \n",edge3dy;

          edge4dy := sum(edge ee where original==4,
           sum(ee.facet ff, (ff.z*ee.x-ff.x*ee.z)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          printf " edge4dy: %g   \n",edge4dy;

          edge5dy := sum(edge ee where original==5,
           sum(ee.facet ff, (ff.z*ee.x-ff.x*ee.z)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          printf " edge5dy: %g   \n",edge5dy;
        }

read "adjoint.cmd"

// Call this to do adjoint transformation!
adj := { unset vertex constraint 1; 
         unset edge constraint 1;
         adjoint;
       }

// Applying constraints after adjointing
frame := {
    unfix vertices; unfix edges;
    aa := minimum(min(vertex,x-y),min(vertex,x+y)); set vertex x x-aa;
    aa := min(vertex,z-x); set vertex z z-aa;
    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 4; set ee constraint 4; };
    foreach edge ee where original==2 do
    { set ee constraint 5; set ee.vertex constraint 5; };
    foreach edge ee where original==3 do
    { set ee constraint 4; set ee.vertex constraint 4; };
    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 4; set ee.vertex constraint 4; };
    foreach edge ee where original==6 do
    { set ee constraint 5; set ee.vertex constraint 5; };
    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 6; set ee.vertex constraint 6; 
      set ee constraint 7; set ee.vertex constraint 7;
      fix ee; fix ee.vertex; 
    };
}



// To get true asize after evolving after adjointing
true_asize := { printf "True asize: %20.15f\n",
   sum(edge where original == 2,length)/sum(edge where original==1,length); }
true_bsize := { printf "True bsize: %20.15f\n",
   sum(edge where original == 3,length)/sum(edge where original==1,length); }
true_csize := { printf "True csize: %20.15f\n",
   sum(edge where original == 4,length)/sum(edge where original==1,length); }
true_dsize := { printf "True dsize: %20.15f\n",
   sum(edge where original == 5,length)/sum(edge where original==1,length); }
true_esize := { printf "True esize: %20.15f\n",
   sum(edge where original == 6,length)/sum(edge where original==1,length); }
true := {true_asize;true_bsize;true_csize;true_dsize;true_esize}

// command to show full cell in rhombic dodecahedron
showrhombic := { transform_expr "abcabcg"; transforms on;                
                 show edge where valence <= 1;
                 show_trans "R";
	      }
rhombic_edges := {
      va := new_vertex(0,0,0);
      vb := new_vertex(1,-1,1);
      newe1 := new_edge(va,vb);
      vc := new_vertex(0,0,2);
      newe2 := new_edge(vb,vc);
      set edge[newe1] bare;
      set edge[newe1] fixed;
      set edge[newe1] no_refine;
      set edge[newe2] bare;
      set edge[newe2] fixed;
      set edge[newe2] no_refine;

}

// command to show full cell in cube
showcube := { transform_expr "abcabc"; transforms on; 
              set edge color black;
              show edge where valence == 1;
              show_trans "R";
	      }

setcolor := { set facet frontcolor yellow }

// To show just one fundamental region, do "transforms off".
// To show tetrahedron outline again, do "show edges".

// A typical evolution
gogo := { gg; adj; frame; g 5; V; V; g 5; show_trans "R"; hessian; hessian;
        }

/* Commands:
   gogo - typical evolution; do before the following commands.
   showcube - display cubic unit cell
   showrhombic - show rhombic dodecahedron unit cell
   rhombic_edges - create edges outlining rhombus
   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.
*/