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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2549 - (show annotations)
Mon Jul 20 06:43:47 2009 UTC (10 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 3811 byte(s)
Remainder of copyright date fixes
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
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__="https://launchpad.net/escript-finley"
21
22 """
23 calculation of the stress distribution around a fault from the slip on the fault
24
25 e.g. use slip_stress_mesh.py to generate mesh
26
27 @var __author__: name of author
28 @var __copyright__: copyrights
29 @var __license__: licence agreement
30 @var __url__: url entry point on documentation
31 @var __version__: version
32 @var __date__: date of the version
33 """
34
35 __author__="Lutz Gross, Louise Kettle"
36
37 from esys.escript import *
38 from esys.escript.pdetools import SaddlePointProblem
39 from esys.escript.linearPDEs import LinearPDE
40 from esys.finley import ReadMesh
41
42
43 rho=0.
44 lam_lmbd=1.7e11
45 lam_mu=1.7e11
46 g=9.81
47 fstart = [50000.0, 40000.0, 10909.09090909091]
48 fend = [50000.0, 60000.0, 19090.909090909092]
49
50
51
52
53 class SlippingFault(SaddlePointProblem):
54 """
55 simple example of saddle point problem
56 """
57 def __init__(self,domain):
58 super(SlippingFault, self).__init__(self)
59 self.domain=domain
60 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
61 self.__pde_u.setSymmetryOn()
62
63 def initialize(self,density=1.,lmbd=1., mu=1., traction=Data(),fixed_u_mask=Data(), slip=0.):
64 d=self.domain.getDim()
65 self.slip=slip
66 A =self.__pde_u.createCoefficientOfGeneralPDE("A")
67 for i in range(self.domain.getDim()):
68 for j in range(self.domain.getDim()):
69 A[i,j,j,i] += mu
70 A[i,j,i,j] += mu
71 A[i,i,j,j] += lmbd
72 self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=-kronecker(Function(self.domain))[d-1]*g*density,y=traction)
73
74 def inner(self,p0,p1):
75 return integrate(inner(p0,p1),FunctionOnContactZero(self.domain))
76
77 def solve_f(self,u,p,tol=1.e-8):
78 self.__pde_u.setTolerance(tol)
79 self.__pde_u.setValue(y_contact=-p)
80 # print "p:",inf(p),sup(p)
81 # print "u:",inf(u),sup(u)
82 self.__pde_u.setValue(y_contact=-p)
83 return self.__pde_u.getSolution()
84
85 def solve_g(self,u,tol=1.e-8):
86 dp=Vector(0.,FunctionOnContactZero(self.domain))
87 h=FunctionOnContactZero(self.domain).getSize()
88 # print jump(u)-self.slip
89 dp[0]=(self.slip[0]-jump(u[0]))*lam_mu/h
90 dp[1]=(self.slip[1]-jump(u[1]))*lam_mu/h
91 dp[2]=(self.slip[2]-jump(u[2]))*lam_mu/h
92 return dp
93
94
95 dom=ReadMesh("meshfault3D.fly",integrationOrder=-1)
96 prop=SlippingFault(dom)
97 d=dom.getDim()
98 x=dom.getX()[0]
99 # x=dom.getX()[d-1]
100 mask=whereZero(x-inf(x))*numpy.ones((d,))
101 x=FunctionOnContactZero(dom).getX()
102 s=numpy.array([-100000.,1.,1.])
103 for i in range(3):
104 d=fend[i]-fstart[i]
105 if d>0:
106 q=(x[i]-fstart[i])/d
107 s=q*(1-q)*4*s
108 elif d<0:
109 q=(x[i]-fend[i])/d
110 s=q*(1-q)*4*s
111 u0=Vector(0.,Solution(dom))
112 p0=Vector(1.,FunctionOnContactZero(dom))
113 prop.initialize(fixed_u_mask=mask,slip=Data(s,FunctionOnContactZero(dom)), density=rho,lmbd=lam_lmbd, mu=lam_mu)
114 u,p=prop.solve(u0,p0,iter_max=50,tolerance=0.13,accepted_reduction=1.)
115 saveVTK("dis.vtu",u=u)
116 saveVTK("fault.vtu",sigma=p,s=jump(u))

  ViewVC Help
Powered by ViewVC 1.1.26