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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1178 - (show annotations)
Wed May 30 05:03:17 2007 UTC (13 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 3907 byte(s)
bug in visualization interface test fixed.
1 # $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 __copyright__=""" Copyright (c) 2006, 2007 by ACcESS MNRF
15 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 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):
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 # fault location:
49
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 s=[0.,1.,0.]
53
54 print "=== generate mesh over %s x %s x %s ==="%(width,width,height)
55 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 "=== prepare PDE coefficients ==="
59 # 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 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)
90
91 print "=== calculate stress ==="
92 g_s=grad(u)
93 sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d)
94
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