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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1178 - (hide annotations)
Wed May 30 05:03:17 2007 UTC (14 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 3907 byte(s)
bug in visualization interface test fixed.
1 gross 896 # $Id$
2     """
3     calculation of the stress distribution around a fault from the slip on the fault
4    
5     @var __author__: name of author
6     @var __copyright__: copyrights
7     @var __license__: licence agreement
8     @var __url__: url entry point on documentation
9     @var __version__: version
10     @var __date__: date of the version
11     """
12    
13     __author__="Lutz Gross, Louise Kettle"
14 gross 1178 __copyright__=""" Copyright (c) 2006, 2007 by ACcESS MNRF
15 gross 896 http://www.access.edu.au
16     Primary Business: Queensland, Australia"""
17     __license__="""Licensed under the Open Software License version 3.0
18     http://www.opensource.org/licenses/osl-3.0.php"""
19     __url__="http://www.iservo.edu.au/esys"
20     __version__="$Revision$"
21     __date__="$Date$"
22    
23     from esys.escript import *
24     from esys.escript.pdetools import SaddlePointProblem
25     from esys.escript.linearPDEs import LinearPDE
26     from esys.finley import Brick
27 gross 1178 from esys.pyvisi import Scene, DataCollector, Contour, Camera, Velocity, Text2D, LocalPosition
28     from esys.pyvisi.constant import *
29     import os
30 gross 896
31 gross 1178 JPG_RENDERER = Renderer.OFFLINE_JPG
32     ne = 20 # number of elements in spatial direction
33    
34 gross 896 def lockToGrid(x,l,n):
35     h=l/n
36     return int(x/h+0.5)*h
37    
38    
39     width = 100000.
40     height = 30000.
41    
42     ne_w=int((ne/height)*width+0.5)
43     rho=0.
44     lmbd=1.7e11
45     mu=1.7e11
46     g=9.81
47    
48 gross 1178 # fault location:
49 gross 896
50     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)]
52 gross 897 s=[0.,1.,0.]
53 gross 896
54 gross 1178 print "=== generate mesh over %s x %s x %s ==="%(width,width,height)
55 gross 896 dom=Brick(l0=width, l1=width, l2=height, n0=ne_w, n1=ne_w, n2=ne)
56 gross 1178 print " total number of elements = ",ne*ne_w*ne_w
57 gross 896
58 gross 1178 print "=== prepare PDE coefficients ==="
59 gross 896 # fixed displacements:
60     d=dom.getDim()
61     x=dom.getX()
62     mask=whereZero(x[d-1]-inf(x[d-1]))*numarray.ones((d,))
63    
64     # map the fault into the unit square:
65     x_hat=x-fstart
66     p=[0.,0.,0.]
67     q=1.
68     for i in range(d):
69     p[i]=fend[i]-fstart[i]
70     if abs(p[i])>0:
71     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
73     else:
74     flip=i
75     q=whereZero(x_hat[i],1.e-3*Lsup(x_hat[i]))*q
76     g_s=grad(q*s)
77     sigma0=(mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d))*(whereNegative(interpolate(x_hat[flip],g_s.getFunctionSpace()))-0.5)
78     pde=LinearPDE(dom)
79     pde.setSymmetryOn()
80     A =Tensor4(0.,Function(dom))
81     for i in range(d):
82     for j in range(d):
83     A[i,j,j,i] += mu
84     A[i,j,i,j] += mu
85     A[i,i,j,j] += lmbd
86 gross 897 pde.setValue(A=A,q=mask,Y=-kronecker(Function(dom))[d-1]*g*rho,X=sigma0)
87 gross 1178
88     print "=== solve pde === "
89 gross 896 u=pde.getSolution(verbose=True)
90 gross 1178
91     print "=== calculate stress ==="
92 gross 896 g_s=grad(u)
93 gross 897 sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d)
94 gross 1178
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 ==="

  ViewVC Help
Powered by ViewVC 1.1.26