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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4154 - (show annotations)
Tue Jan 22 09:30:23 2013 UTC (6 years, 7 months ago) by jfenwick
File MIME type: text/x-python
File size: 4800 byte(s)
Round 1 of copyright fixes
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17 http://www.uq.edu.au
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 """
24 calculation of the stress distribution around a fault from the slip on the fault
25
26 :var __author__: name of author
27 :var __copyright__: copyrights
28 :var __license__: licence agreement
29 :var __url__: url entry point on documentation
30 :var __version__: version
31 :var __date__: date of the version
32 """
33
34 __author__="Lutz Gross, Louise Kettle"
35
36 from esys.escript import *
37 from esys.escript.pdetools import SaddlePointProblem
38 from esys.escript.linearPDEs import LinearPDE
39 from esys.dudley import Brick
40 from esys.pyvisi import Scene, DataCollector, Contour, Camera, Velocity, Text2D, LocalPosition, Legend
41 from esys.pyvisi.constant import *
42 import os
43
44 JPG_RENDERER = Renderer.ONLINE_JPG
45 ne = 15 # number of elements in spatial direction
46 s=[1.,1.,0.]
47
48 XSIZE=900
49 YSIZE=800
50 VIS_RANGE_STRESS=[-3e4, 20.e4]
51
52 def lockToGrid(x,l,n):
53 h=l/n
54 return int(x/h+0.5)*h
55
56
57 width = 100000.
58 height = 30000.
59
60 ne_w=int((ne/height)*width+0.5)
61 rho=0.
62 lmbd=1.7e11
63 mu=1.7e11
64 g=9.81
65
66 # fault location:
67
68 fstart = [lockToGrid(50000.0,width,ne_w), lockToGrid(40000.0,width,ne_w), lockToGrid(8000.,height,ne)]
69 fend = [lockToGrid(50000.0,width,ne_w), lockToGrid(60000.0,width,ne_w), lockToGrid(20000.,height,ne)]
70
71 print("=== generate mesh over %s x %s x %s ==="%(width,width,height))
72 dom=Brick(l0=width, l1=width, l2=height, n0=ne_w, n1=ne_w, n2=ne)
73 print(" total number of elements = ",ne*ne_w*ne_w)
74
75 print("=== prepare PDE coefficients ===")
76 # fixed displacements:
77 d=dom.getDim()
78 x=dom.getX()
79 mask=whereZero(x[d-1]-inf(x[d-1]))*numarray.ones((d,))
80
81 # map the fault into the unit square:
82 x_hat=x-fstart
83 p=[0.,0.,0.]
84 q=1.
85 for i in range(d):
86 p[i]=fend[i]-fstart[i]
87 if abs(p[i])>0:
88 x_hat[i]/=p[i]
89 q=whereNonNegative(x_hat[i])*whereNonPositive(x_hat[i]-1.)*x_hat[i]*(1.-x_hat[i])*4.*q
90 else:
91 flip=i
92 q=whereZero(x_hat[i],1.e-3*Lsup(x_hat[i]))*q
93 g_s=grad(q*s)
94 sigma0=(mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d))*(whereNegative(interpolate(x_hat[flip],g_s.getFunctionSpace()))-0.5)
95 pde=LinearPDE(dom)
96 pde.setSymmetryOn()
97 A =Tensor4(0.,Function(dom))
98 for i in range(d):
99 for j in range(d):
100 A[i,j,j,i] += mu
101 A[i,j,i,j] += mu
102 A[i,i,j,j] += lmbd
103 pde.setValue(A=A,q=mask,Y=-kronecker(Function(dom))[d-1]*g*rho,X=sigma0)
104
105 print("=== solve pde === ")
106 u=pde.getSolution(verbose=True)
107
108 print("=== calculate stress ===")
109 g_s=grad(u)
110 sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d)
111
112 print("=== start rendering ===")
113 # Create a Scene.
114 s = Scene(renderer = JPG_RENDERER, num_viewport = 1, x_size=XSIZE, y_size=YSIZE )
115 dc1 = DataCollector(source = Source.ESCRIPT)
116 dc1.setData(disp=u,sigma=sigma,cfs=sigma[0,1]-0.4*sigma[0,0])
117
118 # Create a Contour.
119 ctr1 = Contour(scene = s, data_collector = dc1, viewport = Viewport.SOUTH_WEST,
120 lut = 'color', cell_to_point = True, outline = True)
121 ctr1.generateContours(contours = 5, lower_range=VIS_RANGE_STRESS[0], upper_range=VIS_RANGE_STRESS[1])
122 ctr1.setScalarRange(lower_range=VIS_RANGE_STRESS[0], upper_range=VIS_RANGE_STRESS[1])
123 ctr1.setOpacity(1.0)
124
125 # Create a Camera.
126 cam1 = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
127 cam1.elevation(angle =-50)
128 cam1.roll(angle = 10)
129 cam1.dolly(6500)
130
131 lg=Legend(scene = s, data_collector = dc1, legend='scalar')
132 lg.setOrientationToVertical()
133 lg.setScalarRange(lower_range=VIS_RANGE_STRESS[0], upper_range=VIS_RANGE_STRESS[1])
134
135 # add some text
136 t2 = Text2D(scene = s, text = "Coulomb failure stress distribution around vertical fault in the Earth crust.")
137 t2.setPosition(LocalPosition(40,YSIZE-50))
138 t2.setColor(color = Color.BLACK)
139 t2.setFontSize(size = 22)
140 t2.setFontToArial()
141 t2.boldOn()
142 t2.shadowOn()
143 # add some text
144 t3 = Text2D(scene = s, text = "Earth Systems Science Computational Center at The University of Queensland.")
145 t3.setPosition(LocalPosition(30,30))
146 t3.setColor(color = Color.BLACK)
147 t3.setFontSize(size = 20)
148 t3.setFontToArial()
149 t3.shadowOn()
150
151 # Render the object.
152 s.render("stress_contour.jpg")
153 print("=== finished rendering ===")

  ViewVC Help
Powered by ViewVC 1.1.26