24 |
|
|
25 |
require 'check.ph'; |
require 'check.ph'; |
26 |
|
|
27 |
# Empty Povray scene |
Run('','elasticity.ff'); |
|
WriteFile('test.pov',<< 'EOF'); |
|
|
box{ |
|
|
<0,-0.5,-0.5>,<5,0.5,0.5> |
|
|
pigment{color rgb <1,0,0>} |
|
|
} |
|
|
EOF |
|
|
|
|
|
Run(<< 'EOF'); |
|
|
vector n = (50,10,10); |
|
|
vector a = (0,-0.5,-0.5); |
|
|
vector b = (5, 0.5, 0.5); |
|
|
|
|
|
mesh M = structured(n,a,b); |
|
|
scene S = pov("test.pov"); |
|
|
domain Omega = domain(S,inside(<1,0,0>)); |
|
|
|
|
|
double L = 1000; // Lambda |
|
|
double Mu = 500; |
|
|
|
|
|
solve (u1,u2,u3) in Omega by M |
|
|
krylov(precond=diagonal),memory(matrix=none) |
|
|
{ |
|
|
test(v1,v2,v3) |
|
|
|
|
|
int(Mu*dx(u1)*dx(v1))+int(Mu*dx(u2)*dx(v2))+int(Mu*dx(u3)*dx(v3)) |
|
|
+int(Mu*dy(u1)*dy(v1))+int(Mu*dy(u2)*dy(v2))+int(Mu*dy(u3)*dy(v3)) |
|
|
+int(Mu*dz(u1)*dz(v1))+int(Mu*dz(u2)*dz(v2))+int(Mu*dz(u3)*dz(v3)) |
|
|
|
|
|
+int(Mu*dx(u1)*dx(v1))+int(Mu*dx(u2)*dy(v1))+int(Mu*dx(u3)*dz(v1)) |
|
|
+int(Mu*dy(u1)*dx(v2))+int(Mu*dy(u2)*dy(v2))+int(Mu*dy(u3)*dz(v2)) |
|
|
+int(Mu*dz(u1)*dx(v3))+int(Mu*dz(u2)*dy(v3))+int(Mu*dz(u3)*dz(v3)) |
|
|
|
|
|
+int(L*dx(u1)*dx(v1))+int(L*dx(u1)*dy(v2))+int(L*dx(u1)*dz(v3)) |
|
|
+int(L*dy(u2)*dx(v1))+int(L*dy(u2)*dy(v2))+int(L*dy(u2)*dz(v3)) |
|
|
+int(L*dz(u3)*dx(v1))+int(L*dz(u3)*dy(v2))+int(L*dz(u3)*dz(v3)) |
|
|
|
|
|
= int(-v2); // gravity |
|
|
|
|
|
u1=0 on M xmin; |
|
|
u2=0 on M xmin; |
|
|
u3=0 on M xmin; |
|
|
u1=0 on M xmax; |
|
|
u2=0 on M xmax; |
|
|
u3=0 on M xmax; |
|
|
}; |
|
|
|
|
|
save(opendx,"u1.dat",u1,M); |
|
|
save(opendx,"u2.dat",u2,M); |
|
|
save(opendx,"u3.dat",u3,M); |
|
|
EOF |
|
28 |
|
|
29 |
# Compares the resulting output with a reference |
# Compares the resulting output with a reference |
30 |
|
|
31 |
VectorCheck('test.bb','elasticity-ref.bb',1e-6); |
VectorCheck('elasticity-u1.bb','elasticity-u1-ref.bb',1e-6); |
32 |
|
VectorCheck('elasticity-u2.bb','elasticity-u2-ref.bb',1e-6); |
33 |
|
VectorCheck('elasticity-u3.bb','elasticity-u3-ref.bb',1e-6); |
34 |
|
|
35 |
Report(); |
Report(); |