// t0.fe   - bizarre symmetric starting point with s1xs1 topology
// Test of facet_2form_sq_integral method, with corrected symplectic area form
//  Andrew Hanson, Indiana University,  September 2001

space_dimension 4

quantity  mycurv  ENERGY modulus 0 global_method star_sq_mean_curvature

// symplectic area
// Correspondence: z1 = (x1,x2)  z2 = (x3,x4)
#define DENOM ((x1^2+x2^2+x3^2+x4^2)^2)
quantity symplectic_sq energy method facet_2form_sq_integral global
form_integrand:
q1: -2*(x3^2 + x4^2)/DENOM    // dx1 wedge dx2 term
q2:  2*(x2*x3-x1*x4)/DENOM    // dx1 wedge dx3 term
q3:  2*(x1*x3+x2*x4)/DENOM    // dx1 wedge dx4 term
q4: -2*(x1*x3+x2*x4)/DENOM    // dx2 wedge dx3 term
q5:  2*(x2*x3-x1*x4)/DENOM    // dx2 wedge dx4 term
q6: -2*(x1^2 + x2^2)/DENOM    // dx3 wedge dx4 term


#define a  0.5
#define b  0.866025

vertices

//  first tri
1   0.  a   0.   b
2   0. -1.  0.   0.
3   0.  a   0.  -b

//  second tri
4   a   0.  b   0.
5  -1.  0.  0.  0.
6   a   0. -b   0.

//  symmetric point
7   0  0  0  0
8   0  0  0  0
9   0  0  0  0

edges

1   1 4  
2   1 5
3   1 6
4   2 4
5   2 5
6   2 6
7   3 4
8   3 5
9   3 6
10  4 7
11  4 8
12  4 9
13  5 7
14  5 8
15  5 9
16  6 7
17  6 8
18  6 9
19  7 1
20  7 2
21  7 3
22  8 1
23  8 2
24  8 3
25  9 1
26  9 2
27  9 3

faces
1  22 2 14     
2  -2 -19 -13  
3  1  10 19     
4  -4 -20 -10    
5  23 4 11      
6  -14 -5  -23   
7  17 24 9      
8  -3 -22 -17    
9  25 3 18      
10  -12 -1 -25    
11   7 12 27    
12  -24 -11 -7    
13  21 8 13    
14  -16 -9 -21    
15  6 16 20    
16  -26 -18 -6    
17  15 26 5    
18  -8 -27 -15    

read
set facet tension 0

gogo := { r; g 10; r; g 10; r; g 10 }

// to see it evolve by square curvature
go_curv := { mycurv.modulus := 1; symplectic_sq.modulus := 0; recalc; gogo; 
             symplectic_sq.modulus := 1; recalc; v; }

// to see it evolve by symplectic_sq
go_symp := { mycurv.modulus := 0; symplectic_sq.modulus := 1; recalc; 
             scale_limit := 10;
             r; g 10; r; g 10; conj_grad; g 100; r; g 100;
             mycurv.modulus := 1; recalc; v;
           }