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

Revision 897 - (hide annotations)
Thu Nov 9 07:37:40 2006 UTC (16 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 2755 byte(s)
```cfs added
```
 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 19 Primary Business: Queensland, Australia""" 20 __license__="""Licensed under the Open Software License version 3.0 21 22 __url__= 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 gross 897 s=[0.,1.,0.] 53 gross 896 54 gross 897 # 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.,0.] 57 gross 896 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 gross 897 pde.setValue(A=A,q=mask,Y=-kronecker(Function(dom))[d-1]*g*rho,X=sigma0) 93 gross 896 u=pde.getSolution(verbose=True) 94 g_s=grad(u) 95 gross 897 sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d) 96 saveVTK("dis.xml",disp=u,sigma=sigma,cfs=sigma[0,1]-0.4*sigma[0,0])