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

Contents of /branches/symbolic_from_3470/dudley/test/python/slip_stress.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3789 - (show annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 4738 byte(s)
Fast forward to latest trunk revision 3788.

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

  ViewVC Help
Powered by ViewVC 1.1.26