test/python/slip_stress.py

# $Id$

"""
calculation of the stress distribution around a fault from the slip on the fault

e.g. use slip_stress_mesh.py to generate mesh

@var __author__: name of author
@var __copyright__: copyrights
@var __license__: licence agreement
@var __url__: url entry point on documentation
@var __version__: version
@var __date__: date of the version
"""

__author__="Lutz Gross, Louise Kettle"
__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""
__url__="http://www.iservo.edu.au/esys"
__version__="$Revision$"
__date__="$Date$"

from esys.escript import *
from esys.escript.pdetools import SaddlePointProblem
from esys.escript.linearPDEs import LinearPDE
from esys.finley import ReadMesh


rho=0.
lam_lmbd=1.7e11
lam_mu=1.7e11
g=9.81
fstart =  [50000.0, 44444.444444444445, 11666.666666666666]
fend =  [50000.0, 50000.0, -13333.333333333336]
fstart =  [50000.0, 44444.444444444445, 11666.666666666666]
fend =  [50000.0, 55555.555555555555, -11666.666666666668]


class SlippingFault(SaddlePointProblem):
      """
      simple example of saddle point problem
      """
      def __init__(self,domain):
         super(SlippingFault, self).__init__(self)
         self.domain=domain
         self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
         self.__pde_u.setSymmetryOn()

      def initialize(self,density=1.,lmbd=1., mu=1., traction=Data(),fixed_u_mask=Data(), slip=0.):
         d=self.domain.getDim()
         self.slip=slip
         A =self.__pde_u.createCoefficientOfGeneralPDE("A")
         for i in range(self.domain.getDim()):
           for j in range(self.domain.getDim()):
             A[i,j,j,i] += mu
             A[i,j,i,j] += mu
             A[i,i,j,j] += lmbd
         self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=-kronecker(Function(self.domain))[d-1]*g*density,y=traction)

      def inner(self,p0,p1):
         return integrate(inner(p0,p1),FunctionOnContactZero(self.domain))

      def solve_f(self,u,p,tol=1.e-8):
         self.__pde_u.setTolerance(tol)
         self.__pde_u.setValue(y_contact=p)
         return  self.__pde_u.getSolution()

      def solve_g(self,u,tol=1.e-8):
         dp=(self.slip-jump(u))*lam_lmbd/FunctionOnContactZero(self.domain).getX()
         return  dp


dom=ReadMesh("meshfault3D.fly")
prop=SlippingFault(dom)
d=dom.getDim()
x=dom.getX()[d-1]
mask=whereZero(x-inf(x))*numarray.ones((d,))
s=numarray.array([0.,1.,1.])
x=FunctionOnContactZero(dom).getX()
s=numarray.array([0.,1.,1.])
for i in range(3):
     d=fend[i]-fstart[i]
     if d>0:
         q=(x[i]-fstart[i])/d
         s=q*(1-q)*4*s
     elif d<0:
         q=(x[i]-fend[i])/d
         s=q*(1-q)*4*s
u0=Vector(0.,Solution(dom))
p0=Vector(1.,FunctionOnContactZero(dom))
prop.initialize(fixed_u_mask=mask,slip=Data(s,FunctionOnContactZero(dom)), density=rho,lmbd=lam_lmbd, mu=lam_mu)
u,p=prop.solve(u0,p0,iter_max=100,tolerance=0.5,accepted_reduction=1.1 )
saveVTK("dis.xml",u=u)
saveVTK("fault.xml",sigma=p,s=jump(u))
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]gdensity,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))