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

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