// tankex.fe

// Equilibrium shape of liquid in flat-ended cylindrical tank.

// Tank has axis along y-axis and flat bottom in x-z plane.  This
// is so gravity acting vertically draws liquid toward wall.

// Straight edges cannot conform exactly to curved walls.  
// We need to give them an area so that area cannot shrink by straght edges 
// pulling away from the walls. The gaps are also accounted for
// in volume and gravitational energy.

SYMMETRIC_CONTENT     // for volume calculations

// Contact angles, initially for 30 degrees (seems unstable for angles > 45).
PARAMETER ENDT  = cos(30*pi/180)   /* surface tension of uncovered base   */
PARAMETER WALLT = cos(30*pi/180)   /* surface tension of uncovered side wall */

// Gravity components
PARAMETER GY = 0
PARAMETER GZ = -.3

SPRING_CONSTANT 1  // for most accurate gap areas for constraint 2

#define  TR  1.00       /* tank radius */
#define  RHO 1.00       /* fuel density */

constraint 1  // flat base
function:  y = 0
energy:               // for contact energy line integral
e1:   -ENDT*z
e2:   0
e3:   0

#define wstuff (WALLT*TR*y/(x^2 + z^2))   // common wall area term
#define vstuff (TR^2/3*y/(x^2 + z^2))     // common wall volume term
#define gstuff (GY*TR^2*y^2/4/(x^2 + z^2) + GZ*TR^3*y*z/3/(x^2+z^2)**1.5)
                                          // common gap gravity term

constraint 2  CONVEX         // cylindrical wall
function:  x^2 + z^2 = TR^2
energy:                // for contact energy and gravity
e1:  -wstuff*z  + RHO*GY*y^2*z/4 + RHO*GZ*y*z^2/3 - RHO*gstuff*z
e2:    0
e3:   wstuff*x  - RHO*GY*y^2*x/4 - RHO*GZ*y*z*x/3 + RHO*gstuff*x
content:               // so volumes calculated correctly
c1:   vstuff*z - z*y/6 + vstuff*z/2
c2:    0
c3:  -vstuff*x + x*y/6 - vstuff*x/2

// named quantity for arbitrary direction gravity on facets
quantity arb_grav energy method facet_vector_integral global
vector_integrand:
q1:  0
q2:  -RHO*GY*y^2/2 - RHO*GZ*y*z
q3:  0


// Now the specification of the initial shape

vertices
1    0.5  0.0  0.5  constraint 1
2    0.5  0.0 -0.5  constraint 1
3   -0.5  0.0 -0.5  constraint 1
4   -0.5  0.0  0.5  constraint 1
5    1.0  0.5  0.0  constraint 2
6    0.0  0.5 -1.0  constraint 2
7   -1.0  0.5  0.0  constraint 2
8    0.0  0.5  1.0  constraint 2

edges
1    2  1  constraint 1
2    1  4  constraint 1
3    4  3  constraint 1
4    3  2  constraint 1
5    5  6  constraint 2
6    6  7  constraint 2
7    7  8  constraint 2
8    8  5  constraint 2
9    1  8
10   1  5
11   2  5
12   2  6
13   3  6
14   3  7
15   4  7
16   4  8

faces
1   13   6 -14
2    3  14 -15 
3   15   7 -16
4    2  16  -9
5    9   8 -10
6    1  10 -11
7   11   5 -12
8    4  12 -13

bodies
1    1 2 3 4 5 6 7 8   volume 0.6  density 0  /* no default gravity */ 


read

// typical evolution
gogo := { r; g5; r; g 5; U; g 70; t .01; g 10; t .05; g 50; t .05;
	  r; g 20;
	}