/[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 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years, 5 months ago) by ksteube
File MIME type: text/x-python
File size: 3856 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

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

  ViewVC Help
Powered by ViewVC 1.1.26