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

revision 2657 by jfenwick, Mon Sep 7 00:06:15 2009 UTC revision 2658 by ahallam, Thu Sep 10 02:58:44 2009 UTC
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
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()
108    ##########################################################VISUALISATION
109    # calculate gradient of solution for quiver plot
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

# 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\