/[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 896 - (show annotations)
Thu Nov 9 05:46:31 2006 UTC (12 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 2709 byte(s)
a new version
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 = ",ne*ne_w*ne_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-3*Lsup(x_hat[i]))*q
82 g_s=grad(q*s)
83 sigma0=(mu*symmetric(g_s)+lmbd*trace(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]*g*rho,X=-sigma0)
93 u=pde.getSolution(verbose=True)
94 g_s=grad(u)
95 saveVTK("dis.xml",disp=u,sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d))

  ViewVC Help
Powered by ViewVC 1.1.26