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 

15 
\section{The First Steps}\label{FirstSteps} 
16 
In this chapter we give an introduction how to use \escript to solve 
17 
a partial differential equation\index{partial differential equation} (PDE\index{partial differential equation!PDE}). 
18 
We assume you are at least a little familiar with Python. 
19 
The knowledge presented in the Python tutorial at \url{http://docs.python.org/tut/tut.html} is more than sufficient. 
20 

21 
The PDE\index{partial differential equation} we wish to solve is the Poisson equation\index{Poisson equation} 
22 
\begin{equation} 
23 
\Delta u=f 
24 
\label{eq:FirstSteps.1} 
25 
\end{equation} 
26 
for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$, 
27 
is the unit square 
28 
\begin{equation} 
29 
\Omega=[0,1]^2=\{ (x_0;x_1)  0\le x_{0} \le 1 \mbox{ and } 0\le x_{1} \le 1 \} 
30 
\label{eq:FirstSteps.1b} 
31 
\end{equation} 
32 
The domain is shown in \fig{fig:FirstSteps.1}. 
33 
\begin{figure}[ht] 
34 
\centerline{\includegraphics{FirstStepDomain}} 
35 
\caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.} 
36 
\label{fig:FirstSteps.1} 
37 
\end{figure} 
38 

39 
$\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by 
40 
\begin{equation} 
41 
\Delta u = (u_{,0})_{,0}+(u_{,1})_{,1} 
42 
\label{eq:FirstSteps.1.1} 
43 
\end{equation} 
44 
where, for any function $u$ and any direction $i$, $u_{,i}$ 
45 
denotes the partial derivative \index{partial derivative} of $u$ with respect 
46 
to $i$.\footnote{You may be more familiar with the Laplace 
47 
operator\index{Laplace operator} being written as $\nabla^2$, and written in 
48 
the form 
49 
\begin{equation*} 
50 
\nabla^2 u = \nabla^t \cdot \nabla u = \frac{\partial^2 u}{\partial x_0^2} 
51 
+ \frac{\partial^2 u}{\partial x_1^2} 
52 
\end{equation*} 
53 
and \eqn{eq:FirstSteps.1} as 
54 
\begin{equation*} 
55 
\nabla^2 u = f 
56 
\end{equation*} 
57 
} 
58 
Basically, in the subindex of a function, any index to the left of the comma denotes a spatial derivative with respect 
59 
to the index. To get a more compact form we will write $u_{,ij}=(u_{,i})_{,j}$ 
60 
which leads to 
61 
\begin{equation} 
62 
\Delta u = u_{,00}+u_{,11}=\sum_{i=0}^2 u_{,ii} 
63 
\label{eq:FirstSteps.1.1b} 
64 
\end{equation} 
65 
We often find that use 
66 
of nested $\sum$ symbols makes formulas cumbersome, and we use the more 
67 
compact Einstein summation convention\index{summation convention}. This 
68 
drops the $\sum$ sign and assumes that a summation is performed over any repeated index. 
69 
For instance we write 
70 
\begin{eqnarray} 
71 
x_{i}y_{i}=\sum_{i=0}^2 x_{i}y_{i} \\ 
72 
x_{i}u_{,i}=\sum_{i=0}^2 x_{i}u_{,i} \\ 
73 
u_{,ii}=\sum_{i=0}^2 u_{,ii} \\ 
74 
x_{ij}u_{i,j}=\sum_{j=0}^2\sum_{i=0}^2 x_{ij}u_{i,j} \\ 
75 
\label{eq:FirstSteps.1.1c} 
76 
\end{eqnarray} 
77 
With the summation convention we can write the Poisson equation \index{Poisson equation} as 
78 
\begin{equation} 
79 
 u_{,ii} =1 
80 
\label{eq:FirstSteps.1.sum} 
81 
\end{equation} 
82 
where $f=1$ in this example. 
83 

84 
On the boundary of the domain $\Omega$ the normal derivative $n_{i} u_{,i}$ 
85 
of the solution $u$ shall be zero, i.e. $u$ shall fulfill 
86 
the homogeneous Neumann boundary condition\index{Neumann 
87 
boundary condition!homogeneous} 
88 
\begin{equation} 
89 
n_{i} u_{,i}= 0 \;. 
90 
\label{eq:FirstSteps.2} 
91 
\end{equation} 
92 
$n=(n_{i})$ denotes the outer normal field 
93 
of the domain, see \fig{fig:FirstSteps.1}. Remember that we 
94 
are applying the Einstein summation convention \index{summation convention}, i.e. $n_{i} u_{,i}= n_{0} u_{,0} +% 
95 
n_{1} u_{,1}$.\footnote{Some readers may familiar with the 
96 
notation $\frac{\partial u}{\partial n} = n_{i} u_{,i}$ 
97 
for the normal derivative.} 
98 
The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the 
99 
set $\Gamma^N$ which is the top and right edge of the domain: 
100 
\begin{equation} 
101 
\Gamma^N=\{(x_0;x_1) \in \Omega  x_{0}=1 \mbox{ or } x_{1}=1 \} 
102 
\label{eq:FirstSteps.2b} 
103 
\end{equation} 
104 
On the bottom and the left edge of the domain which is defined 
105 
as 
106 
\begin{equation} 
107 
\Gamma^D=\{(x_0;x_1) \in \Omega  x_{0}=0 \mbox{ or } x_{1}=0 \} 
108 
\label{eq:FirstSteps.2c} 
109 
\end{equation} 
110 
the solution shall be identical to zero: 
111 
\begin{equation} 
112 
u=0 \; . 
113 
\label{eq:FirstSteps.2d} 
114 
\end{equation} 
115 
This kind of boundary condition is called a homogeneous Dirichlet boundary 
116 
condition\index{Dirichlet boundary condition!homogeneous}. 
117 
The partial differential equation in \eqn{eq:FirstSteps.1.sum} together 
118 
with the Neumann boundary condition \eqn{eq:FirstSteps.2} and 
119 
Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a socalled 
120 
boundary value 
121 
problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) 
122 
for the unknown function~$u$. 
123 

124 
\begin{figure}[ht] 
125 
\centerline{\includegraphics{FirstStepMesh}} 
126 
\caption{Mesh of $4 \times 4$ elements on a rectangular domain. Here 
127 
each element is a quadrilateral and described by four nodes, namely 
128 
the corner points. The solution is interpolated by a bilinear 
129 
polynomial.} 
130 
\label{fig:FirstSteps.2} 
131 
\end{figure} 
132 

133 
In general the BVP\index{boundary value problem!BVP} cannot be solved 
134 
analytically and numerical methods have to be used to construct an 
135 
approximation of the solution $u$. 
136 
Here we will use the finite element method\index{finite element method} 
137 
(FEM\index{finite element method!FEM}). 
138 
The basic idea is to fill the domain with a set of points called nodes. 
139 
The solution is approximated by its values on the nodes\index{finite element method!nodes}. 
140 
Moreover, the domain is subdivided into smaller subdomains called 
141 
elements\index{finite element method!element}. 
142 
On each element the solution is represented by a polynomial of a certain 
143 
degree through its values at the nodes located in the element. 
144 
The nodes and their connection through elements is called a 
145 
mesh\index{finite element method!mesh}. \fig{fig:FirstSteps.2} shows an 
146 
example of a FEM mesh with four elements in the $x_0$ and four elements 
147 
in the $x_1$ direction over the unit square. 
148 
For more details we refer the reader to the literature, for instance \Ref{Zienc,NumHand}. 
149 

150 
The \escript solver we want to use to solve this problem is embedded into the python interpreter language. 
151 
So you can solve the problem interactively but you will learn quickly that it 
152 
is more efficient to use scripts which you can edit with your favorite editor. 
153 
To enter the escript environment, use the \program{runescript} 
154 
command\footnote{\program{runescript} is not available under Windows yet. 
155 
If you run under Windows you can just use the \program{python} command and the 
156 
\env{OMP_NUM_THREADS} environment variable to control the number of threads.}: 
157 
\begin{verbatim} 
158 
runescript 
159 
\end{verbatim} 
160 
which will pass you on to the python prompt 
161 
\begin{verbatim} 
162 
Python 2.5.2 (r252:60911, Oct 5 2008, 19:29:17) 
163 
[GCC 4.3.2] on linux2 
164 
Type "help", "copyright", "credits" or "license" for more information. 
165 
>>> 
166 
\end{verbatim} 
167 
Here you can use all available python commands and language features, for instance 
168 
\begin{python} 
169 
>>> x=2+3 
170 
>>> print "2+3=",x 
171 
2+3= 5 
172 
\end{python} 
173 
We refer to the python user's guide if you not familiar with python. 
174 

175 
\escript provides the class \Poisson to define a Poisson equation\index{Poisson equation}. 
176 
(We will discuss a more general form of a PDE\index{partial differential equation!PDE} 
177 
that can be defined through the \LinearPDE class later.) 
178 
The instantiation of a \Poisson class object requires the specification of the domain $\Omega$. 
179 
In \escript the \Domain class objects are used to describe the geometry of a 
180 
domain but it also contains information about the discretization methods and 
181 
the actual solver which is used to solve the PDE. 
182 
Here we are using the FEM library \finley\index{finite element method}. 
183 
The following statements create the \Domain object \var{mydomain} from the 
184 
\finley method \method{Rectangle}: 
185 
\begin{python} 
186 
from esys.finley import Rectangle 
187 
mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20) 
188 
\end{python} 
189 
In this case the domain is a rectangle with the lower, left corner at point $(0,0)$ and 
190 
the right, upper corner at $(\var{l0},\var{l1})=(1,1)$. 
191 
The arguments \var{n0} and \var{n1} define the number of elements in $x_{0}$ and 
192 
$x_{1}$direction respectively. For more details on \method{Rectangle} and 
193 
other \Domain generators within the \finley module, 
194 
see \Chap{CHAPTER ON FINLEY}. 
195 

196 
The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and 
197 
the right hand side $f$ of the PDE to constant $1$: 
198 
\begin{python} 
199 
from esys.escript.linearPDEs import Poisson 
200 
mypde = Poisson(mydomain) 
201 
mypde.setValue(f=1) 
202 
\end{python} 
203 
We have not specified any boundary condition but the \Poisson class implicitly 
204 
assumes homogeneous Neuman boundary conditions\index{Neumann boundary condition!homogeneous} defined by \eqn{eq:FirstSteps.2}. 
205 
With this boundary condition the BVP\index{boundary value problem!BVP} we have 
206 
defined has no unique solution. 
207 
In fact, with any solution $u$ and any constant $C$ the function $u+C$ becomes 
208 
a solution as well. 
209 
We have to add a Dirichlet boundary condition\index{Dirichlet boundary condition}. 
210 
This is done by defining a characteristic function\index{characteristic function} 
211 
which has positive values at locations $x=(x_{0},x_{1})$ 
212 
where Dirichlet boundary condition is set and $0$ elsewhere. 
213 
In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c}, we need to 
214 
construct a function \var{gammaD} which is positive for the cases $x_{0}=0$ or $x_{1}=0$. 
215 
To get an object \var{x} which contains the coordinates of the nodes in the domain use 
216 
\begin{python} 
217 
x=mydomain.getX() 
218 
\end{python} 
219 
The method \method{getX} of the \Domain \var{mydomain} gives access to locations 
220 
in the domain defined by \var{mydomain}. 
221 
The object \var{x} is actually a \Data object which will be discussed in 
222 
\Chap{ESCRIPT CHAP} in more detail. 
223 
What we need to know here is that \var{x} has \Rank (number of dimensions) and 
224 
a \Shape (list of dimensions) which can be viewed by calling the \method{getRank} and \method{getShape} methods: 
225 
\begin{python} 
226 
print "rank ",x.getRank(),", shape ",x.getShape() 
227 
\end{python} 
228 
This will print something like 
229 
\begin{python} 
230 
rank 1, shape (2,) 
231 
\end{python} 
232 
The \Data object also maintains type information which is represented by the 
233 
\FunctionSpace of the object. For instance 
234 
\begin{python} 
235 
print x.getFunctionSpace() 
236 
\end{python} 
237 
will print 
238 
\begin{python} 
239 
Function space type: Finley_Nodes on FinleyMesh 
240 
\end{python} 
241 
which tells us that the coordinates are stored on the nodes of (rather than on 
242 
points in the interior of) a Finley mesh. 
243 
To get the $x_{0}$ coordinates of the locations we use the statement 
244 
\begin{python} 
245 
x0=x[0] 
246 
\end{python} 
247 
Object \var{x0} is again a \Data object now with \Rank $0$ and \Shape $()$. 
248 
It inherits the \FunctionSpace from \var{x}: 
249 
\begin{python} 
250 
print x0.getRank(), x0.getShape(), x0.getFunctionSpace() 
251 
\end{python} 
252 
will print 
253 
\begin{python} 
254 
0 () Function space type: Finley_Nodes on FinleyMesh 
255 
\end{python} 
256 
We can now construct a function \var{gammaD} which is only nonzero on the 
257 
bottom and left edges of the domain with 
258 
\begin{python} 
259 
from esys.escript import whereZero 
260 
gammaD=whereZero(x[0])+whereZero(x[1]) 
261 
\end{python} 
262 

263 
\code{whereZero(x[0])} creates a function which equals $1$ where \code{x[0]} is (almost) equal to zero and $0$ elsewhere. 
264 
Similarly, \code{whereZero(x[1])} creates a function which equals $1$ where \code{x[1]} is equal to zero and $0$ elsewhere. 
265 
The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])} 
266 
gives a function on the domain \var{mydomain} which is strictly positive where $x_{0}$ or $x_{1}$ is equal to zero. 
267 
Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. 
268 
So from 
269 
\begin{python} 
270 
print gammaD.getRank(), gammaD.getShape(), gammaD.getFunctionSpace() 
271 
\end{python} 
272 
one gets 
273 
\begin{python} 
274 
0 () Function space type: Finley_Nodes on FinleyMesh 
275 
\end{python} 
276 
An additional parameter \var{q} of the \code{setValue} method of the \Poisson 
277 
class defines the characteristic function\index{characteristic function} of 
278 
the locations of the domain where the homogeneous Dirichlet boundary condition\index{Dirichlet boundary condition!homogeneous} is set. 
279 
The complete definition of our example is now: 
280 
\begin{python} 
281 
from esys.escript.linearPDEs import Poisson 
282 
x = mydomain.getX() 
283 
gammaD = whereZero(x[0])+whereZero(x[1]) 
284 
mypde = Poisson(domain=mydomain) 
285 
mypde.setValue(f=1,q=gammaD) 
286 
\end{python} 
287 
The first statement imports the \Poisson class definition from the \linearPDEs module. 
288 
To get the solution of the Poisson equation defined by \var{mypde} we just have to call its \method{getSolution} method. 
289 

290 
Now we can write the script to solve our Poisson problem 
291 
\begin{python} 
292 
from esys.escript import * 
293 
from esys.escript.linearPDEs import Poisson 
294 
from esys.finley import Rectangle 
295 
# generate domain: 
296 
mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20) 
297 
# define characteristic function of Gamma^D 
298 
x = mydomain.getX() 
299 
gammaD = whereZero(x[0])+whereZero(x[1]) 
300 
# define PDE and get its solution u 
301 
mypde = Poisson(domain=mydomain) 
302 
mypde.setValue(f=1,q=gammaD) 
303 
u = mypde.getSolution() 
304 
\end{python} 
305 
The question is what we do with the calculated solution \var{u}. 
306 
Besides postprocessing, e.g. calculating the gradient or the average value, which will be discussed later, plotting the solution is one of the things you might want to do. 
307 
\escript offers two ways to do this, both based on external modules or packages and so data need to converted to hand over the solution. 
308 
The first option is using the \MATPLOTLIB module which allows plotting 2D results relatively quickly from within the Python script, see~\cite{matplotlib}. 
309 
However, there are limitations when using this tool, especially for large problems and when solving 3dimensional problems. 
310 
Therefore, \escript provides functionality to export data as files which can subsequently be read by thirdparty software packages such as \mayavi\cite{mayavi} or \VisIt~\cite{VisIt}. 
311 

312 
\subsection{Plotting Using \MATPLOTLIB} 
313 
The \MATPLOTLIB module provides a simple and easytouse way to visualize PDE solutions (or other \Data objects). 
314 
To hand over data from \escript to \MATPLOTLIB the values need to mapped onto 
315 
a rectangular grid\footnote{Users of Debian 5 (Lenny) please note: this example 
316 
makes use of the \function{griddata} method in \module{matplotlib.mlab}. 
317 
This method is not part of version 0.98.1 which is available with Lenny. 
318 
If you wish to use contour plots, you may need to install a later version. 
319 
Users of Ubuntu 8.10 or later should be fine.}. We will make use of the \numpy module. 
320 

321 
First we need to create a rectangular grid which is accomplished by the following statements: 
322 
\begin{python} 
323 
import numpy 
324 
x_grid = numpy.linspace(0., 1., 50) 
325 
y_grid = numpy.linspace(0., 1., 50) 
326 
\end{python} 
327 
\var{x_grid} is an array defining the x coordinates of the grid while 
328 
\var{y_grid} defines the y coordinates of the grid. 
329 
In this case we use $50$ points over the interval $[0,1]$ in both directions. 
330 

331 
Now the values created by \escript need to be interpolated to this grid. 
332 
We will use the \MATPLOTLIB \function{mlab.griddata} function to do this. 
333 
Spatial coordinates are easily extracted as a \var{list} by 
334 
\begin{python} 
335 
x=mydomain.getX()[0].toListOfTuples() 
336 
y=mydomain.getX()[1].toListOfTuples() 
337 
\end{python} 
338 
In principle we can apply the same \member{toListOfTuples} method to extract the values from the PDE solution \var{u}. 
339 
However, we have to make sure that the \Data object we extract the values from 
340 
uses the same \FunctionSpace as we have used when extracting \var{x} and \var{y}. 
341 
We apply the \function{interpolation} to \var{u} before extraction to achieve this: 
342 
\begin{python} 
343 
z=interpolate(u, mydomain.getX().getFunctionSpace()) 
344 
\end{python} 
345 
The values in \var{z} are the values at the points with the coordinates given by \var{x} and \var{y}. 
346 
These values are interpolated to the grid defined by \var{x_grid} and \var{y_grid} by using 
347 
\begin{python} 
348 
import matplotlib 
349 
z_grid = matplotlib.mlab.griddata(x, y, z, xi=x_grid, yi=y_grid) 
350 
\end{python} 
351 
Now \var{z_grid} gives the values of the PDE solution \var{u} at the grid which can be plotted using \function{contourf}: 
352 
\begin{python} 
353 
matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5) 
354 
matplotlib.pyplot.savefig("u.png") 
355 
\end{python} 
356 
Here we use 5 contours. The last statement writes the plot to the file \file{u.png} in the PNG format. 
357 
Alternatively, one can use 
358 
\begin{python} 
359 
matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5) 
360 
matplotlib.pyplot.show() 
361 
\end{python} 
362 
which gives an interactive browser window. 
363 

364 
\begin{figure} 
365 
\centerline{\includegraphics[width=\figwidth]{FirstStepResultMATPLOTLIB}} 
366 
\caption{Visualization of the Poisson Equation Solution for $f=1$ using \MATPLOTLIB} 
367 
\label{fig:FirstSteps.3b} 
368 
\end{figure} 
369 

370 
Now we can write the script to solve our Poisson problem 
371 
\begin{python} 
372 
from esys.escript import * 
373 
from esys.escript.linearPDEs import Poisson 
374 
from esys.finley import Rectangle 
375 
import numpy 
376 
import matplotlib 
377 
import pylab 
378 
# generate domain: 
379 
mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20) 
380 
# define characteristic function of Gamma^D 
381 
x = mydomain.getX() 
382 
gammaD = whereZero(x[0])+whereZero(x[1]) 
383 
# define PDE and get its solution u 
384 
mypde = Poisson(domain=mydomain) 
385 
mypde.setValue(f=1,q=gammaD) 
386 
u = mypde.getSolution() 
387 
# interpolate u to a matplotlib grid: 
388 
x_grid = numpy.linspace(0.,1.,50) 
389 
y_grid = numpy.linspace(0.,1.,50) 
390 
x=mydomain.getX()[0].toListOfTuples() 
391 
y=mydomain.getX()[1].toListOfTuples() 
392 
z=interpolate(u,mydomain.getX().getFunctionSpace()) 
393 
z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid ) 
394 
# interpolate u to a rectangular grid: 
395 
matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5) 
396 
matplotlib.pyplot.savefig("u.png") 
397 
\end{python} 
398 
The entire code is available as \file{poisson_matplotlib.py} in the \ExampleDirectory. 
399 
You can run the script using the {\it escript} environment 
400 
\begin{verbatim} 
401 
runescript poisson_matplotlib.py 
402 
\end{verbatim} 
403 
This will create the \file{u.png}, see \fig{fig:FirstSteps.3b}. 
404 
For details on the usage of the \MATPLOTLIB module we refer to the documentation~\cite{matplotlib}. 
405 

406 
As pointed out, \MATPLOTLIB is restricted to the twodimensional case and 
407 
should be used for small problems only. 
408 
It can not be used under \MPI as the \member{toListOfTuples} method is not 
409 
safe under \MPI\footnote{The phrase 'safe under \MPI' means that a program 
410 
will produce correct results when run on more than one processor under \MPI.}. 
411 

412 
\begin{figure} 
413 
\centerline{\includegraphics[width=\figwidth]{FirstStepResult}} 
414 
\caption{Visualization of the Poisson Equation Solution for $f=1$} 
415 
\label{fig:FirstSteps.3} 
416 
\end{figure} 
417 

418 
\subsection{Visualization using export files} 
419 

420 
As an alternative to \MATPLOTLIB, {\it escript} supports exporting data to 
421 
\VTK and \SILO files which can be read by visualization tools such as 
422 
mayavi\cite{mayavi} and \VisIt~\cite{VisIt}. This method is \MPI safe and 
423 
works with large 2D and 3D problems. 
424 

425 
To write the solution \var{u} of the Poisson problem in the \VTK file format 
426 
to the file \file{u.vtu} one needs to add: 
427 
\begin{python} 
428 
saveVTK("u.vtu", sol=u) 
429 
\end{python} 
430 
This file can then be opened in a \VTK compatible visualization tool where the 
431 
solution is accessible by the name {\it sol}. 
432 

433 
The Poisson problem script is now 
434 
\begin{python} 
435 
from esys.escript import * 
436 
from esys.escript.linearPDEs import Poisson 
437 
from esys.finley import Rectangle 
438 
# generate domain: 
439 
mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20) 
440 
# define characteristic function of Gamma^D 
441 
x = mydomain.getX() 
442 
gammaD = whereZero(x[0])+whereZero(x[1]) 
443 
# define PDE and get its solution u 
444 
mypde = Poisson(domain=mydomain) 
445 
mypde.setValue(f=1,q=gammaD) 
446 
u = mypde.getSolution() 
447 
# write u to an external file 
448 
saveVTK("u.vtu",sol=u) 
449 
\end{python} 
450 
The entire code is available as \file{poisson_vtk.py} in the \ExampleDirectory. 
451 

452 
You can run the script using the {\it escript} environment and visualize the 
453 
solution using \mayavi: 
454 
\begin{verbatim} 
455 
runescript poisson_vtk.py 
456 
mayavi2 d u.vtu m SurfaceMap 
457 
\end{verbatim} 
458 
The result is shown in \fig{fig:FirstSteps.3}. 
459 
