1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
// 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;
}
|
