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} 
16 
\label{FirstSteps} 
17 

18 

19 

20 
In this chapter we give an introduction how to use \escript to solve 
21 
a partial differential equation \index{partial differential equation} (PDE \index{partial differential equation!PDE}). We assume you are at least a little familiar with Python. The knowledge presented at the Python tutorial at \url{http://docs.python.org/tut/tut.html} 
22 
is more than sufficient. 
23 

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

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

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

130 

131 
\begin{figure}[ht] 
132 
\centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh}} 
133 
\caption{Mesh of $4 \time 4$ elements on a rectangular domain. Here 
134 
each element is a quadrilateral and described by four nodes, namely 
135 
the corner points. The solution is interpolated by a bilinear 
136 
polynomial.} 
137 
\label{fig:FirstSteps.2} 
138 
\end{figure} 
139 

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

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

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

201 
The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and 
202 
the right hand side $f$ of the PDE to constant $1$: 
203 
\begin{python} 
204 
from esys.escript.linearPDEs import Poisson 
205 
mypde = Poisson(mydomain) 
206 
mypde.setValue(f=1) 
207 
\end{python} 
208 
We have not specified any boundary condition but the 
209 
\Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann 
210 
boundary condition!homogeneous} defined by \eqn{eq:FirstSteps.2}. With this boundary 
211 
condition the BVP\index{boundary value problem!BVP} we have defined has no unique solution. In fact, with any solution $u$ 
212 
and any constant $C$ the function $u+C$ becomes a solution as well. We have to add 
213 
a Dirichlet boundary condition \index{Dirichlet boundary condition}. This is done 
214 
by defining a characteristic function \index{characteristic function} 
215 
which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$ where Dirichlet boundary condition is set 
216 
and $0$ elsewhere. In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c}, 
217 
we need to construct a function \var{gammaD} which is positive for the cases $x\hackscore{0}=0$ or $x\hackscore{1}=0$. To get 
218 
an object \var{x} which contains the coordinates of the nodes in the domain use 
219 
\begin{python} 
220 
x=mydomain.getX() 
221 
\end{python} 
222 
The method \method{getX} of the \Domain \var{mydomain} 
223 
gives access to locations 
224 
in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which will be 
225 
discussed in \Chap{ESCRIPT CHAP} in more detail. What we need to know here is that 
226 

227 
\var{x} has \Rank (number of dimensions) and a \Shape (list of dimensions) which can be viewed by 
228 
calling the \method{getRank} and \method{getShape} methods: 
229 
\begin{python} 
230 
print "rank ",x.getRank(),", shape ",x.getShape() 
231 
\end{python} 
232 
This will print something like 
233 
\begin{python} 
234 
rank 1, shape (2,) 
235 
\end{python} 
236 
The \Data object also maintains type information which is represented by the 
237 
\FunctionSpace of the object. For instance 
238 
\begin{python} 
239 
print x.getFunctionSpace() 
240 
\end{python} 
241 
will print 
242 
\begin{python} 
243 
Function space type: Finley_Nodes on FinleyMesh 
244 
\end{python} 
245 
which tells us that the coordinates are stored on the nodes of (rather than on points in the interior of) a \finley mesh. 
246 
To get the $x\hackscore{0}$ coordinates of the locations we use the 
247 
statement 
248 
\begin{python} 
249 
x0=x[0] 
250 
\end{python} 
251 
Object \var{x0} 
252 
is again a \Data object now with \Rank $0$ and 
253 
\Shape $()$. It inherits the \FunctionSpace from \var{x}: 
254 
\begin{python} 
255 
print x0.getRank(),x0.getShape(),x0.getFunctionSpace() 
256 
\end{python} 
257 
will print 
258 
\begin{python} 
259 
0 () Function space type: Finley_Nodes on FinleyMesh 
260 
\end{python} 
261 
We can now construct a function \var{gammaD} which is only nonzero on the bottom and left edges 
262 
of the domain with 
263 
\begin{python} 
264 
from esys.escript import whereZero 
265 
gammaD=whereZero(x[0])+whereZero(x[1]) 
266 
\end{python} 
267 

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

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

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

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

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

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

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

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

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

419 
\subsection{Visualization using \VTK} 
420 

421 
As an alternative {\it escript} supports the usage of visualization tools which base on \VTK, eg. mayavi \cite{mayavi}, \VisIt~\cite{VisIt}. In this case the solution is written to a file in the \VTK format. This file the can read by the tool of choice. Using \VTK file is \MPI safe. 
422 

423 
To write the solution \var{u} in Poisson problem to the file \file{u.xml} one need to add the line 
424 
\begin{python} 
425 
saveVTK("u.xml",sol=u) 
426 
\end{python} 
427 
The solution \var{u} is now available in the \file{u.xml} tagged with the name "sol". 
428 

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

448 
You can run the script using the {\it escript} environment 
449 
and visualize the solution using \mayavi: 
450 
\begin{verbatim} 
451 
runescript poisson_VTK.py 
452 
mayavi2 d u.xml m SurfaceMap 
453 
\end{verbatim} 
454 
The result is shown in Figure~\fig{fig:FirstSteps.3}. 