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

Contents of /trunk/finley/test/python/slip_stress.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years, 4 months ago) by ksteube
File MIME type: text/x-python
File size: 4328 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 @var __author__: name of author
21 @var __copyright__: copyrights
22 @var __license__: licence agreement
23 @var __url__: url entry point on documentation
24 @var __version__: version
25 @var __date__: date of the version
26 """
27
28 __author__="Lutz Gross, Louise Kettle"
29 __copyright__=""" Copyright (c) 2006, 2007 by ACcESS MNRF
30 http://www.access.edu.au
31 Primary Business: Queensland, Australia"""
32 __license__="""Licensed under the Open Software License version 3.0
33 http://www.opensource.org/licenses/osl-3.0.php"""
34 __url__="http://www.iservo.edu.au/esys"
35 __version__="$Revision$"
36 __date__="$Date$"
37
38 from esys.escript import *
39 from esys.escript.pdetools import SaddlePointProblem
40 from esys.escript.linearPDEs import LinearPDE
41 from esys.finley import Brick
42 from esys.pyvisi import Scene, DataCollector, Contour, Camera, Velocity, Text2D, LocalPosition
43 from esys.pyvisi.constant import *
44 import os
45
46 JPG_RENDERER = Renderer.OFFLINE_JPG
47 ne = 20 # number of elements in spatial direction
48
49 def lockToGrid(x,l,n):
50 h=l/n
51 return int(x/h+0.5)*h
52
53
54 width = 100000.
55 height = 30000.
56
57 ne_w=int((ne/height)*width+0.5)
58 rho=0.
59 lmbd=1.7e11
60 mu=1.7e11
61 g=9.81
62
63 # fault location:
64
65 fstart = [lockToGrid(50000.0,width,ne_w), lockToGrid(40000.0,width,ne_w), lockToGrid(8000.,height,ne)]
66 fend = [lockToGrid(50000.0,width,ne_w), lockToGrid(60000.0,width,ne_w), lockToGrid(20000.,height,ne)]
67 s=[0.,1.,0.]
68
69 print "=== generate mesh over %s x %s x %s ==="%(width,width,height)
70 dom=Brick(l0=width, l1=width, l2=height, n0=ne_w, n1=ne_w, n2=ne)
71 print " total number of elements = ",ne*ne_w*ne_w
72
73 print "=== prepare PDE coefficients ==="
74 # fixed displacements:
75 d=dom.getDim()
76 x=dom.getX()
77 mask=whereZero(x[d-1]-inf(x[d-1]))*numarray.ones((d,))
78
79 # map the fault into the unit square:
80 x_hat=x-fstart
81 p=[0.,0.,0.]
82 q=1.
83 for i in range(d):
84 p[i]=fend[i]-fstart[i]
85 if abs(p[i])>0:
86 x_hat[i]/=p[i]
87 q=whereNonNegative(x_hat[i])*whereNonPositive(x_hat[i]-1.)*x_hat[i]*(1.-x_hat[i])*4.*q
88 else:
89 flip=i
90 q=whereZero(x_hat[i],1.e-3*Lsup(x_hat[i]))*q
91 g_s=grad(q*s)
92 sigma0=(mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d))*(whereNegative(interpolate(x_hat[flip],g_s.getFunctionSpace()))-0.5)
93 pde=LinearPDE(dom)
94 pde.setSymmetryOn()
95 A =Tensor4(0.,Function(dom))
96 for i in range(d):
97 for j in range(d):
98 A[i,j,j,i] += mu
99 A[i,j,i,j] += mu
100 A[i,i,j,j] += lmbd
101 pde.setValue(A=A,q=mask,Y=-kronecker(Function(dom))[d-1]*g*rho,X=sigma0)
102
103 print "=== solve pde === "
104 u=pde.getSolution(verbose=True)
105
106 print "=== calculate stress ==="
107 g_s=grad(u)
108 sigma=mu*symmetric(g_s)+lmbd*trace(g_s)*kronecker(d)
109
110 print "=== start rendering ==="
111 # Create a Scene.
112 s = Scene(renderer = JPG_RENDERER, num_viewport = 1, x_size = 800, y_size = 800)
113 dc1 = DataCollector(source = Source.ESCRIPT)
114 dc1.setData(disp=u,sigma=sigma,cfs=sigma[0,1]-0.4*sigma[0,0])
115
116 # Create a Contour.
117 ctr1 = Contour(scene = s, data_collector = dc1, viewport = Viewport.SOUTH_WEST,
118 lut = Lut.COLOR, cell_to_point = True, outline = True)
119 ctr1.generateContours(contours = 10)
120 ctr1.setOpacity(0.5)
121
122 # Create velocity
123 vopc1 = Velocity(scene = s, data_collector = dc1,
124 color_mode = ColorMode.VECTOR,
125 arrow = Arrow.THREE_D, lut = Lut.COLOR, cell_to_point = False,
126 outline = False)
127 vopc1.setScaleFactor(scale_factor = 10000.)
128 vopc1.setRatio(2)
129 vopc1.randomOn()
130
131 # Create a Camera.
132 cam1 = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
133 cam1.elevation(angle = -50)
134
135 # add some text
136 t2 = Text2D(scene = s, text = "CFS and displacement around vertical fault")
137 t2.setPosition(LocalPosition(200,30))
138 t2.setColor(color = Color.BLACK)
139 t2.setFontSize(size = 20)
140 t2.setFontToArial()
141 t2.shadowOn()
142
143 # Render the object.
144 s.render("stress_contour.jpg")
145 print "=== finished rendering ==="

  ViewVC Help
Powered by ViewVC 1.1.26