/[escript]/trunk/finley/test/python/slip_stress.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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
6  @var __copyright__: copyrights  @var __copyright__: copyrights
7  @var __license__: licence agreement  @var __license__: licence agreement
# 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"
14  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  __copyright__="""  Copyright (c) 2006, 2007 by ACcESS MNRF
15                      http://www.access.edu.au                      http://www.access.edu.au
16                  Primary Business: Queensland, Australia"""                  Primary Business: Queensland, Australia"""
17  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
# Line 27  from esys.escript import * Line 24  from esys.escript import *
24  from esys.escript.pdetools import SaddlePointProblem  from esys.escript.pdetools import SaddlePointProblem
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
86  pde.setValue(A=A,q=mask,Y=-kronecker(Function(dom))[d-1]*g*rho,X=sigma0)  pde.setValue(A=A,q=mask,Y=-kronecker(Function(dom))[d-1]*g*rho,X=sigma0)
87    
88    print "=== solve pde === "
89  u=pde.getSolution(verbose=True)  u=pde.getSolution(verbose=True)
90    
91    print "=== calculate stress ==="
92  g_s=grad(u)  g_s=grad(u)
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    
120    # add some text
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()
126    t2.shadowOn()
127    
128    # Render the object.
129    s.render("stress_contour.jpg")
130    print "=== finished rendering ==="
Legend:
Removed from v.1177  
changed lines
  Added in v.1178
  ViewVC Help
Powered by ViewVC 1.1.26