| 23 |
from esys.finley import Brick |
from esys.finley import Brick |
| 24 |
import time |
import time |
| 25 |
|
|
| 26 |
output=False |
output=True |
| 27 |
n_end=100 |
n_end=10000 |
| 28 |
|
|
| 29 |
resolution=4000. # number of elements per m |
resolution=10000. # number of elements per m |
| 30 |
l=80000. # width and length m |
l=80000. # width and length m |
| 31 |
h=30000. # height in m |
h=30000. # height in m |
| 32 |
|
o=1 # element order |
| 33 |
|
|
| 34 |
rho_bedrock=8e3 |
rho_bedrock=8e3 |
| 35 |
mu_bedrock=1.7e11 |
mu_bedrock=1.7e11 |
| 77 |
|
|
| 78 |
""" |
""" |
| 79 |
global netotal |
global netotal |
| 80 |
ne_l=int(l/resolution+0.5) |
v_p_bedrock=sqrt((2*mu_bedrock+lambda_bedrock)/rho_bedrock) |
| 81 |
ne_h=int(h/resolution+0.5) |
v_p_sand=sqrt((2*mu_sand+lambda_sand)/rho_sand) |
| 82 |
netotal=ne_l*ne_l*ne_h |
v_p_water=sqrt((2*mu_water+lambda_water)/rho_water) |
| 83 |
if output: print "grid : %s x %s x %s"%(ne_l,ne_l,ne_h) |
|
| 84 |
dom=Brick(ne_l,ne_l,ne_h,l0=l,l1=l,l2=h,order=1) |
print v_p_bedrock,v_p_sand,v_p_water |
| 85 |
|
|
| 86 |
|
ne_l_x_bed=int((l-d0_sand-l_x_water)/resolution+0.5) |
| 87 |
|
ne_l_y_bed=int((l-d0_sand-l_y_water)/resolution+0.5) |
| 88 |
|
ne_h_bed=int((h-d_sand-h_water)/resolution+0.5) |
| 89 |
|
|
| 90 |
|
ne_l_sand=int(d0_sand*v_p_bedrock/v_p_sand/resolution+0.5) |
| 91 |
|
ne_h_sand=int(d_sand*v_p_bedrock/v_p_sand/resolution+0.5) |
| 92 |
|
|
| 93 |
|
ne_l_x_water=int(l_x_water*v_p_bedrock/v_p_water/resolution+0.5) |
| 94 |
|
ne_l_y_water=int(l_y_water*v_p_bedrock/v_p_water/resolution+0.5) |
| 95 |
|
ne_h_water=int(h_water*v_p_bedrock/v_p_water/resolution+0.5) |
| 96 |
|
|
| 97 |
|
ne_l_x=ne_l_x_bed+ne_l_sand+ne_l_x_water |
| 98 |
|
ne_l_y=ne_l_y_bed+ne_l_sand+ne_l_y_water |
| 99 |
|
ne_h=ne_h_bed+ne_h_sand+ne_h_water |
| 100 |
|
|
| 101 |
|
print ne_l_x,ne_l_x_bed,ne_l_sand,ne_l_x_water |
| 102 |
|
print ne_l_y,ne_l_y_bed,ne_l_sand,ne_l_y_water |
| 103 |
|
print ne_h,ne_h_bed,ne_h_sand,ne_h_water |
| 104 |
|
|
| 105 |
|
netotal=ne_l_x*ne_l_y*ne_h |
| 106 |
|
if output: print "grid : %s x %s x %s (%s elements)"%(ne_l_x,ne_l_y,ne_h,netotal) |
| 107 |
|
dom=Brick(ne_l_x,ne_l_y,ne_h,l0=o*ne_l_x,l1=o*ne_l_y,l2=o*ne_h,order=o) |
| 108 |
|
x=dom.getX() |
| 109 |
|
|
| 110 |
|
x0=x[0] |
| 111 |
|
x0_new = (l-d0_sand-l_x_water)/(o*ne_l_x_bed)*x0 |
| 112 |
|
msk=whereNonPositive(x0-o*ne_l_x_bed) |
| 113 |
|
x0_new = x0_new * msk + (d0_sand/(o*ne_l_sand)*(x0-o*ne_l_x_bed)+l-d0_sand-l_x_water)*(1.-msk) |
| 114 |
|
msk=whereNonPositive(x0-o*(ne_l_x_bed+ne_l_sand)) |
| 115 |
|
x0_new = x0_new * msk + (l_x_water/(o*ne_l_x_water)*(x0-o*(ne_l_x_bed+ne_l_sand))+l-l_x_water)*(1.-msk) |
| 116 |
|
|
| 117 |
|
x1=x[1] |
| 118 |
|
x1_new = (l-d0_sand-l_y_water)/(o*ne_l_y_bed)*x1 |
| 119 |
|
msk=whereNonPositive(x1-(o*ne_l_y_bed)) |
| 120 |
|
x1_new = x1_new * msk + (d0_sand/(o*ne_l_sand)*(x1-o*ne_l_y_bed)+l-d0_sand-l_y_water)*(1.-msk) |
| 121 |
|
msk=whereNonPositive(x1-o*(ne_l_y_bed+ne_l_sand)) |
| 122 |
|
x1_new = x1_new * msk + (l_y_water/(o*ne_l_y_water)*(x1-o*(ne_l_y_bed+ne_l_sand))+l-l_y_water)*(1.-msk) |
| 123 |
|
|
| 124 |
|
x2=x[2] |
| 125 |
|
x2_new = (h-d_sand-h_water)/(o*ne_h_bed)*x2 |
| 126 |
|
msk=whereNonPositive(x2-(o*ne_h_bed+1)) |
| 127 |
|
x2_new = x2_new * msk + (d_sand/(o*ne_h_sand)*(x2-(o*ne_h_bed+1))+h-d_sand-h_water)*(1.-msk) |
| 128 |
|
msk=whereNonPositive(x2-(o*(ne_h_bed+ne_h_sand)+1)) |
| 129 |
|
x2_new = x2_new * msk + (h_water/(o*ne_h_water)*(x2-(o*(ne_h_bed+ne_h_sand)+1))+h-h_water)*(1.-msk) |
| 130 |
|
|
| 131 |
|
dom.setX(x0_new*[1,0,0]+x1_new*[0,1,0]+x2_new*[0,0,1]) |
| 132 |
|
|
| 133 |
fs=Function(dom) |
fs=Function(dom) |
| 134 |
x=Function(dom).getX() |
x=Function(dom).getX() |
Parent Directory
| 
Revision Log
| 
Patch