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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2588 - (hide annotations)
Wed Aug 5 22:42:36 2009 UTC (11 years, 5 months ago) by ahallam
Original Path: trunk/doc/cookbook/heatrefraction002.py
File MIME type: text/x-python
File size: 5248 byte(s)
Heat diffraction initial committal 2 weeks work.
1 ahallam 2588
2     ########################################################
3     #
4     # Copyright (c) 2003-2009 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7     #
8     # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11     #
12     ########################################################
13    
14     __copyright__="""Copyright (c) 2003-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     http://www.opensource.org/licenses/osl-3.0.php"""
20     __url__="https://launchpad.net/escript-finley"
21    
22     """
23     Author: Antony Hallam antony.hallam@uqconnect.edu.au
24     """
25    
26     # To solve the problem it is necessary to import the modules we require.
27     from esys.escript import * # This imports everything from the escript library
28     from esys.escript.linearPDEs import LinearPDE, Poisson # This defines LinearPDE as LinearPDE
29     from esys.finley import Rectangle, ReadMesh, Domain # This imports the rectangle domain function from finley
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.
32     import numpy as np
33     import pylab as pl
34    
35     from esys.escript.pdetools import *
36    
37     ##ESTABLISHING VARIABLES
38     #DOMAIN
39     mymesh = ReadMesh("heatrefraction_mesh003.fly")
40     #~ print Function(mymesh).getListOfTags()
41    
42     qin=70*Milli*W/(m*m) #our heat source temperature is now zero
43     Ti=290.15*K # Kelvin #the starting temperature of our iron bar
44    
45     # set up kappa (thermal conductivity across domain) using tags
46     kappa=Scalar(0,Function(mymesh))
47     kappa.setTaggedValue("tblockloop",2.0)
48     kappa.setTaggedValue("bblockloopl",3.0)
49     kappa.setTaggedValue("bblockloopr",4.0)
50    
51     #... generate functionspace...
52     #... open PDE ...
53     mypde=LinearPDE(mymesh)
54     #define first coefficient
55     mypde.setValue(A=kappa*kronecker(mymesh))
56    
57     # ... set initial temperature ....
58     x=mymesh.getX()
59    
60     qH=qin*whereZero(x[1]+6000)
61     mypde.setValue(q=Ti*whereZero(x[1]),r=Ti)
62     mypde.setValue(Y=qH,y=17*Celsius)
63    
64     # get steady state solution and export to vtk.
65     T=mypde.getSolution()
66     saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))
67    
68     # rearrage mymesh to suit solution function space
69     oldspacecoords=mymesh.getX()
70     coords=Data(oldspacecoords, T.getFunctionSpace())
71    
72     # calculate gradient of solution for quiver plot
73     qu=-kappa*grad(T)
74    
75     # create quiver locations
76     quivshape = [20,20] #quivers x and quivers y
77     numquiv = quivshape[0]*quivshape[1] # total number of quivers
78     maxx = 5000. # maximum x displacement
79     dx = maxx/quivshape[0]+1. # quiver x spacing
80     maxy = -6000. # maximum y displacement
81     dy = maxy/quivshape[1]+1. # quiver y spacing
82     qulocs=np.zeros([numquiv,2],float) # memory for quiver locations
83     # fill qulocs
84     for i in range(0,quivshape[0]-1):
85     for j in range(0,quivshape[1]-1):
86     qulocs[i*quivshape[0]+j,:] = [dx+dx*i,dy+dy*j]
87     # retreive values for quivers direction and shape from qu
88     quL = Locator(qu.getFunctionSpace(),qulocs.tolist())
89     qu = quL(qu) #array of dx,dy for quivers
90     qu = np.array(qu) #turn into a numpy array
91     qulocs = quL.getX() #returns actual locations from data
92     qulocs = np.array(qulocs) #turn into a numpy array
93    
94     kappaT = kappa.toListOfTuples(scalarastuple=False)
95     coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
96     coordsK = np.array(coordsK.toListOfTuples())
97     coordKX = coordsK[:,0]
98     coordKY = coordsK[:,1]
99    
100     coords = np.array(coords.toListOfTuples())
101     tempT = T.toListOfTuples(scalarastuple=False)
102    
103     coordX = coords[:,0]
104     coordY = coords[:,1]
105    
106     xi = np.linspace(0.0,5000.0,100)
107     yi = np.linspace(-6000.0,0.0,100)
108     # grid the data.
109     zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
110     ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
111     # contour the gridded data, plotting dots at the randomly spaced data points.
112    
113     pl.matplotlib.pyplot.autumn()
114     CKL = pl.contour(xi,yi,ziK,1,linewidths=1.0)
115     CK = pl.contourf(xi,yi,ziK,1)
116     CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
117     pl.clabel(CS, inline=1, fontsize=8)
118     CB = pl.colorbar(CS, shrink=0.8, extend='both')
119     pl.savefig("temp.png")
120    
121     #~ dx = np.zeros(82,float)
122     #~ dy = np.zeros(82,float)
123     #~ vlocs = np.zeros([82,2],float)
124     #~ delx = 5000.0/100
125     #~ dely = 6000.0/100
126     #~ ind = 0
127     #~
128     #~ for i in range(1,10):
129     #~ for j in range(1,10):
130     #~ iloc = i*10
131     #~ jloc = j*10
132     #~ ind += 1
133     #~
134     #~ vlocs[ind,0] = iloc*delx; vlocs[ind,1]=-jloc*dely
135     #~
136     #~ print zi[iloc,jloc], zi[iloc+5,jloc]
137     #~ print zi[iloc,jloc], zi[iloc,jloc+5]
138     #~ dx[ind] = (zi[iloc,jloc] - zi[iloc+5,jloc])/delx
139     #~ dy[ind] = (zi[iloc,jloc] - zi[iloc,jloc+5])/dely
140    
141     #~ zpoints = [zi[iloc,jloc],zi[iloc+2,jloc],zi[iloc+4,jloc]]
142     #~ xpoints = [vlocs[ind,0],vlocs[ind+2,0],vlocs[ind+4,0]]
143     #~ print pl.matplotlib.mlab.slopes(xpoints,zpoints)
144     #~ dx[ind] = pl.matplotlib.mlab.slopes(xpoints,zpoints)
145     #~
146     #~ ypoints = [vlocs[ind,1],vlocs[ind+2,1],vlocs[ind+4,]]
147     #~ zpoints = [zi[iloc,jloc],zi[iloc,jloc+2],zi[iloc,jloc+4]]
148     #~ dy[ind] = pl.matplotlib.mlab.slopes(ypoints,zpoints)
149     #~
150     pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="green")
151     pl.savefig("temp2.png")

  ViewVC Help
Powered by ViewVC 1.1.26