# Diff of /trunk/finley/test/python/seismic_wave.py

revision 843 by gross, Thu Sep 7 05:55:12 2006 UTC revision 844 by gross, Thu Sep 7 08:32:01 2006 UTC
# Line 23  from esys.escript.linearPDEs import Line Line 23  from esys.escript.linearPDEs import Line
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
# Line 76  def getDomain(): Line 77  def getDomain():
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()

Legend:
 Removed from v.843 changed lines Added in v.844