1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032018 by The University of Queensland 
4 
% http://www.uq.edu.au 
5 
% 
6 
% Primary Business: Queensland, Australia 
7 
% Licensed under the Apache License, version 2.0 
8 
% http://www.apache.org/licenses/LICENSE2.0 
9 
% 
10 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
% Development 20122013 by School of Earth Sciences 
12 
% Development from 2014 by Centre for Geoscience Computing (GeoComp) 
13 
% 
14 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 

16 
\section{The First Steps}\label{FirstSteps} 
17 
This chapter is an introduction on how to use \escript to solve 
18 
a partial differential equation\index{partial differential equation} (PDE\index{partial differential equation!PDE}). 
19 
We assume you are at least a little familiar with \PYTHON. 
20 
The knowledge presented in the \PYTHON tutorial at \url{https://docs.python.org/2/tutorial/} is more than sufficient. 
21 

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

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

85 
On the boundary of the domain $\Omega$ the normal derivative $n_{i} u_{,i}$ 
86 
of the solution $u$ shall be zero, i.e. $u$ shall fulfill 
87 
the homogeneous Neumann boundary condition\index{Neumann 
88 
boundary condition!homogeneous} 
89 
\begin{equation} 
90 
n_{i} u_{,i}= 0 \;. 
91 
\label{eq:FirstSteps.2} 
92 
\end{equation} 
93 
$n=(n_{i})$ denotes the outer normal field 
94 
of the domain, see \fig{fig:FirstSteps.1}. Remember that we 
95 
apply the Einstein summation convention \index{summation convention}, i.e. $n_{i} u_{,i}= n_{0} u_{,0} +% 
96 
n_{1} u_{,1}$.\footnote{Some readers may be more familiar with the 
97 
notation $\frac{\partial u}{\partial n} = n_{i} u_{,i}=\mathbf{n}\cdot \nabla u$.} 
98 
The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the 
99 
set $\Gamma^N$, 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, 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 
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 Neumann \eqn{eq:FirstSteps.2} and 
119 
Dirichlet boundary conditions 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 are 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 \Refe{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 can be edited with your favorite editor. 
153 
To enter the escript environment, use the \program{runescript} 
154 
command\footnote{\program{runescript} is not available under Windows. 
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.7.6 (default, Mar 22 2014, 15:40:47) 
163 
[GCC 4.8.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\footnote{Throughout our examples, we use the python 3 form of 
168 
print. That is, print(1) instead of print 1.}, for instance 
169 
\begin{python} 
170 
>>> x=2+3 
171 
>>> print("2+3=",x) 
172 
2+3= 5 
173 
\end{python} 
174 
We refer to the \PYTHON user's guide if you are not familiar with \PYTHON. 
175 

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

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

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

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

319 
\subsection{Plotting Using \MATPLOTLIB} 
320 
The \MATPLOTLIB module provides a simple and easytouse way to visualize PDE 
321 
solutions (or other \Data objects). 
322 
To hand over data from \escript to \MATPLOTLIB the values need to be mapped onto 
323 
a rectangular grid. We will make use of the \numpy module for this. 
324 

325 
First we need to create a rectangular grid which is accomplished by 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 grid while 
332 
\var{y_grid} defines the y coordinates of the grid. 
333 
In this case we use $50$ points over the interval $[0,1]$ in both directions. 
334 

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

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

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

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

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

422 
\subsection{Visualization using export files} 
423 

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

429 
To write the solution \var{u} of the Poisson problem in the \VTK file format 
430 
to the file \file{u.vtu} one needs to add: 
431 
\begin{python} 
432 
from esys.weipa import saveVTK 
433 
saveVTK("u.vtu", sol=u) 
434 
\end{python} 
435 
This file can then be opened in a \VTK compatible visualization tool where the 
436 
solution is accessible by the name {\it sol}. Similarly, 
437 
\begin{python} 
438 
from esys.weipa import saveSilo 
439 
saveSilo("u.silo", sol=u) 
440 
\end{python} 
441 
will write \var{u} to a \SILO file if escript was compiled with support for 
442 
LLNL's \SILO library. 
443 

444 
The Poisson problem script is now 
445 
\begin{python} 
446 
from esys.escript import * 
447 
from esys.escript.linearPDEs import Poisson 
448 
from esys.finley import Rectangle 
449 
from esys.weipa import saveVTK 
450 
# generate domain: 
451 
mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20) 
452 
# define characteristic function of Gamma^D 
453 
x = mydomain.getX() 
454 
gammaD = whereZero(x[0])+whereZero(x[1]) 
455 
# define PDE and get its solution u 
456 
mypde = Poisson(domain=mydomain) 
457 
mypde.setValue(f=1,q=gammaD) 
458 
u = mypde.getSolution() 
459 
# write u to an external file 
460 
saveVTK("u.vtu",sol=u) 
461 
\end{python} 
462 
The entire code is available as \file{poisson_vtk.py} in the \ExampleDirectory. 
463 

464 
You can run the script using the {\it escript} environment and visualize the 
465 
solution using \mayavi: 
466 
\begin{verbatim} 
467 
runescript poisson_vtk.py 
468 
mayavi2 d u.vtu m Surface 
469 
\end{verbatim} 
470 
The result is shown in \fig{fig:FirstSteps.3}. 
471 
