Diff of /trunk/finley/test/python/slip_stress.py

revision 1177 by gross, Thu Nov 9 07:37:40 2006 UTC revision 1178 by gross, Wed May 30 05:03:17 2007 UTC
# Line 1  Line 1
1  # \$Id\$  # \$Id\$

2  """  """
3  calculation of the stress distribution around a fault from the slip on the fault  calculation of the stress distribution around a fault from the slip on the fault
4
e.g. use slip_stress_mesh.py to generate mesh

5  @var __author__: name of author  @var __author__: name of author
# Line 14  e.g. use slip_stress_mesh.py to generate Line 11  e.g. use slip_stress_mesh.py to generate
11  """  """
12
13  __author__="Lutz Gross, Louise Kettle"  __author__="Lutz Gross, Louise Kettle"
15                      http://www.access.edu.au                      http://www.access.edu.au
# Line 27  from esys.escript import * Line 24  from esys.escript import *
25  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
26  from esys.finley import Brick  from esys.finley import Brick
27    from esys.pyvisi import Scene, DataCollector, Contour, Camera, Velocity, Text2D, LocalPosition
28    from esys.pyvisi.constant import *
29    import os
30
31    JPG_RENDERER = Renderer.OFFLINE_JPG
32    ne = 20          # number of elements in spatial direction
33
34  def lockToGrid(x,l,n):  def lockToGrid(x,l,n):
35     h=l/n     h=l/n
print x,l,n,"->",h
36     return int(x/h+0.5)*h     return int(x/h+0.5)*h
37
38
ne = 15
39  width  = 100000.  width  = 100000.
40  height =  30000.  height =  30000.
41
# Line 43  rho=0. Line 44  rho=0.
44  lmbd=1.7e11  lmbd=1.7e11
45  mu=1.7e11  mu=1.7e11
46  g=9.81  g=9.81
print "total ne = ",ne*ne_w*ne_w
47
48  # fault:  # fault location:
49
50  fstart =  [lockToGrid(50000.0,width,ne_w), lockToGrid(40000.0,width,ne_w), lockToGrid(8000.,height,ne)]  fstart =  [lockToGrid(50000.0,width,ne_w), lockToGrid(40000.0,width,ne_w), lockToGrid(8000.,height,ne)]
51  fend =  [lockToGrid(50000.0,width,ne_w), lockToGrid(60000.0,width,ne_w), lockToGrid(20000.,height,ne)]  fend =  [lockToGrid(50000.0,width,ne_w), lockToGrid(60000.0,width,ne_w), lockToGrid(20000.,height,ne)]
52  s=[0.,1.,0.]  s=[0.,1.,0.]
53
54  # fstart =  [lockToGrid(30000.0,width,ne_w), lockToGrid(30000.0,width,ne_w), lockToGrid(20000.,height,ne)]  print "=== generate mesh over %s x %s x %s ==="%(width,width,height)
# fend =  [lockToGrid(70000.0,width,ne_w), lockToGrid(70000.0,width,ne_w), lockToGrid(20000.,height,ne)]
# s=[1.,0.,0.]

55  dom=Brick(l0=width, l1=width, l2=height, n0=ne_w, n1=ne_w, n2=ne)  dom=Brick(l0=width, l1=width, l2=height, n0=ne_w, n1=ne_w, n2=ne)
56    print "  total number of elements = ",ne*ne_w*ne_w
57
58  print "fstart:",fstart  print "=== prepare PDE coefficients ==="
print "fend:",fend
print "slip ",s

59  # fixed displacements:  # fixed displacements:
60  d=dom.getDim()  d=dom.getDim()
61  x=dom.getX()  x=dom.getX()
# Line 72  p=[0.,0.,0.] Line 67  p=[0.,0.,0.]
67  q=1.  q=1.
68  for i in range(d):  for i in range(d):
69       p[i]=fend[i]-fstart[i]       p[i]=fend[i]-fstart[i]
print i,p[i]
70       if abs(p[i])>0:       if abs(p[i])>0:
71           x_hat[i]/=p[i]           x_hat[i]/=p[i]
72           q=whereNonNegative(x_hat[i])*whereNonPositive(x_hat[i]-1.)*x_hat[i]*(1.-x_hat[i])*4.*q           q=whereNonNegative(x_hat[i])*whereNonPositive(x_hat[i]-1.)*x_hat[i]*(1.-x_hat[i])*4.*q
# Line 90  for i in range(d): Line 84  for i in range(d):
84        A[i,j,i,j] += mu        A[i,j,i,j] += mu
85        A[i,i,j,j] += lmbd        A[i,i,j,j] += lmbd
87
88    print "=== solve pde === "
89  u=pde.getSolution(verbose=True)  u=pde.getSolution(verbose=True)
90
91    print "=== calculate stress ==="
93  sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d)  sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d)
94  saveVTK("dis.xml",disp=u,sigma=sigma,cfs=sigma[0,1]-0.4*sigma[0,0])
95    print "=== start rendering ==="
96    # Create a Scene.
97    s = Scene(renderer = JPG_RENDERER, num_viewport = 1, x_size = 800, y_size = 800)
98    dc1 = DataCollector(source = Source.ESCRIPT)
99    dc1.setData(disp=u,sigma=sigma,cfs=sigma[0,1]-0.4*sigma[0,0])
100
101    # Create a Contour.
102    ctr1 = Contour(scene = s, data_collector = dc1, viewport = Viewport.SOUTH_WEST,
103            lut = Lut.COLOR, cell_to_point = True, outline = True)
104    ctr1.generateContours(contours = 10)
105    ctr1.setOpacity(0.5)
106
107    # Create velocity
108    vopc1 = Velocity(scene = s, data_collector = dc1,
109            color_mode = ColorMode.VECTOR,
110            arrow = Arrow.THREE_D, lut = Lut.COLOR, cell_to_point = False,
111            outline = False)
112    vopc1.setScaleFactor(scale_factor = 10000.)
113    vopc1.setRatio(2)
114    vopc1.randomOn()
115
116    # Create a Camera.
117    cam1 = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
118    cam1.elevation(angle = -50)
119
121    t2 = Text2D(scene = s, text = "CFS and displacement around vertical fault")
122    t2.setPosition(LocalPosition(200,30))
123    t2.setColor(color = Color.BLACK)
124    t2.setFontSize(size = 20)
125    t2.setFontToArial()