test/python/slip_stress.py

# $Id$

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

e.g. use slip_stress_mesh.py to generate mesh

@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 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

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


ne = 15
width  = 100000.
height =  30000.

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

# fault:

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.,0.,1.]

fstart =  [lockToGrid(30000.0,width,ne_w), lockToGrid(30000.0,width,ne_w), lockToGrid(20000.,height,ne)]
fend =  [lockToGrid(70000.0,width,ne_w), lockToGrid(70000.0,width,ne_w), lockToGrid(20000.,height,ne)]
s=[1.,0.1,0.]

dom=Brick(l0=width, l1=width, l2=height, n0=ne_w, n1=ne_w, n2=ne)

print "fstart:",fstart
print "fend:",fend
print "slip ",s

# 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]
     print i,p[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)
u=pde.getSolution(verbose=True)
g_s=grad(u)
saveVTK("dis.xml",disp=u,sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d))
1	# $Id$
2
3	"""
4	calculation of the stress distribution around a fault from the slip on the fault
5
6	e.g. use slip_stress_mesh.py to generate mesh
7
8	@var __author__: name of author
9	@var __copyright__: copyrights
10	@var __license__: licence agreement
11	@var __url__: url entry point on documentation
12	@var __version__: version
13	@var __date__: date of the version
14	"""
15
16	__author__="Lutz Gross, Louise Kettle"
17	__copyright__=""" Copyright (c) 2006 by ACcESS MNRF
18	http://www.access.edu.au
19	Primary Business: Queensland, Australia"""
20	__license__="""Licensed under the Open Software License version 3.0
21	http://www.opensource.org/licenses/osl-3.0.php"""
22	__url__="http://www.iservo.edu.au/esys"
23	__version__="$Revision$"
24	__date__="$Date$"
25
26	from esys.escript import *
27	from esys.escript.pdetools import SaddlePointProblem
28	from esys.escript.linearPDEs import LinearPDE
29	from esys.finley import Brick
30
31	def lockToGrid(x,l,n):
32	h=l/n
33	print x,l,n,"->",h
34	return int(x/h+0.5)*h
35
36
37	ne = 15
38	width = 100000.
39	height = 30000.
40
41	ne_w=int((ne/height)*width+0.5)
42	rho=0.
43	lmbd=1.7e11
44	mu=1.7e11
45	g=9.81
46	print "total ne = ",nene_wne_w
47
48	# fault:
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.,0.,1.]
53
54	fstart = [lockToGrid(30000.0,width,ne_w), lockToGrid(30000.0,width,ne_w), lockToGrid(20000.,height,ne)]
55	fend = [lockToGrid(70000.0,width,ne_w), lockToGrid(70000.0,width,ne_w), lockToGrid(20000.,height,ne)]
56	s=[1.,0.1,0.]
57
58	dom=Brick(l0=width, l1=width, l2=height, n0=ne_w, n1=ne_w, n2=ne)
59
60	print "fstart:",fstart
61	print "fend:",fend
62	print "slip ",s
63
64	# fixed displacements:
65	d=dom.getDim()
66	x=dom.getX()
67	mask=whereZero(x[d-1]-inf(x[d-1]))*numarray.ones((d,))
68
69	# map the fault into the unit square:
70	x_hat=x-fstart
71	p=[0.,0.,0.]
72	q=1.
73	for i in range(d):
74	p[i]=fend[i]-fstart[i]
75	print i,p[i]
76	if abs(p[i])>0:
77	x_hat[i]/=p[i]
78	q=whereNonNegative(x_hat[i])whereNonPositive(x_hat[i]-1.)x_hat[i](1.-x_hat[i])4.*q
79	else:
80	flip=i
81	q=whereZero(x_hat[i],1.e-3Lsup(x_hat[i]))q
82	g_s=grad(q*s)
83	sigma0=(musymmetric(g_s)+lmbdtrace(g_s)kronecker(d))(whereNegative(interpolate(x_hat[flip],g_s.getFunctionSpace()))-0.5)
84	pde=LinearPDE(dom)
85	pde.setSymmetryOn()
86	A =Tensor4(0.,Function(dom))
87	for i in range(d):
88	for j in range(d):
89	A[i,j,j,i] += mu
90	A[i,j,i,j] += mu
91	A[i,i,j,j] += lmbd
92	pde.setValue(A=A,q=mask,Y=-kronecker(Function(dom))[d-1]grho,X=-sigma0)
93	u=pde.getSolution(verbose=True)
94	g_s=grad(u)
95	saveVTK("dis.xml",disp=u,sigma=musymmetric(g_s)+lmbdtrace(g_s)*kronecker(d))