/[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 887 - (hide annotations)
Thu Nov 2 07:17:07 2006 UTC (13 years, 7 months ago) by gross
Original Path: trunk/finley/test/python/slip_stress.py
File MIME type: text/x-python
File size: 3218 byte(s)
some more work on the slip distribution
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     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 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 887 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 gross 883
41 gross 887
42    
43 gross 883 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 gross 887 return integrate(inner(p0,p1),FunctionOnContactZero(self.domain))
66 gross 883
67     def solve_f(self,u,p,tol=1.e-8):
68     self.__pde_u.setTolerance(tol)
69 gross 887 self.__pde_u.setValue(y_contact=p)
70 gross 883 return self.__pde_u.getSolution()
71    
72     def solve_g(self,u,tol=1.e-8):
73 gross 887 dp=(self.slip-jump(u))*lam_lmbd/FunctionOnContactZero(self.domain).getX()
74 gross 883 return dp
75    
76    
77 gross 884 dom=ReadMesh("meshfault3D.fly")
78 gross 883 prop=SlippingFault(dom)
79     d=dom.getDim()
80     x=dom.getX()[d-1]
81     mask=whereZero(x-inf(x))*numarray.ones((d,))
82 gross 887 s=numarray.array([0.,1.,1.])
83     x=FunctionOnContactZero(dom).getX()
84     s=numarray.array([0.,1.,1.])
85     for i in range(3):
86     d=fend[i]-fstart[i]
87     if d>0:
88     q=(x[i]-fstart[i])/d
89     s=q*(1-q)*4*s
90     elif d<0:
91     q=(x[i]-fend[i])/d
92     s=q*(1-q)*4*s
93 gross 883 u0=Vector(0.,Solution(dom))
94 gross 884 p0=Vector(1.,FunctionOnContactZero(dom))
95     prop.initialize(fixed_u_mask=mask,slip=Data(s,FunctionOnContactZero(dom)), density=rho,lmbd=lam_lmbd, mu=lam_mu)
96 gross 887 u,p=prop.solve(u0,p0,iter_max=100,tolerance=0.5,accepted_reduction=1.1 )
97 gross 883 saveVTK("dis.xml",u=u)
98     saveVTK("fault.xml",sigma=p,s=jump(u))

  ViewVC Help
Powered by ViewVC 1.1.26