test/python/slip_stress_old.py


##############################################################################
#
# Copyright (c) 2003-2018 by The University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Apache License, version 2.0
# http://www.apache.org/licenses/LICENSE-2.0
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development 2012-2013 by School of Earth Sciences
# Development from 2014 by Centre for Geoscience Computing (GeoComp)
#
##############################################################################

from __future__ import print_function, division

__copyright__="""Copyright (c) 2003-2018 by The University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Apache License, version 2.0
http://www.apache.org/licenses/LICENSE-2.0"""
__url__="https://launchpad.net/escript-finley"

"""
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"

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


rho=0.
lam_lmbd=1.7e11
lam_mu=1.7e11
g=9.81
fstart =  [50000.0, 40000.0, 10909.09090909091]
fend =  [50000.0, 60000.0,   19090.909090909092]


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)
         # print "p:",inf(p),sup(p)
         # print "u:",inf(u),sup(u)
         self.__pde_u.setValue(y_contact=-p)
         return  self.__pde_u.getSolution()

      def solve_g(self,u,tol=1.e-8):
         dp=Vector(0.,FunctionOnContactZero(self.domain))
         h=FunctionOnContactZero(self.domain).getSize()
         # print jump(u)-self.slip
         dp[0]=(self.slip[0]-jump(u[0]))*lam_mu/h
         dp[1]=(self.slip[1]-jump(u[1]))*lam_mu/h
         dp[2]=(self.slip[2]-jump(u[2]))*lam_mu/h
         return  dp


dom=ReadMesh("meshfault3D.fly",integrationOrder=-1)
prop=SlippingFault(dom)
d=dom.getDim()
x=dom.getX()[0]
# x=dom.getX()[d-1]
mask=whereZero(x-inf(x))*numpy.ones((d,))
x=FunctionOnContactZero(dom).getX()
s=numpy.array([-100000.,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=50,tolerance=0.13,accepted_reduction=1.)
saveVTK("dis.vtu",u=u)
saveVTK("fault.vtu",sigma=p,s=jump(u))
1	ksteube	1809
2	jfenwick	3981	##############################################################################
3	ksteube	1312	#
4	jfenwick	6651	# Copyright (c) 2003-2018 by The University of Queensland
5	jfenwick	3981	# http://www.uq.edu.au
6	ksteube	1312	#
7	ksteube	1809	# Primary Business: Queensland, Australia
8	jfenwick	6112	# Licensed under the Apache License, version 2.0
9			# http://www.apache.org/licenses/LICENSE-2.0
10	ksteube	1312	#
11	jfenwick	3981	# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12	jfenwick	4657	# Development 2012-2013 by School of Earth Sciences
13			# Development from 2014 by Centre for Geoscience Computing (GeoComp)
14	jfenwick	3981	#
15			##############################################################################
16	gross	883
17	sshaw	5706	from __future__ import print_function, division
18
19	jfenwick	6651	__copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20	jfenwick	3981	http://www.uq.edu.au
21	ksteube	1809	Primary Business: Queensland, Australia"""
22	jfenwick	6112	__license__="""Licensed under the Apache License, version 2.0
23			http://www.apache.org/licenses/LICENSE-2.0"""
24	jfenwick	2344	__url__="https://launchpad.net/escript-finley"
25	ksteube	1809
26	gross	883	"""
27			calculation of the stress distribution around a fault from the slip on the fault
28
29	gross	884	e.g. use slip_stress_mesh.py to generate mesh
30
31	jfenwick	2625	:var __author__: name of author
32			:var __copyright__: copyrights
33			:var __license__: licence agreement
34			:var __url__: url entry point on documentation
35			:var __version__: version
36			:var __date__: date of the version
37	gross	883	"""
38
39	gross	884	__author__="Lutz Gross, Louise Kettle"
40	gross	883
41			from esys.escript import *
42			from esys.escript.pdetools import SaddlePointProblem
43			from esys.escript.linearPDEs import LinearPDE
44	gross	884	from esys.finley import ReadMesh
45	caltinay	3346	from esys.weipa import saveVTK
46	gross	883
47
48	gross	884	rho=0.
49	gross	887	lam_lmbd=1.7e11
50			lam_mu=1.7e11
51	gross	883	g=9.81
52	gross	893	fstart = [50000.0, 40000.0, 10909.09090909091]
53	gross	894	fend = [50000.0, 60000.0, 19090.909090909092]
54	gross	883
55	gross	887
56
57	gross	893
58	gross	883	class SlippingFault(SaddlePointProblem):
59			"""
60			simple example of saddle point problem
61			"""
62			def __init__(self,domain):
63			super(SlippingFault, self).__init__(self)
64			self.domain=domain
65			self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
66			self.__pde_u.setSymmetryOn()
67
68			def initialize(self,density=1.,lmbd=1., mu=1., traction=Data(),fixed_u_mask=Data(), slip=0.):
69			d=self.domain.getDim()
70			self.slip=slip
71			A =self.__pde_u.createCoefficientOfGeneralPDE("A")
72			for i in range(self.domain.getDim()):
73			for j in range(self.domain.getDim()):
74			A[i,j,j,i] += mu
75			A[i,j,i,j] += mu
76			A[i,i,j,j] += lmbd
77			self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=-kronecker(Function(self.domain))[d-1]gdensity,y=traction)
78
79			def inner(self,p0,p1):
80	gross	887	return integrate(inner(p0,p1),FunctionOnContactZero(self.domain))
81	gross	883
82			def solve_f(self,u,p,tol=1.e-8):
83			self.__pde_u.setTolerance(tol)
84	gross	893	self.__pde_u.setValue(y_contact=-p)
85			# print "p:",inf(p),sup(p)
86			# print "u:",inf(u),sup(u)
87			self.__pde_u.setValue(y_contact=-p)
88	gross	883	return self.__pde_u.getSolution()
89
90			def solve_g(self,u,tol=1.e-8):
91	gross	893	dp=Vector(0.,FunctionOnContactZero(self.domain))
92			h=FunctionOnContactZero(self.domain).getSize()
93			# print jump(u)-self.slip
94			dp[0]=(self.slip[0]-jump(u[0]))*lam_mu/h
95			dp[1]=(self.slip[1]-jump(u[1]))*lam_mu/h
96			dp[2]=(self.slip[2]-jump(u[2]))*lam_mu/h
97	gross	883	return dp
98
99
100	gross	888	dom=ReadMesh("meshfault3D.fly",integrationOrder=-1)
101	gross	883	prop=SlippingFault(dom)
102			d=dom.getDim()
103	gross	888	x=dom.getX()[0]
104			# x=dom.getX()[d-1]
105	gross	2468	mask=whereZero(x-inf(x))*numpy.ones((d,))
106	gross	887	x=FunctionOnContactZero(dom).getX()
107	gross	2468	s=numpy.array([-100000.,1.,1.])
108	gross	887	for i in range(3):
109			d=fend[i]-fstart[i]
110			if d>0:
111			q=(x[i]-fstart[i])/d
112			s=q(1-q)4*s
113			elif d<0:
114			q=(x[i]-fend[i])/d
115			s=q(1-q)4*s
116	gross	883	u0=Vector(0.,Solution(dom))
117	gross	884	p0=Vector(1.,FunctionOnContactZero(dom))
118			prop.initialize(fixed_u_mask=mask,slip=Data(s,FunctionOnContactZero(dom)), density=rho,lmbd=lam_lmbd, mu=lam_mu)
119	gross	893	u,p=prop.solve(u0,p0,iter_max=50,tolerance=0.13,accepted_reduction=1.)
120	caltinay	2534	saveVTK("dis.vtu",u=u)
121			saveVTK("fault.vtu",sigma=p,s=jump(u))