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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 3 months 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 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]*g*density,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))

  ViewVC Help
Powered by ViewVC 1.1.26