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

Revision 883 - (hide annotations)
Mon Oct 30 06:12:54 2006 UTC (15 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 2610 byte(s)
```new example added
```
 1 gross 883 # \$Id\$ 2 3 """ 4 calculation of the stress distribution around a fault from the slip on the fault 5 6 @var __author__: name of author 7 @var __copyright__: copyrights 8 @var __license__: licence agreement 9 @var __url__: url entry point on documentation 10 @var __version__: version 11 @var __date__: date of the version 12 """ 13 14 __author__="Lutz Gross, l.gross@uq.edu.au" 15 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF 16 17 Primary Business: Queensland, Australia""" 18 __license__="""Licensed under the Open Software License version 3.0 19 20 __url__= 21 __version__="\$Revision\$" 22 __date__="\$Date\$" 23 24 from esys.escript import * 25 from esys.escript.pdetools import SaddlePointProblem 26 from esys.escript.linearPDEs import LinearPDE 27 from esys.finley import Rectangle 28 29 30 rho=1. 31 lam_lmbd=1. 32 lam_mu=1. 33 g=9.81 34 35 class SlippingFault(SaddlePointProblem): 36 """ 37 simple example of saddle point problem 38 """ 39 def __init__(self,domain): 40 super(SlippingFault, self).__init__(self) 41 self.domain=domain 42 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 43 self.__pde_u.setSymmetryOn() 44 45 def initialize(self,density=1.,lmbd=1., mu=1., traction=Data(),fixed_u_mask=Data(), slip=0.): 46 d=self.domain.getDim() 47 self.slip=slip 48 A =self.__pde_u.createCoefficientOfGeneralPDE("A") 49 for i in range(self.domain.getDim()): 50 for j in range(self.domain.getDim()): 51 A[i,j,j,i] += mu 52 A[i,j,i,j] += mu 53 A[i,i,j,j] += lmbd 54 self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=-kronecker(Function(self.domain))[d-1]*g*density,y=traction) 55 56 def inner(self,p0,p1): 57 return integrate(p0*p1,FunctionOnContactZero(self.__pde_p.getDomain())) 58 59 def solve_f(self,u,p,tol=1.e-8): 60 self.__pde_u.setTolerance(tol) 61 self.__pde_u.setValue(y_contact=p) 62 return self.__pde_u.getSolution() 63 64 def solve_g(self,u,tol=1.e-8): 65 dp=self.slip-jump(u) 66 return dp 67 68 69 s=numarray.array([0.,1.,1.]) 70 NE=3 71 dom=Rectangle(NE,NE,order=2) 72 prop=SlippingFault(dom) 73 d=dom.getDim() 74 x=dom.getX()[d-1] 75 mask=whereZero(x-inf(x))*numarray.ones((d,)) 76 u0=Vector(0.,Solution(dom)) 77 p0=Vector(0.,FunctionOnContactZero(dom)) 78 prop.initialize(fixed_u_mask=mask,slip=s, density=rho,lmbd=lam_lmbd, mu=lam_mu) 79 u,p=prop.solve(u0,p0,relaxation=1.,iter_max=50,tolerance=0.01) 80 saveVTK("dis.xml",u=u) 81 saveVTK("fault.xml",sigma=p,s=jump(u))