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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 883 by gross, Mon Oct 30 06:12:54 2006 UTC revision 893 by gross, Wed Nov 8 08:20:19 2006 UTC
# Line 3  Line 3 
3  """  """
4  calculation of the stress distribution around a fault from the slip on the fault  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  @var __author__: name of author
9  @var __copyright__: copyrights  @var __copyright__: copyrights
10  @var __license__: licence agreement  @var __license__: licence agreement
# Line 11  calculation of the stress distribution a Line 13  calculation of the stress distribution a
13  @var __date__: date of the version  @var __date__: date of the version
14  """  """
15    
16  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, Louise Kettle"
17  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
18                      http://www.access.edu.au                      http://www.access.edu.au
19                  Primary Business: Queensland, Australia"""                  Primary Business: Queensland, Australia"""
# Line 24  __date__="$Date$" Line 26  __date__="$Date$"
26  from esys.escript import *  from esys.escript import *
27  from esys.escript.pdetools import SaddlePointProblem  from esys.escript.pdetools import SaddlePointProblem
28  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
29  from esys.finley import Rectangle  from esys.finley import ReadMesh
30    
31    
32  rho=1.  rho=0.
33  lam_lmbd=1.  lam_lmbd=1.7e11
34  lam_mu=1.  lam_mu=1.7e11
35  g=9.81  g=9.81
36    fstart =  [50000.0, 40000.0, 10909.09090909091]
37    fend =  [50000.0, 60000.0, 19090.909090909092]
38    
39    
40    
41    
42  class SlippingFault(SaddlePointProblem):  class SlippingFault(SaddlePointProblem):
43        """        """
# Line 54  class SlippingFault(SaddlePointProblem): Line 61  class SlippingFault(SaddlePointProblem):
61           self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=-kronecker(Function(self.domain))[d-1]*g*density,y=traction)           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):        def inner(self,p0,p1):
64           return integrate(p0*p1,FunctionOnContactZero(self.__pde_p.getDomain()))           return integrate(inner(p0,p1),FunctionOnContactZero(self.domain))
65    
66        def solve_f(self,u,p,tol=1.e-8):        def solve_f(self,u,p,tol=1.e-8):
67           self.__pde_u.setTolerance(tol)           self.__pde_u.setTolerance(tol)
68           self.__pde_u.setValue(y_contact=p)           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           return  self.__pde_u.getSolution()           return  self.__pde_u.getSolution()
73    
74        def solve_g(self,u,tol=1.e-8):        def solve_g(self,u,tol=1.e-8):
75           dp=self.slip-jump(u)           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           return  dp           return  dp
82    
83    
84  s=numarray.array([0.,1.,1.])  dom=ReadMesh("meshfault3D.fly",integrationOrder=-1)
 NE=3  
 dom=Rectangle(NE,NE,order=2)  
85  prop=SlippingFault(dom)  prop=SlippingFault(dom)
86  d=dom.getDim()  d=dom.getDim()
87  x=dom.getX()[d-1]  x=dom.getX()[0]
88    # x=dom.getX()[d-1]
89  mask=whereZero(x-inf(x))*numarray.ones((d,))  mask=whereZero(x-inf(x))*numarray.ones((d,))
90    x=FunctionOnContactZero(dom).getX()
91    s=numarray.array([-100000.,1.,1.])
92    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  u0=Vector(0.,Solution(dom))  u0=Vector(0.,Solution(dom))
101  p0=Vector(0.,FunctionOnContactZero(dom))  p0=Vector(1.,FunctionOnContactZero(dom))
102  prop.initialize(fixed_u_mask=mask,slip=s, density=rho,lmbd=lam_lmbd, mu=lam_mu)  prop.initialize(fixed_u_mask=mask,slip=Data(s,FunctionOnContactZero(dom)), density=rho,lmbd=lam_lmbd, mu=lam_mu)
103  u,p=prop.solve(u0,p0,relaxation=1.,iter_max=50,tolerance=0.01)  u,p=prop.solve(u0,p0,iter_max=50,tolerance=0.13,accepted_reduction=1.)
104  saveVTK("dis.xml",u=u)  saveVTK("dis.xml",u=u)
105  saveVTK("fault.xml",sigma=p,s=jump(u))  saveVTK("fault.xml",sigma=p,s=jump(u))

Legend:
Removed from v.883  
changed lines
  Added in v.893

  ViewVC Help
Powered by ViewVC 1.1.26