32 
blocks have the same temperature. 
blocks have the same temperature. 
33 


34 
\subsection{1D Heat Diffusion Equation} 
\subsection{1D Heat Diffusion Equation} 
35 
We can model the heat distribution of this problem over time using one 
We can model the heat distribution of this problem over time using the one 
36 
dimensional heat diffusion equation\footnote{A detailed discussion on how the 
dimensional heat diffusion equation\footnote{A detailed discussion on how the 
37 
heat diffusion equation is derived can be found at 
heat diffusion equation is derived can be found at 
38 
\url{ 
\url{ 
58 
change in temperature with time while $\frac{\partial ^2 T}{\partial x^2}$ is 
change in temperature with time while $\frac{\partial ^2 T}{\partial x^2}$ is 
59 
the spatial change of temperature. As there is only a single spatial dimension 
the spatial change of temperature. As there is only a single spatial dimension 
60 
to our problem, our temperature solution $T$ is only dependent on the time $t$ 
to our problem, our temperature solution $T$ is only dependent on the time $t$ 
61 
and our signed distance from the the blockblock interface $x$. 
and our signed distance from the blockblock interface $x$. 
62 


63 
\subsection{PDEs and the General Form} 
\subsection{PDEs and the General Form} 
64 
It is possible to solve PDE \refEq{eqn:hd} analytically and obtain an exact 
It is possible to solve PDE \refEq{eqn:hd} analytically and obtain an exact 
75 
steady linear PDE which involves spatial derivatives only and needs to be solved 
steady linear PDE which involves spatial derivatives only and needs to be solved 
76 
in each time step to progress in time. \esc can help us here. 
in each time step to progress in time. \esc can help us here. 
77 


78 
For time discretization we use the Backwards Euler approximation 
For time discretisation we use the Backward Euler approximation 
79 
scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is 
scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is 
80 
based on the 
based on the approximation 

approximation 

81 
\begin{equation} 
\begin{equation} 
82 
\frac{\partial T(t)}{\partial t} \approx \frac{T(t)T(th)}{h} 
\frac{\partial T(t)}{\partial t} \approx \frac{T(t)T(th)}{h} 
83 
\label{eqn:beuler} 
\label{eqn:beuler} 
106 
Notice that we evaluate the spatial derivative term at the current time 
Notice that we evaluate the spatial derivative term at the current time 
107 
$t^{(n)}$  therefore the name \textbf{backward Euler} scheme. Alternatively, 
$t^{(n)}$  therefore the name \textbf{backward Euler} scheme. Alternatively, 
108 
one can evaluate the spatial derivative term at the previous time $t^{(n1)}$. 
one can evaluate the spatial derivative term at the previous time $t^{(n1)}$. 
109 
This 
This approach is called the \textbf{forward Euler} scheme. This scheme can 
110 
approach is called the \textbf{forward Euler} scheme. This scheme can provide 
provide some computational advantages, which 

some computational advantages, which 

111 
are not discussed here. However, the \textbf{forward Euler} scheme has a major 
are not discussed here. However, the \textbf{forward Euler} scheme has a major 
112 
disadvantage. Namely, depending on the 
disadvantage. Namely, depending on the 
113 
material parameters as well as the domain discretization of the spatial 
material parameters as well as the domain discretization of the spatial 
114 
derivative term, the time step size $h$ needs to be chosen sufficiently small to 
derivative term, the time step size $h$ needs to be chosen sufficiently small to 
115 
achieve a stable temperature when progressing in time. Stabiliy is achieved if 
achieve a stable temperature when progressing in time. Stability is achieved if 
116 
the temperature does not grow beyond its initial bounds and become 
the temperature does not grow beyond its initial bounds and becomes 
117 
nonphysical. 
nonphysical. 
118 
The backward Euler scheme, which we use here, is unconditionally stable meaning 
The backward Euler scheme, which we use here, is unconditionally stable meaning 
119 
that under the assumption of a 
that under the assumption of a 
120 
physically correct problem setup the temperature approximation remains physical 
physically correct problem setup the temperature approximation remains physical 
121 
for all time steps. 
for all time steps. 
122 
The user needs to keep in mind that the discretization error introduced by 
The user needs to keep in mind that the discretisation error introduced by 
123 
\refEq{eqn:beuler} 
\refEq{eqn:beuler} 
124 
is sufficiently small, thus a good approximation of the true temperature is 
is sufficiently small, thus a good approximation of the true temperature is 
125 
computed. It is 
computed. It is therefore very important that any results are viewed with 
126 
therefore very important that any results are viewed with caution. For example, 
caution. For example, one may compare the results for different time and 
127 
one may compare the results for different time and spatial step sizes. 
spatial step sizes. 
128 


129 
To get the temperature $T^{(n)}$ at time $t^{(n)}$ we need to solve the linear 
To get the temperature $T^{(n)}$ at time $t^{(n)}$ we need to solve the linear 
130 
differential equation \refEq{eqn:hddisc} which only includes spatial 
differential equation \refEq{eqn:hddisc} which only includes spatial 
131 
derivatives. To solve this problem 
derivatives. To solve this problem we want to use \esc. 

we want to use \esc. 

132 


133 
In \esc any given PDE can be described by the general form. For the purpose of 
In \esc any given PDE can be described by the general form. For the purpose of 
134 
this introduction we illustrate a simpler version of the general form for full 
this introduction we illustrate a simpler version of the general form for full 
185 
A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to 
A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to 
186 
prescribe a known value to the unknown solution (in our example the temperature) 
prescribe a known value to the unknown solution (in our example the temperature) 
187 
on parts of the boundary or on the entire boundary of the region of interest. 
on parts of the boundary or on the entire boundary of the region of interest. 
188 
We discuss Dirichlet boundary condition in our second example presented in 
We discuss the Dirichlet boundary condition in our second example presented in 
189 
Section~\ref{Sec:1DHDv0}. 
Section~\ref{Sec:1DHDv0}. 
190 


191 
However, for this example we have made the model assumption that the system is 
However, for this example we have made the model assumption that the system is 
192 
insulated, so we need 
insulated, so we need to add an appropriate boundary condition to prevent 

to add an appropriate boundary condition to prevent 

193 
any loss or inflow of energy at the boundary of our domain. Mathematically this 
any loss or inflow of energy at the boundary of our domain. Mathematically this 
194 
is expressed by prescribing 
is expressed by prescribing 
195 
the heat flux $\kappa \frac{\partial T}{\partial x}$ to zero. In our simplified 
the heat flux $\kappa \frac{\partial T}{\partial x}$ to zero. In our simplified 
204 
\end{equation} 
\end{equation} 
205 
where $n$ is the outer normal field \index{outer normal field} at the surface 
where $n$ is the outer normal field \index{outer normal field} at the surface 
206 
of the domain. 
of the domain. 
207 
The $\cdot$ (dot) refers to the dot product of the vectors $\nabla T$ and $n$. 
The $\cdot$ (dot) refers to the dot product of the vectors $\nabla T$ and $n$. 
208 
In fact, the term $\nabla T \cdot n$ is the normal derivative of 
In fact, the term $\nabla T \cdot n$ is the normal derivative of 
209 
the temperature $T$. Other notations used here are\footnote{The \esc notation 
the temperature $T$. Other notations used here are\footnote{The \esc notation 
210 
for the normal 
for the normal 
212 
\begin{equation} 
\begin{equation} 
213 
\nabla T \cdot n = \frac{\partial T}{\partial n} \; . 
\nabla T \cdot n = \frac{\partial T}{\partial n} \; . 
214 
\end{equation} 
\end{equation} 
215 
A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neuman boundary 
A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neumann boundary 
216 
condition} for the PDE. 
condition} for the PDE. 
217 


218 
The PDE \refEq{eqn:hdgenf} 
The PDE \refEq{eqn:hdgenf} 
219 
and the Neuman boundary condition~\ref{eqn:hdgenf} (potentially together with 
and the Neumann boundary condition~\ref{eqn:hdgenf} (potentially together with 
220 
the Dirichlet boundary conditions) define a \textbf{boundary value problem}. 
the Dirichlet boundary conditions) define a \textbf{boundary value problem}. 
221 
It is the nature of a boundary value problem to allow making statements about 
It is the nature of a boundary value problem to allow making statements about 
222 
the solution in the 
the solution in the 
223 
interior of the domain from information known on the boundary only. In most 
interior of the domain from information known on the boundary only. In most 
224 
cases 
cases we use the term partial differential equation but in fact it is a 
225 
we use the term partial differential equation but in fact it is a boundary value 
boundary value problem. 

problem. 

226 
It is important to keep in mind that boundary conditions need to be complete and 
It is important to keep in mind that boundary conditions need to be complete and 
227 
consistent in the sense that 
consistent in the sense that 
228 
at any point on the boundary either a Dirichlet or a Neuman boundary condition 
at any point on the boundary either a Dirichlet or a Neumann boundary condition 
229 
must be set. 
must be set. 
230 


231 
Conveniently, \esc makes default assumption on the boundary conditions which the 
Conveniently, \esc makes a default assumption on the boundary conditions which 
232 
user may modify where appropriate. 
the user may modify where appropriate. 
233 
For a problem of the form in~\refEq{eqn:commonform nabla} the default 
For a problem of the form in~\refEq{eqn:commonform nabla} the default 
234 
condition\footnote{In the form of the \esc users guide which is using the 
condition\footnote{In the \esc user guide which uses the Einstein convention 
235 
Einstein convention is written as 
this is written as 
236 
$n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is; 
$n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is; 
237 
\begin{equation}\label{NEUMAN 2} 
\begin{equation}\label{NEUMAN 2} 
238 
n\cdot A \cdot\nabla u = 0 
n\cdot A \cdot\nabla u = 0 
247 
\begin{equation}\label{NEUMAN 2b} 
\begin{equation}\label{NEUMAN 2b} 
248 
\kappa \frac{\partial T}{\partial x} = 0 
\kappa \frac{\partial T}{\partial x} = 0 
249 
\end{equation} 
\end{equation} 
250 
for the boundary of the domain. This is identical to the Neuman boundary 
for the boundary of the domain. This is identical to the Neumann boundary 
251 
condition we want to set. \esc will take care of this condition for us. We 
condition we want to set. \esc will take care of this condition for us. We 
252 
discuss the Dirichlet boundary condition later. 
discuss the Dirichlet boundary condition later. 
253 


263 
each time step to get the updated temperature. Thirdly, we need to define the 
each time step to get the updated temperature. Thirdly, we need to define the 
264 
coefficients of the PDE and finally we need to solve the PDE. The last two steps 
coefficients of the PDE and finally we need to solve the PDE. The last two steps 
265 
need to be repeated until the final time marker has been reached. The work flow 
need to be repeated until the final time marker has been reached. The work flow 
266 
described in \reffig{fig:wf}. 
is described in \reffig{fig:wf}. 
267 
% \begin{enumerate} 
% \begin{enumerate} 
268 
% \item create domain 
% \item create domain 
269 
% \item create PDE 
% \item create PDE 
279 
\begin{figure}[h!] 
\begin{figure}[h!] 
280 
\centering 
\centering 
281 
\includegraphics[width=1in]{figures/workflow.png} 
\includegraphics[width=1in]{figures/workflow.png} 
282 
\caption{Workflow for developing an \esc model and solution.} 
\caption{Workflow for developing an \esc model and solution} 
283 
\label{fig:wf} 
\label{fig:wf} 
284 
\end{figure} 
\end{figure} 
285 


305 
\subsection{The Domain Constructor in \esc} 
\subsection{The Domain Constructor in \esc} 
306 
\label{ss:domcon} 
\label{ss:domcon} 
307 
Whilst it is not strictly relevant or necessary, a better understanding of 
Whilst it is not strictly relevant or necessary, a better understanding of 
308 
how values are spatially distributed (\textit{e.g.} Temperature) and how PDE 
how values are spatially distributed (\textit{e.g.} Temperature) and how PDE 
309 
coefficients are interpreted in \esc can be helpful. 
coefficients are interpreted in \esc can be helpful. 
310 


311 
There are various ways to construct domain objects. The simplest form is a 
There are various ways to construct domain objects. The simplest form is a 
319 
\url{http://en.wikipedia.org/wiki/Finite_element_method} for details.}. 
\url{http://en.wikipedia.org/wiki/Finite_element_method} for details.}. 
320 
Therefore we will have access to an approximation of the true PDE solution 
Therefore we will have access to an approximation of the true PDE solution 
321 
only. 
only. 
322 
The quality of the approximation depends  besides other factors mainly on the 
The quality of the approximation depends  besides other factors  mainly on the 
323 
number of elements being used. In fact, the 
number of elements being used. In fact, the 
324 
approximation becomes better when more elements are used. However, computational 
approximation becomes better when more elements are used. However, computational 
325 
cost grows with the number of 
cost grows with the number of 
330 
elements. 
elements. 
331 
As shown in \reffig{fig:fs}, an element is defined by the nodes that are used to 
As shown in \reffig{fig:fs}, an element is defined by the nodes that are used to 
332 
describe its vertices. 
describe its vertices. 
333 
To represent spatial distributed values the user can use 
To represent spatially distributed values the user can use 
334 
the values at the nodes, at the elements in the interior of the domain or at the 
the values at the nodes, at the elements in the interior of the domain or at the 
335 
elements located at the surface of the domain. 
elements located on the surface of the domain. 
336 
The different approach used to represent values is called \textbf{function 
The different approach used to represent values is called \textbf{function 
337 
space} and is attached to all objects 
space} and is attached to all objects 
338 
in \esc representing a spatial distributed value such as the solution of a PDE. 
in \esc representing a spatially distributed value such as the solution of 
339 
The three 
a PDE. The three function spaces we use at the moment are; 

function spaces we use at the moment are; 

340 
\begin{enumerate} 
\begin{enumerate} 
341 
\item the nodes, called by \verbContinuousFunction(domain) ; 
\item the nodes, called by \verbContinuousFunction(domain) ; 
342 
\item the elements/cells, called by \verbFunction(domain) ; and 
\item the elements/cells, called by \verbFunction(domain) ; and 
343 
\item the boundary, called by \verbFunctionOnBoundary(domain) . 
\item the boundary, called by \verbFunctionOnBoundary(domain). 
344 
\end{enumerate} 
\end{enumerate} 
345 
A function space object such as \verbContinuousFunction(domain) has the method 
A function space object such as \verbContinuousFunction(domain) has the method 
346 
\verbgetX attached to it. This method returns the 
\verbgetX attached to it. This method returns the 
348 
particular function space. So the 
particular function space. So the 
349 
call \verbContinuousFunction(domain).getX() will return the coordinates of the 
call \verbContinuousFunction(domain).getX() will return the coordinates of the 
350 
nodes used to describe the domain while 
nodes used to describe the domain while 
351 
the \verbFunction(domain).getX() returns the coordinates of numerical 
\verbFunction(domain).getX() returns the coordinates of numerical 
352 
integration points within elements, see 
integration points within elements, see \reffig{fig:fs}. 

\reffig{fig:fs}. 

353 


354 
This distinction between different representations of spatially distributed 
This distinction between different representations of spatially distributed 
355 
values 
values 
359 
\verbFunction() type. 
\verbFunction() type. 
360 
On the other hand a temperature distribution must be continuous and needs to be 
On the other hand a temperature distribution must be continuous and needs to be 
361 
represented with a \verbContinuousFunction() function space. 
represented with a \verbContinuousFunction() function space. 
362 
An influx may only be defined at the boundary and is therefore a \verb 
An influx may only be defined at the boundary and is therefore a 
363 
FunctionOnBoundary() object. 
\verbFunctionOnBoundary() object. 
364 
\esc allows certain transformations of the function spaces. A \verb 
\esc allows certain transformations of the function spaces. A 
365 
ContinuousFunction() can be transformed into a \verbFunctionOnBoundary() 
\verbContinuousFunction() can be transformed into a 
366 
or \verbFunction(). On the other hand there is not enough information in a 
\verbFunctionOnBoundary() or \verbFunction(). On the other hand there is 
367 
\verb FunctionOnBoundary() to transform it to a \verb ContinuousFunction() . 
not enough information in a \verbFunctionOnBoundary() to transform it to a 
368 

\verbContinuousFunction(). 
369 
These transformations, which are called \textbf{interpolation} are invoked 
These transformations, which are called \textbf{interpolation} are invoked 
370 
automatically by \esc if needed. 
automatically by \esc if needed. 
371 


372 
Later in this introduction we discuss how 
Later in this introduction we discuss how 
373 
to define specific areas of geometry with different materials which are 
to define specific areas of geometry with different materials which are 
374 
represented by different material coefficients such the 
represented by different material coefficients such as the 
375 
thermal conductivities $\kappa$. A very powerful technique to define these types 
thermal conductivities $\kappa$. A very powerful technique to define these types 
376 
of PDE 
of PDE 
377 
coefficients is tagging. Blocks of materials and boundaries can be named and 
coefficients is tagging. Blocks of materials and boundaries can be named and 
384 
\subsection{A Clarification for the 1D Case} 
\subsection{A Clarification for the 1D Case} 
385 
\label{SEC: 1D CLARIFICATION} 
\label{SEC: 1D CLARIFICATION} 
386 
It is necessary for clarification that we revisit our general PDE from 
It is necessary for clarification that we revisit our general PDE from 
387 
\refeq{eqn:commonform nabla} for two dimensional domain. \esc is inherently 
\refeq{eqn:commonform nabla} for a two dimensional domain. \esc is inherently 
388 
designed to solve problems that are greater than one dimension and so 
designed to solve problems that are multidimensional and so 
389 
\refEq{eqn:commonform nabla} needs to be read as a higher dimensional problem. 
\refEq{eqn:commonform nabla} needs to be read as a higher dimensional problem. 
390 
In the case of two spatial dimensions the \textit{Nabla operator} has in fact 
In the case of two spatial dimensions the \textit{Nabla operator} has in fact 
391 
two components $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial 
two components $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial 
434 
if the functions and classes required are known they can be specified 
if the functions and classes required are known they can be specified 
435 
individually. The function \verbLinearPDE has been imported explicitly for 
individually. The function \verbLinearPDE has been imported explicitly for 
436 
ease of use later in the script. \verbRectangle is going to be our type of 
ease of use later in the script. \verbRectangle is going to be our type of 
437 
domain. The module \verb unitsSI provides support for SI unit definitions with 
domain. The module \verbunitsSI provides support for SI unit definitions with 
438 
our variables. 
our variables. 
439 


440 
Once our library dependencies have been established, defining the problem 
Once our library dependencies have been established, defining the problem 
460 
the number of calculations that must be solved per time step. This means more 
the number of calculations that must be solved per time step. This means more 
461 
computational time is required to produce a solution. In this \textit{1D} 
computational time is required to produce a solution. In this \textit{1D} 
462 
problem, the bar is defined as being 1 metre long. An appropriate step size 
problem, the bar is defined as being 1 metre long. An appropriate step size 
463 
\verbndx would be 1 to 10\% of the length. Our \verbndy need only be 1, this 
\verbndx would be 1 to 10\% of the length. Our \verbndy needs only be 1, 
464 
is because our problem stipulates no partial derivatives in the $y$ direction. 
this is because our problem stipulates no partial derivatives in the $y$ 
465 

direction. 
466 
Thus the temperature does not vary with $y$. Hence, the model parameters can be 
Thus the temperature does not vary with $y$. Hence, the model parameters can be 
467 
defined as follows; note we have used the \verb unitsSI convention to make sure 
defined as follows; note we have used the \verbunitsSI convention to make sure 
468 
all our input units are converted to SI. 
all our input units are converted to SI. 
469 
\begin{python} 
\begin{python} 
470 
mx = 500.*m #meters  model length 
mx = 500.*m #meters  model length 
497 
i=0 #loop counter 
i=0 #loop counter 
498 
\end{python} 
\end{python} 
499 
Now that we know our inputs we will build a domain using the 
Now that we know our inputs we will build a domain using the 
500 
\verb Rectangle() function from \verb finley . The four arguments allow us to define our domain 
\verbRectangle() function from \FINLEY. The four arguments allow us to 
501 
\verb model as: 
define our domain \verbmodel as: 
502 
\begin{python} 
\begin{python} 
503 
#generate domain using rectangle 
#generate domain using rectangle 
504 
blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
505 
\end{python} 
\end{python} 
506 
\verb blocks now describes a domain in the manner of Section \ref{ss:domcon}. 
\verbblocks now describes a domain in the manner of Section \ref{ss:domcon}. 
507 


508 
With a domain and all our required variables established, it is now possible to 
With a domain and all the required variables established, it is now possible to 
509 
set up our PDE so that it can be solved by \esc. The first step is to define the 
set up our PDE so that it can be solved by \esc. The first step is to define the 
510 
type of PDE that we are trying to solve in each time step. In this example it is 
type of PDE that we are trying to solve in each time step. In this example it is 
511 
a single linear PDE\footnote{in comparison to a system of PDEs which we discuss 
a single linear PDE\footnote{in contrast to a system of PDEs which we discuss 
512 
later.}. We also need to state the values of our general form variables. 
later.}. We also need to state the values of our general form variables. 
513 
\begin{python} 
\begin{python} 
514 
mypde=LinearPDE(blocks) 
mypde=LinearPDE(blocks) 
516 
A[0,0]=kappa 
A[0,0]=kappa 
517 
mypde.setValue(A=A, D=rhocp/h) 
mypde.setValue(A=A, D=rhocp/h) 
518 
\end{python} 
\end{python} 
519 
In a many cases it may be possible to decrease the computational time of the 
In many cases it may be possible to decrease the computational time of the 
520 
solver if the PDE is symmetric. 
solver if the PDE is symmetric. 
521 
Symmetry of a PDE is defined by; 
Symmetry of a PDE is defined by; 
522 
\begin{equation}\label{eqn:symm} 
\begin{equation}\label{eqn:symm} 
524 
\end{equation} 
\end{equation} 
525 
Symmetry is only dependent on the $A$ coefficient in the general form and the 
Symmetry is only dependent on the $A$ coefficient in the general form and the 
526 
other coefficients $D$ as well as the right hand side $Y$. From the above 
other coefficients $D$ as well as the right hand side $Y$. From the above 
527 
definition we can see that our PDE is symmetric. The \verb LinearPDE class 
definition we can see that our PDE is symmetric. The \verbLinearPDE class 
528 
provides the method \method{checkSymmetry} to check if the given PDE is 
provides the method \method{checkSymmetry} to check if the given PDE is 
529 
symmetric. As our PDE is symmetrical we enable symmetry via; 
symmetric. As our PDE is symmetrical we enable symmetry via; 
530 
\begin{python} 
\begin{python} 
531 
myPDE.setSymmetryOn() 
myPDE.setSymmetryOn() 
532 
\end{python} 
\end{python} 
533 
Next we need to establish the initial temperature distribution \verbT. We need 
Next we need to establish the initial temperature distribution \verbT. We need 
534 
to 
to 
535 
assign the value \verbT1 to all sample points left to the contact interface at 
assign the value \verbT1 to all sample points left to the contact interface at 
536 
$x\hackscore{0}=\frac{mx}{2}$ 
$x\hackscore{0}=\frac{mx}{2}$ 
537 
and the value \verbT2 right to the contact interface. \esc 
and the value \verbT2 right to the contact interface. \esc 
538 
provides the \verbwhereNegative function to construct this. In fact, 
provides the \verbwhereNegative function to construct this. More 
539 
\verbwhereNegative returns the value $1$ at those sample points where the 
specifically, \verbwhereNegative returns the value $1$ at those sample points 
540 
argument 
where the argument has a negative value. Otherwise zero is returned. 
541 
has a negative value. Otherwise zero is returned. If \verbx are the 
If \verbx are the $x\hackscore{0}$ 

$x\hackscore{0}$ 

542 
coordinates of the sample points used to represent the temperature distribution 
coordinates of the sample points used to represent the temperature distribution 
543 
then \verbx[0]boundloc gives us a negative value for 
then \verbx[0]boundloc gives us a negative value for 
544 
all sample points left to the interface and nonnegative value to 
all sample points left to the interface and nonnegative value to 
548 
T= T1*whereNegative(x[0]boundloc)+T2*(1whereNegative(x[0]boundloc)) 
T= T1*whereNegative(x[0]boundloc)+T2*(1whereNegative(x[0]boundloc)) 
549 
\end{python} 
\end{python} 
550 
we get the desired temperature distribution. To get the actual sample points 
we get the desired temperature distribution. To get the actual sample points 
551 
\verbx we use 
\verbx we use the \verbgetX() method of the function space 
552 
the \verbgetX() method of the function space \verbSolution(blocks) 
\verbSolution(blocks) which is used to represent the solution of a PDE; 

which is used to represent the solution of a PDE; 

553 
\begin{python} 
\begin{python} 
554 
x=Solution(blocks).getX() 
x=Solution(blocks).getX() 
555 
\end{python} 
\end{python} 
558 
the initial temperature \verbT is using these sample points for 
the initial temperature \verbT is using these sample points for 
559 
representation. 
representation. 
560 
Although \esc is trying to be forgiving with the choice of sample points and to 
Although \esc is trying to be forgiving with the choice of sample points and to 
561 
convert 
convert 
562 
where necessary the adjustment of the function space is not always possible. So 
where necessary the adjustment of the function space is not always possible. So 
563 
it is 
it is advisable to make a careful choice on the function space used. 

advisable to make a careful choice on the function space used. 

564 


565 
Finally we initialise an iteration loop to solve our PDE for all the time steps 
Finally we initialise an iteration loop to solve our PDE for all the time steps 
566 
we specified in the variable section. As the right hand side of the general form 
we specified in the variable section. As the right hand side of the general form 
576 
totE = integrate(rhocp*T) #get the total heat (energy) in the system 
totE = integrate(rhocp*T) #get the total heat (energy) in the system 
577 
\end{python} 
\end{python} 
578 
The last statement in this script calculates the total energy in the system as 
The last statement in this script calculates the total energy in the system as 
579 
volume integral 
the volume integral of $\rho c\hackscore{p} T$ over the block. 
580 
of $\rho c\hackscore{p} T$ over the block. As the blocks are insulated no energy 
As the blocks are insulated no energy should be lost or added. 

should be get lost or added. 

581 
The total energy should stay constant for the example discussed here. 
The total energy should stay constant for the example discussed here. 
582 


583 
\subsection{Running the Script} 
\subsection{Running the Script} 
584 
The script presented so far is available under 
The script presented so far is available under 
585 
\verbexample01a.py. You can edit this file with your favourite text editor. 
\verbexample01a.py. You can edit this file with your favourite text editor. 
586 
On most operating systems\footnote{The you can use \texttt{runescript} launcher 
On most operating systems\footnote{The \texttt{runescript} launcher is not 
587 
is not supported under {\it MS Windows} yet.} you can use the 
supported under {\it MS Windows} yet.} you can use the 
588 
\program{runescript} command 
\program{runescript} command 
589 
to launch {\it escript} scripts. For the example script use; 
to launch {\it escript} scripts. For the example script use; 
590 
\begin{verbatim} 
\begin{verbatim} 
595 
\begin{verbatim} 
\begin{verbatim} 
596 
python example01a.py 
python example01a.py 
597 
\end{verbatim} 
\end{verbatim} 
598 
if the system is configured correctly (Please talk to your system 
if the system is configured correctly (please talk to your system 
599 
administrator). 
administrator). 
600 


601 
\begin{figure} 
\begin{figure} 
602 
\begin{center} 
\begin{center} 
603 
\includegraphics[width=4in]{figures/ttblockspyplot150} 
\includegraphics[width=4in]{figures/ttblockspyplot150} 
604 
\caption{Example 1b: Total Energy in the Blocks over Time (in seconds).} 
\caption{Example 1b: Total Energy in the Blocks over Time (in seconds)} 
605 
\label{fig:onedheatout1} 
\label{fig:onedheatout1} 
606 
\end{center} 
\end{center} 
607 
\end{figure} 
\end{figure} 
613 
use a variety of free \pyt packages for visualisation. 
use a variety of free \pyt packages for visualisation. 
614 
Two types will be demonstrated in this cookbook; 
Two types will be demonstrated in this cookbook; 
615 
\mpl\footnote{\url{http://matplotlib.sourceforge.net/}} and 
\mpl\footnote{\url{http://matplotlib.sourceforge.net/}} and 
616 
\verb VTK \footnote{\url{http://www.vtk.org/}} visualisation. 
\verbVTK\footnote{\url{http://www.vtk.org/}}. 
617 
The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} 
The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} 
618 
and is good for basic graphs and plots. 
and is good for basic graphs and plots. 
619 
For more complex visualisation tasks in particular, two and three dimensional 
For more complex visualisation tasks, in particular two and three dimensional 
620 
problems we recommend the use of more advanced tools. For instance, \mayavi 
problems we recommend the use of more advanced tools. For instance, \mayavi 
621 
\footnote{\url{http://code.enthought.com/projects/mayavi/}} 
\footnote{\url{http://code.enthought.com/projects/mayavi/}} 
622 
which is based upon the \verbVTK toolkit. The usage of \verbVTK based 
which is based upon the \verbVTK toolkit. The usage of \verbVTK based 
623 
visualization is discussed in Chapter~\ref{Sec:2DHD} which focusses on a two 
visualisation is discussed in Chapter~\ref{Sec:2DHD} which focuses on a two 
624 
dimensional PDE. 
dimensional PDE. 
625 


626 
For our simple granite block problem, we have two plotting tasks. Firstly, we 
For our simple granite block problem, we have two plotting tasks. Firstly, we 
627 
are interested in showing the 
are interested in showing the 
628 
behaviour of the total energy over time and secondly, how the temperature 
behaviour of the total energy over time and secondly, how the temperature 
629 
distribution within the block is 
distribution within the block is developing over time. 
630 
developing over time. Lets start with the first task. 
Let us start with the first task. 
631 


632 
The trick is to create a record of the time marks and the corresponding total 
The idea is to create a record of the time marks and the corresponding total 
633 
energies observed. 
energies observed. 
634 
\pyt provides the concept of lists for this. Before 
\pyt provides the concept of lists for this. Before 
635 
the time loop is opened we create empty lists for the time marks \verbt_list 
the time loop is opened we create empty lists for the time marks \verbt_list 
651 
E_list.append(totE) # add current total energy to record 
E_list.append(totE) # add current total energy to record 
652 
\end{python} 
\end{python} 
653 
To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs 
To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs 
654 
to be loaded before used; 
to be loaded before use; 
655 
\begin{python} 
\begin{python} 
656 
import pylab as pl # plotting package. 
import pylab as pl # plotting package. 
657 
\end{python} 
\end{python} 
658 
Here we are not using the \verbfrom pylab import * in order to avoid name 
Here we are not using \verbfrom pylab import * in order to avoid name 
659 
clashes for function names 
clashes for function names within \esc. 

within \esc. 

660 


661 
The following statements are added to the script after the time loop has been 
The following statements are added to the script after the time loop has been 
662 
completed; 
completed; 
667 
pl.savefig("totE.png") 
pl.savefig("totE.png") 
668 
\end{python} 
\end{python} 
669 
The first statement hands over the time marks and corresponding total energies 
The first statement hands over the time marks and corresponding total energies 
670 
to the plotter. 
to the plotter. 
671 
The second statment is setting the title for the plot. The third statement 
The second statement sets the title for the plot. The third statement 
672 
sets the axis ranges. In most cases these are set appropriately by the plotter. 
sets the axis ranges. In most cases these are set appropriately by the plotter. 
673 


674 
The last statement renders the plot and writes the 
The last statement generates the plot and writes the result into the file 
675 
result into the file \verbtotE.png which can be displayed by (almost) any 
\verbtotE.png which can be displayed by (almost) any image viewer. 

image viewer. 

676 
As expected the total energy is constant over time, see 
As expected the total energy is constant over time, see 
677 
\reffig{fig:onedheatout1}. 
\reffig{fig:onedheatout1}. 
678 


680 
\label{sec: plot T} 
\label{sec: plot T} 
681 
\sslist{example01c.py} 
\sslist{example01c.py} 
682 
For plotting the spatial distribution of the temperature we need to modify the 
For plotting the spatial distribution of the temperature we need to modify the 
683 
strategy we have used 
strategy we have used for the total energy. 
684 
for the total energy. Instead of producing a final plot at the end we will 
Instead of producing a final plot at the end we will generate a 

generate a 

685 
picture at each time step which can be browsed as a slide show or composed into 
picture at each time step which can be browsed as a slide show or composed into 
686 
a movie. 
a movie. 
687 
The first problem we encounter is that if we produce an image at each time step 
The first problem we encounter is that if we produce an image at each time step 
688 
we need 
we need to make sure that the images previously generated are not overwritten. 

to make sure that the images previously generated are not overwritten. 

689 


690 
To develop an incrementing file name we can use the following convention. It is 
To develop an incrementing file name we can use the following convention. It is 
691 
convenient to 
convenient to put all image files showing the same variable  in our case the 
692 
put all image file showing the same variable  in our case the temperature 
temperature distribution  into a separate directory. 
693 
distribution  
As part of the \verbos module\footnote{The \texttt{os} module provides 

into a separate directory. As part of the \verbos module\footnote{The 


\texttt{os} module provides 

694 
a powerful interface to interact with the operating system, see 
a powerful interface to interact with the operating system, see 
695 
\url{http://docs.python.org/library/os.html}.} \pyt 
\url{http://docs.python.org/library/os.html}.} \pyt 
696 
provides the \verbos.path.join command to build file and 
provides the \verbos.path.join command to build file and directory names in a 
697 
directory names in a platform independent way. Assuming that 
platform independent way. Assuming that 
698 
\verbsave_path is name of directory we want to put the results the command 
\verbsave_path is the name of the directory we want to put the results in the 
699 
is; 
command is; 
700 
\begin{python} 
\begin{python} 
701 
import os 
import os 
702 
os.path.join(save_path, "tempT%03d.png"%i ) 
os.path.join(save_path, "tempT%03d.png"%i ) 
703 
\end{python} 
\end{python} 
704 
where \verbi is the time step counter. 
where \verbi is the time step counter. 
705 
There are two arguments to the \verb join command. The 
There are two arguments to the \verbjoin command. The \verbsave_path 
706 
\verb save_path variable is a predefined string pointing to the directory we want to save our 
variable is a predefined string pointing to the directory we want to save our 
707 
data, for example a single subfolder called \verb data would be defined by; 
data, for example a single subfolder called \verbdata would be defined by; 
708 
\begin{verbatim} 
\begin{verbatim} 
709 
save_path = "data" 
save_path = "data" 
710 
\end{verbatim} 
\end{verbatim} 
711 
while a subfolder of \verb data called \verb example01 would be defined by; 
while a subfolder of \verbdata called \verbexample01 would be defined by; 
712 
\begin{verbatim} 
\begin{verbatim} 
713 
save_path = os.path.join("data","example01") 
save_path = os.path.join("data","example01") 
714 
\end{verbatim} 
\end{verbatim} 
715 
The second argument of \verb join \xspace contains a string which is the file 
The second argument of \verbjoin contains a string which is the file 
716 
name or subdirectory name. We can use the operator \verb% to use the value of 
name or subdirectory name. We can use the operator \verb% to use the value of 
717 
\verbi as part of our filename. The substring \verb %03d indicates that we 
\verbi as part of our filename. The substring \verb%03d indicates that we 
718 
want to substitude a value into the name; 
want to substitute a value into the name; 
719 
\begin{itemize} 
\begin{itemize} 
720 
\item \verb 0 means that small numbers should have leading zeroes; 
\item \verb 0 means that small numbers should have leading zeroes; 
721 
\item \verb 3 means that numbers should be written using at least 3 digits; 
\item \verb 3 means that numbers should be written using at least 3 digits; 
722 
and 
and 
723 
\item \verb d means that the value to substitute will be an integer. 
\item \verb d means that the value to substitute will be a decimal integer. 
724 
\end{itemize} 
\end{itemize} 
725 
To actually substitute the value of \verbi into the name write \verb %i after 
To actually substitute the value of \verbi into the name write \verb%i after 
726 
the string. 
the string. 
727 
When done correctly, the output files from this command would be place in the 
When done correctly, the output files from this command will be placed in the 
728 
directory defined by \verb save_path as; 
directory defined by \verb save_path as; 
729 
\begin{verbatim} 
\begin{verbatim} 
730 
blockspyplot001.png 
blockspyplot001.png 
738 
\begin{verbatim} 
\begin{verbatim} 
739 
mkDir(save_path) 
mkDir(save_path) 
740 
\end{verbatim} 
\end{verbatim} 
741 
will check for the existence of \verb save_path and if missing, make the 
will check for the existence of \verb save_path and if missing, create the 
742 
required directories. 
required directories. 
743 


744 
We start by modifying our solution script. 
We start by modifying our solution script. 
745 
Prior to the \verbwhile loop we will need to extract our finite solution 
Prior to the \verbwhile loop we need to extract our finite solution 
746 
points to a data object that is compatible with \mpl. First we create the node 
points to a data object that is compatible with \mpl. First we create the node 
747 
coordinates of the sample points used to represent 
coordinates of the sample points used to represent 
748 
the temperature as a \pyt list of tuples or a \numpy array as requested by the 
the temperature as a \pyt list of tuples or a \numpy array as requested by the 
765 
\includegraphics[width=4in]{figures/blockspyplot050} 
\includegraphics[width=4in]{figures/blockspyplot050} 
766 
\includegraphics[width=4in]{figures/blockspyplot200} 
\includegraphics[width=4in]{figures/blockspyplot200} 
767 
\caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps 
\caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps 
768 
$1$, $50$ and $200$.} 
$1$, $50$ and $200$} 
769 
\label{fig:onedheatout} 
\label{fig:onedheatout} 
770 
\end{center} 
\end{center} 
771 
\end{figure} 
\end{figure} 
774 
energy over time. 
energy over time. 
775 
For each time step we generate a plot of the temperature distribution and save 
For each time step we generate a plot of the temperature distribution and save 
776 
each to a file. 
each to a file. 
777 
The following is appended to the end of the \verb while loop and creates one 
The following is appended to the end of the \verbwhile loop and creates one 
778 
figure of the temperature distribution. We start by converting the solution to a 
figure of the temperature distribution. We start by converting the solution to a 
779 
tuple and then plotting this against our \textit{x coordinates} \verb plx we 
tuple and then plotting this against our \textit{x coordinates} \verbplx we 
780 
have generated before. We add a title to the diagram before it is rendered into 
have generated before. We add a title to the diagram before it is rendered into 
781 
a file. 
a file. 
782 
Finally, the figure is saved to a \verb*.png file and cleared for the 
Finally, the figure is saved to a \verb*.png file and cleared for the 
797 
\end{python} 
\end{python} 
798 
Some results are shown in \reffig{fig:onedheatout}. 
Some results are shown in \reffig{fig:onedheatout}. 
799 


800 
\subsection{Make a video} 
\subsection{Making a Video} 
801 
Our saved plots from the previous section can be cast into a video using the 
Our saved plots from the previous section can be cast into a video using the 
802 
following command appended to the end of the script. The \verb mencoder command 
following command appended to the end of the script. The \verb mencoder command 
803 
is Linux only, so other platform users need to use an alternative video encoder. 
is not available on every platform, so some users need to use an alternative 
804 

video encoder. 
805 
\begin{python} 
\begin{python} 
806 
# compile the *.png files to create a *.avi videos that show T change 
# compile the *.png files to create a *.avi video that shows T change 
807 
# with time. This operation uses Linux mencoder. For other operating 
# with time. This operation uses Linux mencoder. For other operating 
808 
# systems it is possible to use your favourite video compiler to 
# systems it is possible to use your favourite video compiler to 
809 
# convert image files to videos. 
# convert image files to videos. 