/[escript]/trunk/finley/test/python/slip_stress_old.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 883 - (hide annotations)
Mon Oct 30 06:12:54 2006 UTC (13 years, 6 months ago) by gross
Original Path: trunk/finley/test/python/slip_stress.py
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     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