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

Annotation of /trunk/finley/test/python/slip_stress.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 896 - (hide annotations)
Thu Nov 9 05:46:31 2006 UTC (13 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 2709 byte(s)
a new version
1 gross 896 # $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