1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032010 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/osl3.0.php 
11 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 

14 
\section{Steadystate Heat Refraction} 
15 
\label{STEADYSTATE HEAT REFRACTION} 
16 

17 
In this chapter we demonstrate how to handle more complex geometries. 
18 

19 
Steadystate heat refraction will give us an opportunity to investigate some of the richer features that the \esc package has to offer. One of these is \pycad . The advantage of using \pycad is that it offers an easy method for developing and manipulating complex domains. In conjunction with \gmsh we can generate finite element meshes that conform to our domain's shape providing accurate modelling of interfaces and boundaries. Another useful function of \pycad is that we can tag specific areas of our domain with labels as we construct them. These labels can then be used in \esc to define properties like material constants and source locations. 
20 

21 
We proceed in this chapter by first looking at a very simple geometry. In fact, we solve 
22 
the steadystate heat equation over a rectangular domain. This is not a very interesting problem but then we extend our script by introducing an internal, curved interface. 
23 

24 
\begin{figure}[ht] 
25 
\centerline{\includegraphics[width=4.in]{figures/pycadrec}} 
26 
\caption{Example 4: Rectangular Domain for \pycad.} 
27 
\label{fig:pycad rec} 
28 
\end{figure} 
29 

30 
\section{Example 4: Creating the domain with \pycad} 
31 
\sslist{example04a.py} 
32 

33 
We modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways. Firstly, we look at the steady state 
34 
case with slightly modified boundary conditions and then we use a more flexible tool 
35 
to generate the geometry. Lets look at the geometry first. 
36 

37 
We want to define a rectangular domain of width $5 km$ and depth $6 km$ below the surface of the Earth. The domain is subject to a few conditions. The temperature is known at the surface and the basement has a known heat flux. Each side of the domain is insulated and the aim is to calculate the final temperature distribution. 
38 

39 
In \pycad there are a few primary constructors that build upon each other to define domains and boundaries; 
40 
the ones we use are: 
41 
\begin{python} 
42 
from esys.pycad import * 
43 
Point() #Create a point in space. 
44 
Line() #Creates a line from a number of points. 
45 
CurveLoop() #Creates a closed loop from a number of lines. 
46 
PlaneSurface() #Creates a surface based on a CurveLoop 
47 
\end{python} 
48 
So to construct our domain as shown in \reffig{fig:pycad rec}, we first need to create 
49 
the corner points. From the corner points we build the four edges of the rectangle. The four edges 
50 
then form a closed loop which defines our domain as a surface. 
51 
We start by inputting the variables we need to construct the model. 
52 
\begin{python} 
53 
width=5000.0*m #width of model 
54 
depth=6000.0*m #depth of model 
55 
\end{python} 
56 
The variables are then used to construct the four corners of our domain, which from the origin has the dimensions of 5000 meters width and 6000 meters depth. This is done with the \verb Point() function which accepts x, y and z coordinates. Our domain is in two dimensions so z should always be zero. 
57 
\begin{python} 
58 
# Overall Domain 
59 
p0=Point(0.0, 0.0, 0.0) 
60 
p1=Point(0.0, depth, 0.0) 
61 
p2=Point(width, depth, 0.0) 
62 
p3=Point(width, 0.0, 0.0) 
63 
\end{python} 
64 
Now lines are defined using our points. This forms a rectangle around our domain; 
65 
\begin{python} 
66 
l01=Line(p0, p1) 
67 
l12=Line(p1, p2) 
68 
l23=Line(p2, p3) 
69 
l30=Line(p3, p0) 
70 
\end{python} 
71 
Notice that lines have a direction. These lines form the basis for our domain boundary, which is a closed loop. 
72 
\begin{python} 
73 
c=CurveLoop(l01, l12, l23, l30) 
74 
\end{python} 
75 
Be careful to define the curved loop in an \textbf{anticlockwise} manner otherwise the meshing algorithm may fail. 
76 
Finally we can define the domain as 
77 
\begin{python} 
78 
rec = PlaneSurface(c) 
79 
\end{python} 
80 
At this point the introduction of the curved loop seems to be unnecessary but this concept plays an important role 
81 
if holes are introduced. 
82 

83 
Now we are ready to handover the domain \verbrec to a mesher which subdivides the domain into triangles (or tetrahedron in 3D). In our case we use \gmsh. We create 
84 
an instance of the \verbDesign class which will handle the interface to \gmsh: 
85 
\begin{python} 
86 
from esys.pycad.gmsh import Design 
87 
d=Design(dim=2, element_size=200*m) 
88 
\end{python} 
89 
The argument \verbdim defines the spatial dimension of the domain\footnote{If \texttt{dim}=3 the rectangle would be interpreted as a surface in the three dimensional space.}. The second argument \verbelement_size defines element size which is the maximum length of a triangle edge in the mesh. The element size needs to be chosen with care in order to avoid very large meshes if you don't want to. In our case with an element size of $200$m 
90 
and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$ triangles in each spatial direction. So we end up with about $30 \times 30 = 900$ triangles which is a size that can be handled easily. 
91 
We can easily add our domain \verbrec to the \verbDesign; 
92 
\begin{python} 
93 
d.addItem(rec) 
94 
\end{python} 
95 
We have the plan to set a heat flux on the bottom of the domain. One can use the masking technique to do this 
96 
but \pycad offers a more convenient technique called tagging. With this technique items in the domain are 
97 
named using the \verbPropertySet class. We can then later use this name to set values secifically for 
98 
those sample points located on the named items. Here we name the bottom face of the 
99 
domain where we will set the heat influx: 
100 
\begin{python} 
101 
ps=PropertySet("linebottom",l12)) 
102 
d.addItem(ps) 
103 
\end{python} 
104 
Now we are ready to hand over the \verbDesign to \FINLEY: 
105 
\begin{python} 
106 
from esys.finley import MakeDomain 
107 
domain=MakeDomain(d) 
108 
\end{python} 
109 
The \verbdomain object can now be used in the same way like the return object of the \verbRectangle 
110 
object we have used previously to generate a mesh. It is common practice to separate the 
111 
mesh generation from the PDE solution. The main reason for this is that mesh generation can be computationally very expensive in particular in 3D. So it is more efficient to generate the mesh once and write it to a file. The mesh 
112 
can then be read in every time a new simulation is run. \FINLEY supports this in the following 
113 
way~\footnote{An alternative are the \texttt{dump} and \texttt{load} functions. They using a binary format and tend to be much smaller.} 
114 
\begin{python} 
115 
# write domain to an text file 
116 
domain.write("example04.fly") 
117 
\end{python} 
118 
and then for reading in another script; 
119 
\begin{python} 
120 
# read domain from text file 
121 
from esys.finley import ReadMesh 
122 
domain =ReadMesh("example04.fly") 
123 
\end{python} 
124 

125 
\begin{figure}[ht] 
126 
\centerline{\includegraphics[width=4.in]{figures/simplemesh}} 
127 
\caption{Example 4a: Mesh over rectangular domain, see \reffig{fig:pycad rec}.} 
128 
\label{fig:pycad rec mesh} 
129 
\end{figure} 
130 

131 
Before we discuss how we solve the PDE for the 
132 
problem it is useful to present two additional options of the \verbDesign. 
133 
They allow the user accessing the script which is used by \gmsh to generate the mesh as well as 
134 
the mesh as it has been generated by \gmsh. This is done by setting specific names for these files: 
135 
\begin{python} 
136 
d.setScriptFileName("example04.geo") 
137 
d.setMeshFileName("example04.msh") 
138 
\end{python} 
139 
Usually the extension \texttt{geo} is used for the script file of the \gmsh geometry and 
140 
the extension \texttt{msh} for the mesh file. Normally these files are deleted after usage. 
141 
Accessing these file can be helpful to debug the generation of more complex geometries. The geometry and the mesh can be visualised from the command line using 
142 
\begin{verbatim} 
143 
gmsh example04.geo # show geometry 
144 
gmsh example04.msh # show mesh 
145 
\end{verbatim} 
146 
The mesh is shown in \reffig{fig:pycad rec mesh}. 
147 
\begin{figure}[ht] 
148 
\centerline{\includegraphics[width=4.in]{figures/simpleheat}} 
149 
\caption{Example 4b: Result of simple steady state heat problem.} 
150 
\label{fig:steady state heat} 
151 
\end{figure} 
152 

153 

154 
\section{The Steadystate Heat Equation} 
155 
\sslist{example04b.py, cblib} 
156 

157 
Steady state is reached in the temperature when it is not changing in time. So to calculate the 
158 
the steady state we set the time derivative term in \refEq{eqn:Tform nabla} to zero; 
159 
\begin{equation}\label{eqn:Tform nabla steady} 
160 
\nabla \cdot \kappa \nabla T = q\hackscore H 
161 
\end{equation} 
162 
This PDE is easier to solve, than the PDE in \refEq{eqn:hdgenf2} in each time step. We can just drop 
163 
the \verbD term; 
164 
\begin{python} 
165 
mypde=LinearPDE(domain) 
166 
mypde.setValue(A=kappa*kronecker(model), Y=qH) 
167 
\end{python} 
168 
The temperature at the top face of the domain is known as \verbTtop (=$20 C$). In \refSec{Sec:1DHDv0} we have 
169 
already discussed how this constraint is added to the PDE: 
170 
\begin{python} 
171 
x=Solution(domain).getX() 
172 
mypde.setValue(q=whereZero(x[1]sup(x[1])),r=Ttop) 
173 
\end{python} 
174 
Notice that we use the \verbsup function to calculate the maximum of $y$ coordinates of the relevant sample points. 
175 

176 
In all cases so far we have assumed that the domain is insulated which translates 
177 
into a zero normal flux $n \cdot \kappa \nabla T$, see \refEq{eq:hom flux}. In the modeling 
178 
setup of this chapter we want to set the normal heat flux at the bottom to \verbqin while we still 
179 
maintain insulation at the left and right face. Mathematically we can express this as the format 
180 
\begin{equation} 
181 
n \cdot \kappa \nabla T = q\hackscore{S} 
182 
\label{eq:inhom flux} 
183 
\end{equation} 
184 
where $q\hackscore{S}$ is a function of its location on the boundary. Its value becomes zero 
185 
for locations on the left or right face of the domain while it has the value \verbqin at the bottom face. 
186 
Notice that the value of $q\hackscore{S}$ at the top face is not relevant as we prescribe the temperature here. 
187 
We could define $q\hackscore{S}$ by using the masking techniques demonstrated earlier. The tagging mechanism provides an alternative and in many cases more convenient way of defining piecewise 
188 
constant functions such as $q\hackscore{S}$. You need to remember now that we 
189 
marked the bottom face with the name \verblinebottom when we defined the domain. 
190 
We can use this now to create $q\hackscore{S}$; 
191 
\begin{python} 
192 
qS=Scalar(0,FunctionOnBoundary(domain)) 
193 
qS.setTaggedValue("linebottom",qin) 
194 
\end{python} 
195 
In the first line \verbqS is defined as a scalar value over the sample points on the boundary of the domain. It is 
196 
initialized to zero for all sample points. In the second statement the values for those sample points which 
197 
on the line marked by \verblinebottom are set to \verbqin. 
198 

199 
The Neuman boundary condition assumed by \esc has in fact the form 
200 
\begin{equation}\label{NEUMAN 2b} 
201 
n\cdot A \cdot\nabla u = y 
202 
\end{equation} 
203 
In comparison to the version in \refEq{NEUMAN 2} we have used so far the right hand side is now 
204 
the new PDE coefficient $y$. As we have not specific $y$ in our previous examples \esc has assumed 
205 
the value zero for $y$. A comparison of \refEq{NEUMAN 2b} and \refEq{eq:inhom flux} reveals that one need to 
206 
choose $y=q\hackscore{S}$; 
207 
\begin{python} 
208 
qS=Scalar(0,FunctionOnBoundary(domain)) 
209 
qS.setTaggedValue("linebottom",qin) 
210 
mypde.setValue(y=qS) 
211 
\end{python} 
212 
To plot the results we are using the \modmpl library as shown \refSec{Sec:2DHD plot}. For convenience 
213 
the interpolation of the temperature to a rectangular grid for contour plotting is made available 
214 
via the \verbtoRegGrid function in the \verbcblib module. Your result should look similar to 
215 
\reffig{fig:steady state heat}. 
216 

217 
\begin{figure}[ht] 
218 
\centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}} 
219 
\caption{Example 5a: Heat refraction model with point and line labels.} 
220 
\label{fig:anticlinehrmodel} 
221 
\end{figure} 
222 

223 
\begin{figure}[ht] 
224 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction}} 
225 
\caption{Example 5a: Temperature Distribution in the Heat Refraction Model.} 
226 
\label{fig:anticlinetemp} 
227 
\end{figure} 
228 

229 
\section{Example 5: A Heat Refraction model} 
230 
\sslist{example05a.py and cblib.py} 
231 

232 
Our heat refraction model will be a large anticlinal structure that is experiencing a constant temperature at the surface and a steady heat flux at it's base. Our aim is to show that the temperature flux across the surface is not linear from bottom to top but is in fact warped by the structure of the model and is heavily dependent upon the material properties and shape. 
233 

234 

235 
We modify the example of the previous section by subdividing the block into two parts. The curve 
236 
separating the two blocks is given as a spline, see \reffig{fig:anticlinehrmodel}. The data points 
237 
used to define the curve may be measurement but for simplicity we assume here that there coordinates are 
238 
known in analytic form. 
239 

240 
There are two modes available in our example. When \verb modal=1 this indicates to the script that we wish to model an anticline. Otherwise when \verb modal=1 this will model a syncline. The modal operator simply changes the orientation of the boundary function so that it is either upwards or downwards curving. A \verb save_path has also been defined so that we can easily separate our data from other examples and our scripts. 
241 

242 
It is now possible to start defining our domain and boundaries. 
243 

244 
The curve defining our clinal structure is located approximately in the middle of the domain and has a sinusoidal shape. We define the curve by generating points at discrete intervals; $51$ in this case, and then create a smooth curve through the points using the \verb Spline() function. 
245 
\begin{python} 
246 
# Material Boundary 
247 
x=[ Point(i*dsp\ 
248 
,dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\ 
249 
for i in range(0,sspl)\ 
250 
] 
251 
mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple 
252 
\end{python} 
253 
The start and end points of the spline can be returned to help define the material boundaries. 
254 
\begin{python} 
255 
x1=mysp.getStartPoint() 
256 
x2=mysp.getEndPoint() 
257 
\end{python} 
258 
The top block or material above the clinal/spline boundary is defined in a \textbf{anticlockwise} manner by creating lines and then a closed loop. As we will be meshing the subdomain we also need to assign it a planar surface. 
259 
\begin{python} 
260 
# TOP BLOCK 
261 
tbl1=Line(p0,x1) 
262 
tbl2=mysp 
263 
tbl3=Line(x2,p3) 
264 
l30=Line(p3, p0) 
265 
tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30) 
266 
tblock = PlaneSurface(tblockloop) 
267 
\end{python} 
268 
We must repeat the above for every other subdomain. In this example there is only one other, the bottom block. The process is fairly similar to the top block but with a few differences. The spline points must be reversed by setting the spline as negative. 
269 
\begin{python} 
270 
bbl4=mysp 
271 
\end{python} 
272 
This reverse spline option unfortunately does not work for the getLoopCoords command, however, the \modmpl polygon tool will accept clockwise oriented points so we can define a new curve. 
273 
\begin{python} 
274 
#clockwise check 
275 
bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 
276 
\end{python} 
277 
The last few steps in creating the model are: take the previously defined domain and subdomain points and generate a mesh that can be imported into \esc. 
278 
To initialise the mesh it first needs some design parameters. In this case we have 2 dimensions \verb dim and a specified number of finite elements that need to be applied to the domain \verb element_size . It then becomes a simple task of adding the subdomains and flux boundaries to the design. Each element of our model can be given an identifier which makes it easier to define the subdomain properties in the solution script. This is done using the \verb PropertySet() function. The geometry and mesh are then saved so the \esc domain can be created. 
279 
\begin{python} 
280 
# Create a Design which can make the mesh 
281 
d=Design(dim=2, element_size=200) 
282 
# Add the subdomains and flux boundaries. 
283 
d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 
284 
PropertySet("linebottom",l12)) 
285 
# Create the geometry, mesh and \esc domain 
286 
d.setScriptFileName(os.path.join(save_path,"example05.geo")) 
287 
d.setMeshFileName(os.path.join(save_path,"example05.msh")) 
288 
domain=MakeDomain(d,optimizeLabeling=True) 
289 
\end{python} 
290 
The creation of our domain and its mesh is complete. 
291 

292 
With the mesh imported it is now possible to use our tagging property to set up our PDE coefficients. In this case $\kappa$ is set via the \verb setTaggedValue() function which takes two arguments, the name of the tagged points and the value to assign to them. 
293 
\begin{python} 
294 
# set up kappa (thermal conductivity across domain) using tags 
295 
kappa=Scalar(0,Function(domain)) 
296 
kappa.setTaggedValue("top",2.0*W/m/K) 
297 
kappa.setTaggedValue("bottom",4.0*W/m/K) 
298 
\end{python} 
299 
No further changes are required to the PDE solution step, see \reffig{fig:anticlinetemp} for the result. 
300 

301 
\begin{figure} 
302 
\centering 
303 
\subfigure[Temperature Depth Profile]{\includegraphics[width=3in]{figures/heatrefractiontdp.png}\label{fig:tdp}} 
304 
\subfigure[Temperature Gradient Depth Profile]{\includegraphics[width=3in]{figures/heatrefractiontgdp.png}\label{fig:tgdp}} 
305 
\subfigure[Thermal Conductivity Profile]{\includegraphics[width=3in]{figures/heatrefractiontcdp.png}\label{fig:tcdp}} 
306 
\subfigure[Heat Flow Depth Profile]{\includegraphics[width=3in]{figures/heatrefractionhf.png}\label{fig:hf}} 
307 
\caption{Example 5b: Depth profiles down centre of model.} 
308 
\label{figs:dps} 
309 
\end{figure} 
310 

311 
\section{Line profiles of 2D data} 
312 
\sslist{example05b.py and cblib.py} 
313 

314 
We want to investigate the profile of the data of the last example. 
315 
We are particularly interested in the depth profile of the heat flux which is 
316 
the second component of $\kappa \nabla T$. We extend the script developed in the 
317 
previous section to show how for instance vertical profile can be plotted. 
318 

319 
The first important information is that \esc assumes that $\kappa \nabla T$ is not smooth and 
320 
will in fact represent the values at numerical interpolation points. This assumption is reasonable as 
321 
the flux is the product of the piecewise constant function $\kappa$ and 
322 
the gradient of the temperature $T$ which has a kink across the rock interface. 
323 
Before plotting this function we need to smooth it using the 
324 
\verbProjector() class; 
325 
\begin{python} 
326 
from esys.escript.pdetools import Projector 
327 
proj=Projector(domain) 
328 
qu=proj(kappa*grad(T)) 
329 
\end{python} 
330 
The \verbproj object provides a mechanism to distribute values given at the numerical interpolation points 
331 
 in this case the heat flux  to the nodes of the FEM mesh. \verbqu has the same attached function space 
332 
like the temperature \verbT. The smoothed flux is interpolated 
333 
to a regular $200\times 200$ grid; 
334 
\begin{python} 
335 
xiq,yiq,ziq = toRegGrid(qu[1],200,200) 
336 
\end{python} 
337 
using the \verbtoRegGrid function from the cookbook library which we are using for the contour plot. 
338 
At return \verbziq[j,i] is the value of vertical heat flux at point 
339 
(\verbxiq[i],\verbyiq[j]). We can easily create deep profiles now by 
340 
plotting slices \verbziq[:,i] over \verbyiq. The following script 
341 
creates a deep profile at $x_{0}=\frac{width}{2}$; 
342 
\begin{python} 
343 
cut=int(len(xiq)/2) 
344 
pl.plot(ziq[:,cut]*1000.,yiq) 
345 
pl.title("Vertical Heat Flow Depth Profile") 
346 
pl.xlabel("Heat Flow (mW/m^2)") 
347 
pl.ylabel("Depth (m)") 
348 
pl.savefig(os.path.join(save_path,"hf.png")) 
349 
\end{python} 
350 
This process can be repeated for various variations of the solution. In this case we have temperature, temperature gradient, thermal conductivity and heat flow \reffig{figs:dps}. 
351 

352 
\begin{figure}[ht] 
353 
\centerline{\includegraphics[width=5.in]{figures/heatrefractionflux}} 
354 
\caption{Example 5c: Heat refraction model with gradient indicated by vectors.} 
355 
\label{fig:hr001qumodel} 
356 
\end{figure} 
357 

358 
\section{Arrow plots in \mpl} 
359 
\sslist{example05c.py and cblib.py} 
360 
We would like to visualise the distribution of the flux $\kappa \nabla T$ over the domain 
361 
and produce a plot like shown in \reffig{fig:hr001qumodel}. 
362 
The plot puts together three components. A contour plot of the temperature, 
363 
a colored representation of the two subdomains where colour represents the thermal conductivity 
364 
in the particular region and finally the arrows representing the direction of flux. 
365 

366 
Contours have already been discussed in \refSec{Sec:2DHD plot}. To show subdomains, 
367 
we need to go back to \pycad data to get the points used to describe the boundary of the 
368 
subdomains. We have created the \verbCurveLoop class object 
369 
\verbtblockloop to define the boundary of the upper subdomain. 
370 
We use the \verbgetPolygon() method of \verbCurveLoop to get 
371 
access to the \verbPoints used top define the boundary. The statement 
372 
\begin{python} 
373 
[ p.getCoordinates() for p in tblockloop.getPolygon() ]) 
374 
\end{python} 
375 
creates a list of the node coordinates of all the points in question. In order 
376 
to simplify the selection of the $x$ and $y$ coordinates the list is converted 
377 
into \modnumpy array. To add the area colored in brown to the plot we use; 
378 
\begin{python} 
379 
import pylab as pl 
380 
import numarray as np 
381 
tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ]) 
382 
pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=1000) 
383 
\end{python} 
384 
We can apply the same code to \verbbblockloop to a red area for this subdomain to the block. 
385 

386 
To plot vectors representing the flux orientation we use the 
387 
\verbquiver function in \pylab. The function places vectors at locations in the domain. 
388 
For instance one can plot vectors at the locations of the sample points used by \esc 
389 
to represent the flux \verbkappa*grad(T). As a vector is plotted at each sample point one typically ends 
390 
up with two many vectors. So one needs to select a subset of points: 
391 
First we create a coarse grid of points on a rectangular mesh, e.g. $20 \times 20$ points. Here we choose a grid of points which are located at the center of a \verbnx $\times$ \verbny grid; 
392 
\begin{python} 
393 
dx = width/nx # x spacing 
394 
dy = depth/ny # y spacing 
395 
grid = [ ] # the grid points 
396 
for j in xrange(0,ny1): 
397 
for i in xrange(0,nx1): 
398 
grid.append([dx/2+dx*i,dy/2+dy*j]) 
399 
\end{python} 
400 
With the \verbLocator \esc provides a mechanism to identify sample points that are closest 
401 
to the the grid points we have selected and to get the data at these data points; 
402 
\begin{python} 
403 
from esys.escript.pdetools import Locator 
404 
flux=kappa*grad(T) 
405 
fluxLoc = Locator(flux.getFunctionSpace(),grid) 
406 
subflux= fluxLoc(flux) 
407 
\end{python} 
408 
\verbsubflux gives now a list of flux component at certain sample points. To get the 
409 
list of the sample point coordinates one can use the \verbgetX() method of the 
410 
\verbLocator; 
411 
\begin{python} 
412 
subfluxloc = fluxLoc.getX() 
413 
\end{python} 
414 
To simplify the selection of $x$ and $y$ components it is convenient 
415 
to transform \verbsubflux and \verbsubfluxloc to \numpy arrays 
416 
\verbxflux, \verbflux. 
417 
This function is implemented in the \verbsubsample 
418 
in the \file{clib.py} file so we can use it in other examples. One can easily use this function 
419 
to create a vector plot of the flux; 
420 
\begin{python} 
421 
from cblib import subsample 
422 
xflux, flux=subsample(kappa*grad(T), nx=20, ny=20) 
423 
pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white") 
424 
\end{python} 
425 
We add title and labels; 
426 
\begin{python} 
427 
pl.title("Heat Refraction across a clinal structure.") 
428 
pl.xlabel("Horizontal Displacement (m)") 
429 
pl.ylabel("Depth (m)") 
430 
pl.title("Heat Refraction across a clinal structure \n with gradient quivers.") 
431 
pl.savefig(os.path.join(saved_path,"flux.png")) 
432 
\end{python} 
433 
to get the desired result, see \reffig{fig:hr001qumodel}. 
434 

435 
\begin{figure}[ht] 
436 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction2flux}} 
437 
\caption{Example 6: Heat refraction model with three blocks and heat flux.} 
438 
\label{fig:hr002qumodel} 
439 
\end{figure} 
440 

441 
\section{Example 6:Fault and Overburden Model} 
442 
\sslist{example06.py and cblib.py} 
443 
A slightly more complicated model can be found in the examples \textit{heatrefraction2_solver.py} where three blocks are used within the model, see~\reffig{fig:hr002qumodel}. It is left to the reader to work through this example. 
444 

445 
