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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2588 - (show 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
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