/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 4009 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2018 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Apache License, version 2.0
9 # http://www.apache.org/licenses/LICENSE-2.0
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 from __future__ import print_function, division
18
19 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Apache License, version 2.0
23 http://www.apache.org/licenses/LICENSE-2.0"""
24 __url__="https://launchpad.net/escript-finley"
25
26 """
27 calculation of the stress distribution around a fault from the slip on the fault
28
29 e.g. use slip_stress_mesh.py to generate mesh
30
31 :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 """
38
39 __author__="Lutz Gross, Louise Kettle"
40
41 from esys.escript import *
42 from esys.escript.pdetools import SaddlePointProblem
43 from esys.escript.linearPDEs import LinearPDE
44 from esys.finley import ReadMesh
45 from esys.weipa import saveVTK
46
47
48 rho=0.
49 lam_lmbd=1.7e11
50 lam_mu=1.7e11
51 g=9.81
52 fstart = [50000.0, 40000.0, 10909.09090909091]
53 fend = [50000.0, 60000.0, 19090.909090909092]
54
55
56
57
58 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]*g*density,y=traction)
78
79 def inner(self,p0,p1):
80 return integrate(inner(p0,p1),FunctionOnContactZero(self.domain))
81
82 def solve_f(self,u,p,tol=1.e-8):
83 self.__pde_u.setTolerance(tol)
84 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 return self.__pde_u.getSolution()
89
90 def solve_g(self,u,tol=1.e-8):
91 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 return dp
98
99
100 dom=ReadMesh("meshfault3D.fly",integrationOrder=-1)
101 prop=SlippingFault(dom)
102 d=dom.getDim()
103 x=dom.getX()[0]
104 # x=dom.getX()[d-1]
105 mask=whereZero(x-inf(x))*numpy.ones((d,))
106 x=FunctionOnContactZero(dom).getX()
107 s=numpy.array([-100000.,1.,1.])
108 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 u0=Vector(0.,Solution(dom))
117 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 u,p=prop.solve(u0,p0,iter_max=50,tolerance=0.13,accepted_reduction=1.)
120 saveVTK("dis.vtu",u=u)
121 saveVTK("fault.vtu",sigma=p,s=jump(u))

  ViewVC Help
Powered by ViewVC 1.1.26