// ringblob.fe

// Toroidal liquid ring on a rotating rod in weightlessness.
// Half of full torus
// Using second periodic constraint surface intersecting rod to
// confine vertices on rod to vertical motion. 

// This permits drawing both halves of the ring
view_transforms 1
1 0 0 0
0 1 0 0
0 0 -1 0 
0 0 0 1

// Basic parameters.  These may be adjusted at runtime with the
// 'A' command.  Only spin is being adjusted in these experiments.
parameter rodr = 1  // rod radius
parameter spin = 0.5 // angular velocity
parameter angle = 90 // internal contact angle with rod
#define rode (-cos(angle*pi/180))  // liquid-rod contact energy
parameter dens = 1  // density of liquid, negative for bubble

// spin centripetal energy
quantity centrip energy global_method facet_vector_integral
vector_integrand:
q1: 0
q2: 0
q3: -0.5*dens*spin*spin*(x^2+y^2)*z

// y moment, for detecting instability
quantity ymoment info_only global_method facet_vector_integral
vector_integrand:
q1: 0
q2: 0
q3: y*z

// Constraint for vertices and edges confined to rod surface,
// with integral for blob area on rod
constraint 1
formula: x^2 + y^2 = rodr^2
energy: 
e1: rode*z*y/rodr
e2: -rode*z*x/rodr
e3: 0

// Horizontal symmetry plane
constraint 2 
formula: z = 0

// Rod surface as one-sided constraint, to keep stuff from caving in
constraint 3 nonnegative 
formula: x^2 + y^2 = rodr^2

// Constraint to force vertices on rod to move only vertically.
// Expressed in periodic form, so one constraint fits arbitrarily
// many vertices. Note offset to pi/6 to avoid difficulties with
// modulus discontinuity at 0.
parameter pp = pi/2    /* to be halved each refinement */
parameter qq = pi/6    /* to be halved each refinement */
constraint 4 
formula:  (atan2(y,x)+pi/6) % pp = qq

//initial dimensions
#define ht 2
#define wd 3

vertices
1  0   -wd 0  constraints 2    // equatorial vertices
2  wd    0 0  constraints 2
3  0    wd 0  constraints 2
4  -wd   0 0  constraint 2
5  0   -wd ht                 // upper outer corners
6  wd    0 ht
7  0    wd ht
8  -wd   0 ht 
9  0 -rodr ht constraints 1,4   // vertices on rod
10 rodr  0 ht constraints 1,4
11 0  rodr ht constraints 1,4
12 -rodr  0 ht constraints 1,4

edges
1   1 2 constraint 2  // equatorial edges
2   2 3 constraint 2
3   3 4 constraint 2
4   4 1 constraint 2
5   5 6               // upper outer edges 
6   6 7
7   7 8
8   8 5
9   9  10 constraint 1,4   // edges on rod
10  10 11 constraint 1,4
11  11 12 constraint 1,4
12  12  9 constraint 1,4
13   1  5        // vertical outer edges
14   2  6
15   3  7
16   4  8
17   5  9 // constraint 3       // cutting up top face
18   6 10 // constraint 3
19   7 11 // constraint 3
20   8 12 // constraint 3

faces  /* given by oriented edge loop */
1   1 14 -5 -13  // side faces
2   2 15 -6 -14
3   3 16 -7 -15
4   4 13 -8 -16
5   5 18 -9 -17  // constraint 3 // top faces, must stay out of rod
6   6 19 -10 -18 // constraint 3
7   7 20 -11 -19 // constraint 3
8   8 17 -12 -20 // constraint 3

bodies  /* one body, defined by its oriented faces */
1   1 2 3 4 5 6 7 8  volume 25.28

read  // some initializations
hessian_normal    // so normal perturbations only
linear_metric     // so eigenvalues independent of refinement

// special refinement command
r :::= { autorecalc off;  pp := pp/2; qq := qq % pp;  r; autorecalc on; }

// Typical evolution.  This is for spin .5, which is unstable, but
// the initial datafile is symmetric enough the instabilities don't grow
// fast enough to matter before hessian wipes them out.
gogo := { pp := pp/2; qq := qq % pp;
          refine edge where on_constraint 1 or on_constraint 2;
          u; t .5; g 5; u; g 5;
          r; g 5; V; u; V; u; V; g 5; 
          hessian; // warning about indefinite hessian; harmless here
          r; u; V; g 20; V; hessian; hessian;
        }
// After above gogo, you can use hessian_menu to look at the modes of 
// instability.  Can't do a script here, since hessian_menu has its own
// prompt, but here's a sequence of commands:
// Enter command: hessian_menu
// Choice(?,1,2,3,4,7,9,B,C,E,F,L,R,Z,X,P,V,S,Y,U,M,G,0,q): 1
// Choice(?,1,2,3,4,7,9,B,C,E,F,L,R,Z,X,P,V,S,Y,U,M,G,0,q): z
// Ritz probe value: -1
// Number of eigenvalues: 5
// Eigencounts:   0 <,  0 ==,  287 >
//   1.   -0.9127209565686
//   2.   -0.9127209565686
//   3.   -0.4572157022779
//   4.   -0.4555002348136
//   5.    0.2566767756229
// Iterations: 13. Total eigenvalue changes in last iteration: 5.9012795e-011
// Choice(?,1,2,3,4,7,9,B,C,E,F,L,R,Z,X,P,V,S,Y,U,M,G,0,q): x
// Eigenvector number (1 to 5):
// Choice(?,1,2,3,4,7,9,B,C,E,F,L,R,Z,X,P,V,S,Y,U,M,G,0,q): 4
// Step size: 1
// 1. area: 36.914711203537244 energy: 23.145324165369388
// Choice(?,1,2,3,4,7,9,B,C,E,F,L,R,Z,X,P,V,S,Y,U,M,G,0,q): 7
// Choice(?,1,2,3,4,7,9,B,C,E,F,L,R,Z,X,P,V,S,Y,U,M,G,0,q): x
// Eigenvector number (1 to 5): 3
// Choice(?,1,2,3,4,7,9,B,C,E,F,L,R,Z,X,P,V,S,Y,U,M,G,0,q): 4
// Step size: 1
// 1. area: 37.152556797943689 energy: 23.362763236279037
// Choice(?,1,2,3,4,7,9,B,C,E,F,L,R,Z,X,P,V,S,Y,U,M,G,0,q): 7
// Choice(?,1,2,3,4,7,9,B,C,E,F,L,R,Z,X,P,V,S,Y,U,M,G,0,q): q
//   1. area:  37.0104558459429 energy:  23.5761866980545
// Enter command: