// sphere.fe // Spherical tank partially filled with liquid with contact angle. SYMMETRIC_CONTENT // natural for a sphere parameter rad = 1.00 // sphere radius parameter angle = 45 // contact angle, in degrees // cosine of contact angle #define WALLT cos(angle*pi/180) // density of body #define RHO 1.0 // Various expressions used in constraints #define wstuff (rad^2*z/(x^2 + y^2)/sqrt(x^2+y^2+z^2)) #define gstuff (G*RHO/8*rad^4*z^4/(x^2 + y^2)*(x^2+y^2+z^2)^2) #define gapstuff (G*RHO*rad^4/8*z^2/(x^2+y^2+z^2)^2) constraint 1 convex // the sphere, as wall for liquid formula: sqrt(x^2 + y^2 + z^2) = rad energy: e1: WALLT*wstuff*y - gstuff*y + G*RHO/8*y*z^2 - gapstuff*y e2: -WALLT*wstuff*x + gstuff*x - G*RHO/8*x*z^2 + gapstuff*x e3: 0 content: c1: -rad*wstuff*y/3 c2: rad*wstuff*x/3 c3: 0 constraint 2 // the sphere, for display only formula: sqrt(x^2 + y^2 + z^2) = rad vertices: 1 rad 0 0 constraint 1 2 0 rad 0 constraint 1 3 -rad 0 0 constraint 1 4 0 -rad 0 constraint 1 5 0 0 rad constraint 2 FIXED // to show back hemisphere 6 0 0 -rad constraint 2 FIXED 7 0 rad 0 constraint 2 FIXED 8 -rad 0 0 constraint 2 FIXED 9 0 -rad 0 constraint 2 FIXED edges: 1 1 2 constraint 1 2 2 3 constraint 1 3 3 4 constraint 1 4 4 1 constraint 1 5 5 9 constraint 2 FIXED // to show back hemisphere 6 9 6 constraint 2 FIXED 7 6 7 constraint 2 FIXED 8 7 5 constraint 2 FIXED 9 5 8 constraint 2 FIXED 10 9 8 constraint 2 FIXED 11 6 8 constraint 2 FIXED 12 7 8 constraint 2 FIXED faces: 1 1 2 3 4 2 5 10 -9 constraint 2 density 0 FIXED // to show back hemisphere 3 -10 6 11 constraint 2 density 0 FIXED 4 -11 7 12 constraint 2 density 0 FIXED 5 8 9 -12 constraint 2 density 0 FIXED bodies: // start with sphere half full 1 1 volume 2*pi*rad^3/3 volconst 2*pi*rad^3/3 density RHO read // Typical short evolution gogo := { g 5; r; g 5; r; g 5; r; g 10; conj_grad; g 20; } // Command to print out coordinates of points along cross-section // of the surface. This starts at vertex 3 and follows edges that // point in the positive x direction. section := { thisv := 3; nextv := thisv; do { printf "%9.6f %9.6f %9.6f\n",vertex[thisv].x,vertex[thisv].y, vertex[thisv].z; foreach vertex[thisv].edge ee do { if ee.x > 10*abs(ee.y) then { nextv := ee.vertex[2].id; break; } }; if nextv == thisv then break; // done thisv := nextv; } while 1; }