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

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

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

revision 2657 by jfenwick, Mon Sep 7 00:06:15 2009 UTC revision 2658 by ahallam, Thu Sep 10 02:58:44 2009 UTC
# Line 23  __url__="https://launchpad.net/escript-f Line 23  __url__="https://launchpad.net/escript-f
23  Author: Antony Hallam antony.hallam@uqconnect.edu.au  Author: Antony Hallam antony.hallam@uqconnect.edu.au
24  """  """
25    
26  # To solve the problem it is necessary to import the modules we require.  ############################################################FILE HEADER
27  from esys.escript import * # This imports everything from the escript library  # heatrefraction002.py
28  from esys.escript.linearPDEs import LinearPDE, Poisson # This defines LinearPDE as LinearPDE  # Model steady state temperature distribution in three block model,
29  from esys.finley import Rectangle, ReadMesh, Domain # This imports the rectangle domain function from finley  # mesh from heatrefraction_mesher002.py
30  import os #This package is necessary to handle saving our data.  
31  from esys.escript.unitsSI import * # A useful unit handling package which will make sure all our units match up in the equations.  #######################################################EXTERNAL MODULES
32    # To solve the problem it is necessary to import the modules we
33    # require.
34    # This imports everything from the escript library
35    from esys.escript import *
36    # This defines LinearPDE as LinearPDE
37    from esys.escript.linearPDEs import LinearPDE, Poisson
38    # This imports the rectangle domain function from finley
39    from esys.finley import Rectangle, ReadMesh, Domain
40    # This package is necessary to handle saving our data.
41    import os
42    # A useful unit handling package which will make sure all our units
43    # match up in the equations.
44    from esys.escript.unitsSI import *
45    # numpy for array handling
46  import numpy as np  import numpy as np
47    import matplotlib
48    #For interactive use, you can comment out the next two lines
49    import matplotlib
50    matplotlib.use('agg') #It's just here for automated testing
51    # pylab for matplotlib and plotting
52  import pylab as pl  import pylab as pl
53    # cblib functions
54    from cblib import toQuivLocs, toXYTuple, needdirs
55    
56  from esys.escript.pdetools import *  ########################################################MPI WORLD CHECK
57  from cblib import needdirs  if getMPISizeWorld() > 1:
58        import sys
59        print "This example will not run in an MPI world."
60        sys.exit(0)
61    
62  ##ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
63  qin=300.*Milli*W/(m*m) #our heat source temperature is now zero  qin=300.*Milli*W/(m*m) #our heat source temperature is now zero
64  Ti=290.15*K # Kelvin #the starting temperature of our iron bar  Ti=290.15*K # Kelvin #the starting temperature of our iron bar
65    
66  # the folder to gett our outputs from, leave blank "" for script path -  # the folder to gett our outputs from, leave blank "" for script path -
67  # note these depen. are generated from heatrefraction_mesher001.py  # note these depen. are generated from heatrefraction_mesher001.py
68  saved_path = "data/heatrefrac002"  saved_path = save_path= os.path.join("data","heatrefrac002")
69  needdirs([saved_path])  needdirs([saved_path])
70    
71  ###### 2 BLOCK MODEL #########  ################################################ESTABLISHING PARAMETERS
72  ## DOMAIN  ## DOMAIN
 ## Anticline  
73  mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh003.fly"))  mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh003.fly"))
74  tpg = np.loadtxt(os.path.join(saved_path,"toppg"))  tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
75  tpgx = tpg[:,0]  tpgx = tpg[:,0]
# Line 79  qH.setTaggedValue("linebottom",qin) Line 102  qH.setTaggedValue("linebottom",qin)
102  mypde.setValue(q=whereZero(x[1]),r=Ti)  mypde.setValue(q=whereZero(x[1]),r=Ti)
103  mypde.setValue(y=qH)  mypde.setValue(y=qH)
104    
105  # get steady state solution and export to vtk.  ###########################################################GET SOLUTION
106  T=mypde.getSolution()  T=mypde.getSolution()
107  #saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))  
108    ##########################################################VISUALISATION
109    # calculate gradient of solution for quiver plot
110    qu=-kappa*grad(T)
111    
112  # rearrage mymesh to suit solution function space        # rearrage mymesh to suit solution function space      
113  oldspacecoords=mymesh.getX()  oldspacecoords=mymesh.getX()
114  coords=Data(oldspacecoords, T.getFunctionSpace())  coords=Data(oldspacecoords, T.getFunctionSpace())
115    
 # calculate gradient of solution for quiver plot  
 qu=-kappa*grad(T)  
   
 # create quiver locations  
116  quivshape = [20,20] #quivers x and quivers y  quivshape = [20,20] #quivers x and quivers y
117  numquiv = quivshape[0]*quivshape[1] # total number of quivers  # function to calculate quiver locations
118  maxx = 5000. # maximum x displacement  qu,qulocs = toQuivLocs(quivshape,width,depth,qu)
 dx = maxx/quivshape[0]+1. # quiver x spacing  
 maxy = -6000. # maximum y displacement  
 dy = maxy/quivshape[1]+1. # quiver y spacing  
 qulocs=np.zeros([numquiv,2],float) # memory for quiver locations  
 # fill qulocs  
 for i in range(0,quivshape[0]-1):  
     for j in range(0,quivshape[1]-1):  
         qulocs[i*quivshape[0]+j,:] = [dx+dx*i,dy+dy*j]  
 # retreive values for quivers direction and shape from qu  
 quL = Locator(qu.getFunctionSpace(),qulocs.tolist())  
 qu = quL(qu) #array of dx,dy for quivers  
 qu = np.array(qu) #turn into a numpy array  
 qulocs = quL.getX() #returns actual locations from data  
 qulocs = np.array(qulocs) #turn into a numpy array  
119    
120  kappaT = kappa.toListOfTuples(scalarastuple=False)  kappaT = kappa.toListOfTuples(scalarastuple=False)
121  coordsK = Data(oldspacecoords, kappa.getFunctionSpace())  coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
122  coordsK = np.array(coordsK.toListOfTuples())  coordKX, coordKY = toXYTuple(coordsK)
 coordKX = coordsK[:,0]  
 coordKY = coordsK[:,1]  
123                
 coords = np.array(coords.toListOfTuples())  
124  tempT = T.toListOfTuples(scalarastuple=False)  tempT = T.toListOfTuples(scalarastuple=False)
125    coordX, coordY = toXYTuple(coords)
126    
127  coordX = coords[:,0]  xi = np.linspace(0.0,width,100)
128  coordY = coords[:,1]  yi = np.linspace(depth,0.0,100)
   
 xi = np.linspace(0.0,5000.0,100)  
 yi = np.linspace(-6000.0,0.0,100)  
129  # grid the data.  # grid the data.
130  zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)  zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
131  ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)  ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
132  # contour the gridded data, plotting dots at the randomly spaced data points.  # contour the gridded data, plotting dots at the randomly
133    # spaced data points.
134    
135  pl.matplotlib.pyplot.autumn()  pl.matplotlib.pyplot.autumn()
136  CKL = pl.fill(tpgx,tpgy,'brown',bpglx,bpgly,'red',bpgrx,bpgry,'orange',zorder=-1000)  CKL = pl.fill(tpgx,tpgy,'brown',bpglx,bpgly,'red',\
137  #~ CK = pl.contourf(xi,yi,ziK,2)                bpgrx,bpgry,'orange',zorder=-1000)
138  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
139  pl.clabel(CS, inline=1, fontsize=8)  pl.clabel(CS, inline=1, fontsize=8)
140  pl.title("Heat Refraction across an anisotropic structure.")  pl.title("Heat Refraction across an anisotropic structure.")
141  pl.xlabel("Horizontal Displacement (m)")  pl.xlabel("Horizontal Displacement (m)")
142  pl.ylabel("Depth (m)")  pl.ylabel("Depth (m)")
143  #~ CB = pl.colorbar(CS, shrink=0.8, extend='both')  if getMPIRankWorld() == 0: #check for MPI processing
144  pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))      pl.savefig(os.path.join(saved_path,"heatrefraction002_cont.png"))
145    
146  QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white")  QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],\
147  pl.title("Heat Refraction across an anisotropic structure \n with gradient quivers.")                  angles='xy',color="white")
148  pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))  pl.title("Heat Refraction across an anisotropic structure\
149                \n with gradient quivers.")
150    if getMPIRankWorld() == 0: #check for MPI processing
151        pl.savefig(os.path.join(saved_path,"heatrefraction002_contqu.png"))

Legend:
Removed from v.2657  
changed lines
  Added in v.2658

  ViewVC Help
Powered by ViewVC 1.1.26