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

Revision 895 - (hide annotations)
Thu Nov 9 05:45:55 2006 UTC (14 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 3436 byte(s)
```moved to make plave for better version
```
 1 gross 883 # \$Id\$ 2 3 """ 4 calculation of the stress distribution around a fault from the slip on the fault 5 6 gross 884 e.g. use slip_stress_mesh.py to generate mesh 7 8 gross 883 @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 gross 884 __author__="Lutz Gross, Louise Kettle" 17 gross 883 __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 gross 884 from esys.finley import ReadMesh 30 gross 883 31 32 gross 884 rho=0. 33 gross 887 lam_lmbd=1.7e11 34 lam_mu=1.7e11 35 gross 883 g=9.81 36 gross 893 fstart = [50000.0, 40000.0, 10909.09090909091] 37 gross 894 fend = [50000.0, 60000.0, 19090.909090909092] 38 gross 883 39 gross 887 40 41 gross 893 42 gross 883 class SlippingFault(SaddlePointProblem): 43 """ 44 simple example of saddle point problem 45 """ 46 def __init__(self,domain): 47 super(SlippingFault, self).__init__(self) 48 self.domain=domain 49 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 50 self.__pde_u.setSymmetryOn() 51 52 def initialize(self,density=1.,lmbd=1., mu=1., traction=Data(),fixed_u_mask=Data(), slip=0.): 53 d=self.domain.getDim() 54 self.slip=slip 55 A =self.__pde_u.createCoefficientOfGeneralPDE("A") 56 for i in range(self.domain.getDim()): 57 for j in range(self.domain.getDim()): 58 A[i,j,j,i] += mu 59 A[i,j,i,j] += mu 60 A[i,i,j,j] += lmbd 61 self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=-kronecker(Function(self.domain))[d-1]*g*density,y=traction) 62 63 def inner(self,p0,p1): 64 gross 887 return integrate(inner(p0,p1),FunctionOnContactZero(self.domain)) 65 gross 883 66 def solve_f(self,u,p,tol=1.e-8): 67 self.__pde_u.setTolerance(tol) 68 gross 893 self.__pde_u.setValue(y_contact=-p) 69 # print "p:",inf(p),sup(p) 70 # print "u:",inf(u),sup(u) 71 self.__pde_u.setValue(y_contact=-p) 72 gross 883 return self.__pde_u.getSolution() 73 74 def solve_g(self,u,tol=1.e-8): 75 gross 893 dp=Vector(0.,FunctionOnContactZero(self.domain)) 76 h=FunctionOnContactZero(self.domain).getSize() 77 # print jump(u)-self.slip 78 dp[0]=(self.slip[0]-jump(u[0]))*lam_mu/h 79 dp[1]=(self.slip[1]-jump(u[1]))*lam_mu/h 80 dp[2]=(self.slip[2]-jump(u[2]))*lam_mu/h 81 gross 883 return dp 82 83 84 gross 888 dom=ReadMesh("meshfault3D.fly",integrationOrder=-1) 85 gross 883 prop=SlippingFault(dom) 86 d=dom.getDim() 87 gross 888 x=dom.getX()[0] 88 # x=dom.getX()[d-1] 89 gross 883 mask=whereZero(x-inf(x))*numarray.ones((d,)) 90 gross 887 x=FunctionOnContactZero(dom).getX() 91 gross 893 s=numarray.array([-100000.,1.,1.]) 92 gross 887 for i in range(3): 93 d=fend[i]-fstart[i] 94 if d>0: 95 q=(x[i]-fstart[i])/d 96 s=q*(1-q)*4*s 97 elif d<0: 98 q=(x[i]-fend[i])/d 99 s=q*(1-q)*4*s 100 gross 883 u0=Vector(0.,Solution(dom)) 101 gross 884 p0=Vector(1.,FunctionOnContactZero(dom)) 102 prop.initialize(fixed_u_mask=mask,slip=Data(s,FunctionOnContactZero(dom)), density=rho,lmbd=lam_lmbd, mu=lam_mu) 103 gross 893 u,p=prop.solve(u0,p0,iter_max=50,tolerance=0.13,accepted_reduction=1.) 104 gross 883 saveVTK("dis.xml",u=u) 105 saveVTK("fault.xml",sigma=p,s=jump(u))

 ViewVC Help Powered by ViewVC 1.1.26