test/python/slip_stress.py

# $Id$
"""
calculation of the stress distribution around a fault from the slip on the fault

@var __author__: name of author
@var __copyright__: copyrights
@var __license__: licence agreement
@var __url__: url entry point on documentation
@var __version__: version
@var __date__: date of the version
"""

__author__="Lutz Gross, Louise Kettle"
__copyright__="""  Copyright (c) 2006, 2007 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""
__url__="http://www.iservo.edu.au/esys"
__version__="$Revision$"
__date__="$Date$"

from esys.escript import *
from esys.escript.pdetools import SaddlePointProblem
from esys.escript.linearPDEs import LinearPDE
from esys.finley import Brick
from esys.pyvisi import Scene, DataCollector, Contour, Camera, Velocity, Text2D, LocalPosition
from esys.pyvisi.constant import *
import os

JPG_RENDERER = Renderer.OFFLINE_JPG
ne = 20          # number of elements in spatial direction

def lockToGrid(x,l,n):
   h=l/n
   return int(x/h+0.5)*h


width  = 100000.
height =  30000.

ne_w=int((ne/height)*width+0.5)
rho=0.
lmbd=1.7e11
mu=1.7e11
g=9.81

# fault location:

fstart =  [lockToGrid(50000.0,width,ne_w), lockToGrid(40000.0,width,ne_w), lockToGrid(8000.,height,ne)]
fend =  [lockToGrid(50000.0,width,ne_w), lockToGrid(60000.0,width,ne_w), lockToGrid(20000.,height,ne)]
s=[0.,1.,0.]

print "=== generate mesh over %s x %s x %s ==="%(width,width,height)
dom=Brick(l0=width, l1=width, l2=height, n0=ne_w, n1=ne_w, n2=ne)
print "  total number of elements = ",ne*ne_w*ne_w

print "=== prepare PDE coefficients ==="
# fixed displacements:
d=dom.getDim()
x=dom.getX()
mask=whereZero(x[d-1]-inf(x[d-1]))*numarray.ones((d,))

# map the fault into the unit square:
x_hat=x-fstart
p=[0.,0.,0.]
q=1.
for i in range(d):
     p[i]=fend[i]-fstart[i]
     if abs(p[i])>0: 
         x_hat[i]/=p[i]
         q=whereNonNegative(x_hat[i])*whereNonPositive(x_hat[i]-1.)*x_hat[i]*(1.-x_hat[i])*4.*q
     else:
         flip=i
         q=whereZero(x_hat[i],1.e-3*Lsup(x_hat[i]))*q
g_s=grad(q*s)
sigma0=(mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d))*(whereNegative(interpolate(x_hat[flip],g_s.getFunctionSpace()))-0.5)
pde=LinearPDE(dom)
pde.setSymmetryOn()
A =Tensor4(0.,Function(dom))
for i in range(d):
   for j in range(d):
      A[i,j,j,i] += mu
      A[i,j,i,j] += mu
      A[i,i,j,j] += lmbd
pde.setValue(A=A,q=mask,Y=-kronecker(Function(dom))[d-1]*g*rho,X=sigma0)

print "=== solve pde === "
u=pde.getSolution(verbose=True)

print "=== calculate stress ==="
g_s=grad(u)
sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d)

print "=== start rendering ==="
# Create a Scene.
s = Scene(renderer = JPG_RENDERER, num_viewport = 1, x_size = 800, y_size = 800)
dc1 = DataCollector(source = Source.ESCRIPT)
dc1.setData(disp=u,sigma=sigma,cfs=sigma[0,1]-0.4*sigma[0,0])

# Create a Contour.
ctr1 = Contour(scene = s, data_collector = dc1, viewport = Viewport.SOUTH_WEST,
        lut = Lut.COLOR, cell_to_point = True, outline = True)
ctr1.generateContours(contours = 10)
ctr1.setOpacity(0.5)

# Create velocity
vopc1 = Velocity(scene = s, data_collector = dc1,
        color_mode = ColorMode.VECTOR,
        arrow = Arrow.THREE_D, lut = Lut.COLOR, cell_to_point = False,
        outline = False) 
vopc1.setScaleFactor(scale_factor = 10000.)
vopc1.setRatio(2)
vopc1.randomOn()

# Create a Camera.
cam1 = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
cam1.elevation(angle = -50)

# add some text
t2 = Text2D(scene = s, text = "CFS and displacement around vertical fault")
t2.setPosition(LocalPosition(200,30))
t2.setColor(color = Color.BLACK)
t2.setFontSize(size = 20)
t2.setFontToArial()
t2.shadowOn()

# Render the object.
s.render("stress_contour.jpg")
print "=== finished rendering ==="
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 = ",nene_wne_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-3Lsup(x_hat[i]))q
76	g_s=grad(q*s)
77	sigma0=(musymmetric(g_s)+lmbdtrace(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]grho,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=musymmetric(g_s)+lmbdtrace(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 ==="