/[escript]/trunk/doc/examples/cookbook/example10a.py
ViewVC logotype

Diff of /trunk/doc/examples/cookbook/example10a.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3134 by ahallam, Sun Aug 22 04:35:15 2010 UTC revision 3135 by ahallam, Thu Sep 2 03:28:39 2010 UTC
# Line 32  Author: Antony Hallam antony.hallam@uqco Line 32  Author: Antony Hallam antony.hallam@uqco
32  from esys.escript import * # This imports everything from the escript library  from esys.escript import * # This imports everything from the escript library
33  from esys.escript.unitsSI import *  from esys.escript.unitsSI import *
34  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
 from esys.escript.linearPDEs import Poisson # This defines LinearPDE as LinearPDE  
35  from esys.finley import Rectangle # This imports the rectangle domain function from finley  from esys.finley import Rectangle # This imports the rectangle domain function from finley
 #For interactive use, you can comment out the next two lines  
 import matplotlib  
 matplotlib.use('agg') #It's just here for automated testing  
 import pylab as pl #Plotting package.  
 import numpy as np #Array package.  
36  import os, sys #This package is necessary to handle saving our data.  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    import pylab as pl #Plotting package
42    import numpy as np
43    
44  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
45  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
# Line 49  if getMPISizeWorld() > 1: Line 49  if getMPISizeWorld() > 1:
49    
50  #################################################ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
51  #Domain related.  #Domain related.
52  mx = 500*m #meters - model length  mx = 5000*m #meters - model length
53  my = 500*m #meters - model width  my = -5000*m #meters - model width
54  ndx = 500 # mesh steps in x direction  ndx = 100 # mesh steps in x direction
55  ndy = 500 # mesh steps in y direction - one dimension means one element  ndy = 100 # mesh steps in y direction - one dimension means one element
56  #PDE related  #PDE related
57  rho=100.0  rho=200.0
58  rholoc=[250,250]  rholoc=[2500,-2500]
59  G=6.67300*10E-11  G=6.67300*10E-11
60    
61  ################################################ESTABLISHING PARAMETERS  ################################################ESTABLISHING PARAMETERS
62  #the folder to put our outputs in, leave blank "" for script path  #the folder to put our outputs in, leave blank "" for script path
63  save_path= os.path.join("data","example10a")  save_path= os.path.join("data","example10")
64  #ensure the dir exists  #ensure the dir exists
65  mkDir(save_path)  mkDir(save_path)
66    
67  ####################################################DOMAIN CONSTRUCTION  ####################################################DOMAIN CONSTRUCTION
68  domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
69  x=Solution(domain).getX()  x=Solution(domain).getX()
70  #mask=whereZero(x[0]-rholoc[0])*whereZero(x[1]-rholoc[1])  mask=wherePositive(10-length(x-rholoc))
 mask=wherePositive(100-length(x-rholoc))  
 #rho=Scalar(rho,ContinuousFunction(domain))  
71  rho=rho*mask  rho=rho*mask
 y=Scalar(0.0,FunctionOnBoundary(domain))  
72  kro=kronecker(domain)  kro=kronecker(domain)
73    
74  q=whereZero(x[0])+\  q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
   whereZero(x[1])+\  
   whereZero(x[0]-sup(x[0]))+\  
   whereZero(x[1]-sup(x[1]))  
75  ###############################################ESCRIPT PDE CONSTRUCTION  ###############################################ESCRIPT PDE CONSTRUCTION
76    
77  mypde=LinearPDE(domain)  mypde=LinearPDE(domain)
78  mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0.)  mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0)
79  #mypde.setSymmetryOn()  mypde.setSymmetryOn()
 #mypde=Poisson(domain)  
 #mypde.setValue(f=rho*4.*3.1415*G,q=q)  
80  sol=mypde.getSolution()  sol=mypde.getSolution()
81  saveVTK("ex10a.vtu",field_strength=sol)#,field=-grad(sol))  
82    g_field=grad(sol) #The graviational accelleration g.
83    g_fieldz=g_field*[0,1] #
84    gz=length(g_fieldz)
85    # Save the output to file.
86    saveVTK(os.path.join(save_path,"ex10a.vtu"),\
87            grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
88    
89    ##################################################REGRIDDING & PLOTTING
90    
91    
92    xi, yi, zi = toRegGrid(sol, nx=50, ny=50)
93    pl.matplotlib.pyplot.autumn()
94    pl.contourf(xi,yi,zi,10)
95    pl.xlabel("Horizontal Displacement (m)")
96    pl.ylabel("Depth (m)")
97    pl.savefig(os.path.join(save_path,"Ucontour.png"))
98    print "Solution has been plotted  ..."
99    
100    cut=int(len(xi)/2)
101    
102    pl.clf()
103    
104    r=np.linspace(0,mx/2,100)
105    m=2*pl.pi*10*10*200*-G/(r*r)
106    
107    pl.plot(xi,zi[:,cut])
108    #pl.plot(r+2500,m)
109    pl.title("Potential Profile")
110    pl.xlabel("Horizontal Displacement (m)")
111    pl.ylabel("Potential")
112    pl.savefig(os.path.join(save_path,"Upot00.png"))
113    
114    out=np.array([xi,zi[:,cut]])
115    pl.savetxt('profile1.asc',out.transpose())
116    pl.clf()

Legend:
Removed from v.3134  
changed lines
  Added in v.3135

  ViewVC Help
Powered by ViewVC 1.1.26