/[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 883 - (show annotations)
Mon Oct 30 06:12:54 2006 UTC (13 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 2610 byte(s)
new example added
1 # $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 http://www.access.edu.au
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="http://www.iservo.edu.au/esys"
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))

  ViewVC Help
Powered by ViewVC 1.1.26