# Contents of /trunk/doc/examples/cookbook/example10a.py

Revision 3148 - (show annotations)
Fri Sep 3 02:09:47 2010 UTC (8 years, 6 months ago) by jfenwick
File MIME type: text/x-python
File size: 3823 byte(s)
```Another attempt to patch the X issue

```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2009 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2009 by University of Queensland 15 Earth Systems Science Computational Center (ESSCC) 16 http://www.uq.edu.au/esscc 17 Primary Business: Queensland, Australia""" 18 __license__="""Licensed under the Open Software License version 3.0 19 20 __url__= 21 22 """ 23 Author: Antony Hallam antony.hallam@uqconnect.edu.au 24 """ 25 26 ############################################################FILE HEADER 27 # example10a.py 28 # Model of gravitational Potential. 29 30 #######################################################EXTERNAL MODULES 31 # To solve the problem it is necessary to import the modules we require. 32 from esys.escript import * # This imports everything from the escript library 33 from esys.escript.unitsSI import * 34 from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE 35 from esys.finley import Rectangle # This imports the rectangle domain function from finley 36 import os, sys #This package is necessary to handle saving our data. 37 from math import pi, sqrt, sin, cos 38 39 from esys.escript.pdetools import Projector 40 from cblib import toRegGrid 41 42 import matplotlib 43 matplotlib.use('agg') #It's just here for automated testing 44 45 import pylab as pl #Plotting package 46 import numpy as np 47 48 ########################################################MPI WORLD CHECK 49 if getMPISizeWorld() > 1: 50 import sys 51 print "This example will not run in an MPI world." 52 sys.exit(0) 53 54 #################################################ESTABLISHING VARIABLES 55 #Domain related. 56 mx = 5000*m #meters - model length 57 my = -5000*m #meters - model width 58 ndx = 100 # mesh steps in x direction 59 ndy = 100 # mesh steps in y direction - one dimension means one element 60 #PDE related 61 rho=200.0 62 rholoc=[2500,-2500] 63 G=6.67300*10E-11 64 65 ################################################ESTABLISHING PARAMETERS 66 #the folder to put our outputs in, leave blank "" for script path 67 save_path= os.path.join("data","example10") 68 #ensure the dir exists 69 mkDir(save_path) 70 71 ####################################################DOMAIN CONSTRUCTION 72 domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 73 x=Solution(domain).getX() 74 mask=wherePositive(10-length(x-rholoc)) 75 rho=rho*mask 76 kro=kronecker(domain) 77 78 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx) 79 ###############################################ESCRIPT PDE CONSTRUCTION 80 81 mypde=LinearPDE(domain) 82 mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0) 83 mypde.setSymmetryOn() 84 sol=mypde.getSolution() 85 86 g_field=grad(sol) #The graviational accelleration g. 87 g_fieldz=g_field*[0,1] # 88 gz=length(g_fieldz) 89 # Save the output to file. 90 saveVTK(os.path.join(save_path,"ex10a.vtu"),\ 91 grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz) 92 93 ##################################################REGRIDDING & PLOTTING 94 95 96 xi, yi, zi = toRegGrid(sol, nx=50, ny=50) 97 pl.matplotlib.pyplot.autumn() 98 pl.contourf(xi,yi,zi,10) 99 pl.xlabel("Horizontal Displacement (m)") 100 pl.ylabel("Depth (m)") 101 pl.savefig(os.path.join(save_path,"Ucontour.png")) 102 print "Solution has been plotted ..." 103 104 cut=int(len(xi)/2) 105 106 pl.clf() 107 108 r=np.linspace(0,mx/2,100) 109 m=2*pl.pi*10*10*200*-G/(r*r) 110 111 pl.plot(xi,zi[:,cut]) 112 #pl.plot(r+2500,m) 113 pl.title("Potential Profile") 114 pl.xlabel("Horizontal Displacement (m)") 115 pl.ylabel("Potential") 116 pl.savefig(os.path.join(save_path,"Upot00.png")) 117 118 out=np.array([xi,zi[:,cut]]) 119 pl.savetxt('profile1.asc',out.transpose()) 120 pl.clf()