/[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 892 - (show annotations)
Wed Nov 8 06:14:29 2006 UTC (13 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 3389 byte(s)


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 ReadMesh
30
31
32 rho=0.
33 lam_lmbd=1.7e11
34 lam_mu=1.7e11
35 g=9.81
36 fstart = [50000.0, 44444.444444444445, 11666.666666666666]
37 fend = [50000.0, 50000.0, -13333.333333333336]
38 fstart = [50000.0, 44444.444444444445, 11666.666666666666]
39 fend = [50000.0, 55555.555555555555, -11666.666666666668]
40
41
42
43 class SlippingFault(SaddlePointProblem):
44 """
45 simple example of saddle point problem
46 """
47 def __init__(self,domain):
48 super(SlippingFault, self).__init__(self)
49 self.domain=domain
50 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
51 self.__pde_u.setSymmetryOn()
52
53 def initialize(self,density=1.,lmbd=1., mu=1., traction=Data(),fixed_u_mask=Data(), slip=0.):
54 d=self.domain.getDim()
55 self.slip=slip
56 A =self.__pde_u.createCoefficientOfGeneralPDE("A")
57 for i in range(self.domain.getDim()):
58 for j in range(self.domain.getDim()):
59 A[i,j,j,i] += mu
60 A[i,j,i,j] += mu
61 A[i,i,j,j] += lmbd
62 self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=-kronecker(Function(self.domain))[d-1]*g*density,y=traction)
63
64 def inner(self,p0,p1):
65 return integrate(inner(p0,p1),FunctionOnContactZero(self.domain))
66
67 def solve_f(self,u,p,tol=1.e-8):
68 self.__pde_u.setTolerance(tol)
69 self.__pde_u.setValue(y_contact=p)
70 print "p:",inf(p),sup(p)
71 print "u:",inf(u),sup(u)
72 self.__pde_u.setValue(y_contact=p)
73 return self.__pde_u.getSolution()
74
75 def solve_g(self,u,tol=1.e-8):
76 dp=-(self.slip-jump(u))*lam_lmbd/FunctionOnContactZero(self.domain).getSize()
77 print dp
78 return dp
79
80
81 dom=ReadMesh("meshfault3D.fly",integrationOrder=-1)
82 prop=SlippingFault(dom)
83 d=dom.getDim()
84 x=dom.getX()[0]
85 # x=dom.getX()[d-1]
86 mask=whereZero(x-inf(x))*numarray.ones((d,))
87 s=numarray.array([0.,1.,1.])
88 x=FunctionOnContactZero(dom).getX()
89 s=numarray.array([0.,1.,1.])
90 for i in range(3):
91 d=fend[i]-fstart[i]
92 if d>0:
93 q=(x[i]-fstart[i])/d
94 s=q*(1-q)*4*s
95 elif d<0:
96 q=(x[i]-fend[i])/d
97 s=q*(1-q)*4*s
98 u0=Vector(0.,Solution(dom))
99 p0=Vector(1.,FunctionOnContactZero(dom))
100 prop.initialize(fixed_u_mask=mask,slip=Data(s,FunctionOnContactZero(dom)), density=rho,lmbd=lam_lmbd, mu=lam_mu)
101 u,p=prop.solve(u0,p0,iter_max=50,tolerance=0.1,accepted_reduction=0.99)
102 saveVTK("dis.xml",u=u)
103 saveVTK("fault.xml",sigma=p,s=jump(u))

  ViewVC Help
Powered by ViewVC 1.1.26