/[escript]/trunk/doc/user/firststep.tex
ViewVC logotype

Diff of /trunk/doc/user/firststep.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 993 by gross, Fri Feb 23 06:39:38 2007 UTC revision 2923 by jfenwick, Thu Feb 4 04:05:36 2010 UTC
# Line 1  Line 1 
1  % $Id$  
2    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3    %
4    % Copyright (c) 2003-2010 by University of Queensland
5    % Earth Systems Science Computational Center (ESSCC)
6    % http://www.uq.edu.au/esscc
7  %  %
8  %           Copyright © 2006 by ACcESS MNRF  % Primary Business: Queensland, Australia
9  %               \url{http://www.access.edu.au  % Licensed under the Open Software License version 3.0
10  %         Primary Business: Queensland, Australia.  % http://www.opensource.org/licenses/osl-3.0.php
 %   Licensed under the Open Software License version 3.0  
 %      http://www.opensource.org/licenses/osl-3.0.php  
11  %  %
12    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13    
14    
15  \section{The First Steps}  \section{The First Steps}
16  \label{FirstSteps}  \label{FirstSteps}
17    
 \begin{figure}  
 \centerline{\includegraphics[width=\figwidth]{figures/FirstStepDomain}}  
 \caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}  
 \label{fig:FirstSteps.1}  
 \end{figure}  
18    
19  In this chapter we will give an introduction how to use \escript to solve  
20  a partial differential equation \index{partial differential equation} (PDE \index{partial differential equation!PDE}). The reader should be familiar with Python. The knowledge presented at the Python tutorial at \url{http://docs.python.org/tut/tut.html}  In this chapter we give an introduction how to use \escript to solve
21  is sufficient. It is helpful if the reader has some basic knowledge of PDEs \index{partial differential equation}.  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}  The PDE \index{partial differential equation} we wish to solve is the Poisson equation \index{Poisson equation}
25  \begin{equation}  \begin{equation}
26  -\Delta u =f  -\Delta u =f
27  \label{eq:FirstSteps.1}  \label{eq:FirstSteps.1}
28  \end{equation}  \end{equation}
29  for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$  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  is the unit square
31  \begin{equation}  \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 \}  \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}  \label{eq:FirstSteps.1b}
34  \end{equation}  \end{equation}
35  The domain is shown in \fig{fig:FirstSteps.1}.  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  $\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by
43  \begin{equation}  \begin{equation}
44  \Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}  \Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}
45  \label{eq:FirstSteps.1.1}  \label{eq:FirstSteps.1.1}
46  \end{equation}  \end{equation}
47  where, for any function $w$ and any direction $i$, $u\hackscore{,i}$  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$.    denotes the partial derivative \index{partial derivative} of $u$ with respect to $i$.  
49  \footnote{Some readers  \footnote{You
50  may be more familiar with the Laplace operator\index{Laplace operator} being written  may be more familiar with the Laplace operator\index{Laplace operator} being written
51  as $\nabla^2$, and written in the form  as $\nabla^2$, and written in the form
52  \begin{equation*}  \begin{equation*}
# Line 54  and \eqn{eq:FirstSteps.1} as Line 59  and \eqn{eq:FirstSteps.1} as
59  \end{equation*}  \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  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 $w\hackscore{,ij}=(w\hackscore {,i})\hackscore{,j}$  to the index. To get a more compact form we will write $u\hackscore{,ij}=(u\hackscore {,i})\hackscore{,j}$
63  which leads to  which leads to
64  \begin{equation}  \begin{equation}
65  \Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii}  \Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii}
66  \label{eq:FirstSteps.1.1b}  \label{eq:FirstSteps.1.1b}
67  \end{equation}  \end{equation}
68  In some cases, and we will see examples for this in the next chapter,  We often find that use
69  the usage of the nested $\sum$ symbols blows up the formulas and therefore  of nested $\sum$ symbols makes formulas cumbersome, and we use the more
70  it is convenient to use the Einstein summation convention \index{summation convention}. This  convenient Einstein summation convention \index{summation convention}. This
71  drops the $\sum$ sign and assumes that a summation over a repeated index is performed  drops the $\sum$ sign and assumes that a summation is performed over any repeated index.
72  ("repeated index means summation"). For instance we write  For instance we write
73  \begin{eqnarray}  \begin{eqnarray}
74  x\hackscore{i}y\hackscore{i}=\sum\hackscore{i=0}^2 x\hackscore{i}y\hackscore{i}   \\  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}   \\  x\hackscore{i}u\hackscore{,i}=\sum\hackscore{i=0}^2 x\hackscore{i}u\hackscore{,i}   \\
# Line 93  are applying the Einstein summation conv Line 98  are applying the Einstein summation conv
98  $n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} +  $n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} +
99  n\hackscore{1} u\hackscore{,1}$.  n\hackscore{1} u\hackscore{,1}$.
100  \footnote{Some readers may familiar with the notation  \footnote{Some readers may familiar with the notation
101  \begin{equation*}  $
102  \frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i}  \frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i}
103  \end{equation*}  $
104  for the normal derivative.}  for the normal derivative.}
105  The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the  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:  set $\Gamma^N$ which is the top and right edge of the domain:
# Line 120  with the Neumann boundary condition \eqn Line 125  with the Neumann boundary condition \eqn
125  Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so  Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so
126  called boundary value  called boundary value
127  problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for  problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for
128  the unknown  the unknown function~$u$.
 function $u$.  
129    
130    
131  \begin{figure}  \begin{figure}[ht]
132  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh.eps}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh}}
133  \caption{Mesh of $4 \time 4$ elements on a rectangular domain.  Here  \caption{Mesh of $4 \time 4$ elements on a rectangular domain.  Here
134  each element is a quadrilateral and described by four nodes, namely  each element is a quadrilateral and described by four nodes, namely
135  the corner points. The solution is interpolated by a bi-linear  the corner points. The solution is interpolated by a bi-linear
# Line 149  elements is called a mesh\index{finite e Line 153  elements is called a mesh\index{finite e
153  method!mesh}. \fig{fig:FirstSteps.2} shows an  method!mesh}. \fig{fig:FirstSteps.2} shows an
154  example of a FEM mesh with four elements in the $x_0$ and four elements  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.    in the $x_1$ direction over the unit square.  
156  For more details we refer the reader to the literature, for instance  For more details we refer the reader to the literature, for instance \Ref{Zienc,NumHand}.
157  \Ref{Zienc,NumHand}.  
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{run-escript} command\footnote{\program{run-escript} 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      run-escript
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}.  \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}  (We will discuss a more general form of a PDE \index{partial differential equation!PDE}
# Line 162  to solve the PDE. Here we are using the Line 188  to solve the PDE. Here we are using the
188  method}. The following statements create the \Domain object \var{mydomain} from the  method}. The following statements create the \Domain object \var{mydomain} from the
189  \finley method \method{Rectangle}  \finley method \method{Rectangle}
190  \begin{python}  \begin{python}
191  from esys.finley import Rectangle    from esys.finley import Rectangle
192  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)    mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
193  \end{python}  \end{python}
194  In this case the domain is a rectangle with the lower, left corner at point $(0,0)$ and  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)$.  the right, upper corner at $(\var{l0},\var{l1})=(1,1)$.
# Line 175  see \Chap{CHAPTER ON FINLEY}. Line 201  see \Chap{CHAPTER ON FINLEY}.
201  The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and  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$:  the right hand side $f$ of the PDE to constant $1$:
203  \begin{python}  \begin{python}
204  from esys.escript.linearPDEs import Poisson    from esys.escript.linearPDEs import Poisson
205  mypde = Poisson(mydomain)    mypde = Poisson(mydomain)
206  mypde.setValue(f=1)    mypde.setValue(f=1)
207  \end{python}  \end{python}
208  We have not specified any boundary condition but the  We have not specified any boundary condition but the
209  \Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann  \Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann
# Line 189  by defining a characteristic function \i Line 215  by defining a characteristic function \i
215  which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$ where Dirichlet boundary condition is set  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},  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  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 represents locations in the domain one uses  an object \var{x} which contains the coordinates of the nodes in the domain use
219  \begin{python}  \begin{python}
220  x=mydomain.getX()    x=mydomain.getX()
221  \end{python}  \end{python}
222  The method \method{getX} of the \Domain \var{mydomain}  The method \method{getX} of the \Domain \var{mydomain}
223  gives access to locations  gives access to locations
224  in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which is  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 details. What we need to know here is that  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 (=tuple of dimensions) which can be checked by  \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:  calling the \method{getRank} and \method{getShape} methods:
229  \begin{python}  \begin{python}
230  print "rank ",x.getRank(),", shape ",x.getShape()    print "rank ",x.getRank(),", shape ",x.getShape()
231  \end{python}  \end{python}
232  will print something like  This will print something like
233  \begin{python}  \begin{python}
234  rank 1, shape (2,)    rank 1, shape (2,)
235  \end{python}  \end{python}
236  The \Data object also maintains type information which is represented by the  The \Data object also maintains type information which is represented by the
237  \FunctionSpace of the object. For instance  \FunctionSpace of the object. For instance
238  \begin{python}  \begin{python}
239  print x.getFunctionSpace()    print x.getFunctionSpace()
240  \end{python}  \end{python}
241  will print  will print
242  \begin{python}  \begin{python}
243  Function space type: Finley_Nodes on FinleyMesh    Function space type: Finley_Nodes on FinleyMesh
244  \end{python}  \end{python}
245  which tells us that the coordinates are stored on the nodes of a \finley mesh.  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  To get the  $x\hackscore{0}$ coordinates of the locations we use the
247  statement  statement
248  \begin{python}  \begin{python}
249  x0=x[0]    x0=x[0]
250  \end{python}  \end{python}
251  Object \var{x0}  Object \var{x0}
252  is again a \Data object now with \Rank $0$ and  is again a \Data object now with \Rank $0$ and
253  \Shape $()$. It inherits the \FunctionSpace from \var{x}:  \Shape $()$. It inherits the \FunctionSpace from \var{x}:
254  \begin{python}  \begin{python}
255  print x0.getRank(),x0.getShape(),x0.getFunctionSpace()    print x0.getRank(),x0.getShape(),x0.getFunctionSpace()
256  \end{python}  \end{python}
257  will print  will print
258  \begin{python}  \begin{python}
259  0 () Function space type: Finley_Nodes on FinleyMesh    0 () Function space type: Finley_Nodes on FinleyMesh
260  \end{python}  \end{python}
261  We can now construct the function \var{gammaD} by  We can now construct a function \var{gammaD} which is only non-zero on the bottom and left edges
262    of the domain with
263  \begin{python}  \begin{python}
264  from esys.escript import whereZero    from esys.escript import whereZero
265  gammaD=whereZero(x[0])+whereZero(x[1])    gammaD=whereZero(x[0])+whereZero(x[1])
266  \end{python}  \end{python}
267  where  
268  \code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (allmost) equal to zero  \code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (almost) equal to zero
269  and $0$ elsewhere.  and $0$ elsewhere.
270  Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is  Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is
271  equal to zero and $0$ elsewhere.  equal to zero and $0$ elsewhere.
272  The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])}  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 exactly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero.  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  Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from
275  \begin{python}  \begin{python}
276  print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace()    print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace()
277  \end{python}  \end{python}
278  one gets  one gets
279  \begin{python}  \begin{python}
280  0 () Function space type: Finley_Nodes on FinleyMesh    0 () Function space type: Finley_Nodes on FinleyMesh
281  \end{python}  \end{python}
282  The additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the  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  characteristic function \index{characteristic function} of the locations
284  of the domain where homogeneous Dirichlet boundary condition \index{Dirichlet boundary condition!homogeneous}  of the domain where homogeneous Dirichlet boundary condition \index{Dirichlet boundary condition!homogeneous}
285  are set. The complete definition of our example is now:  are set. The complete definition of our example is now:
286  \begin{python}  \begin{python}
287  from esys.linearPDEs import Poisson    from esys.linearPDEs import Poisson
288  x = mydomain.getX()    x = mydomain.getX()
289  gammaD = whereZero(x[0])+whereZero(x[1])    gammaD = whereZero(x[0])+whereZero(x[1])
290  mypde = Poisson(domain=mydomain)    mypde = Poisson(domain=mydomain)
291  mypde.setValue(f=1,q=gammaD)    mypde.setValue(f=1,q=gammaD)
292  \end{python}  \end{python}
293  The first statement imports the \Poisson class definition from the \linearPDEs module \escript package.  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  To get the solution of the Poisson equation defined by \var{mypde} we just have to call its
295  \method{getSolution}.  \method{getSolution}.
296    
297  Now we can write the script to solve our test problem (Remember that  Now we can write the script to solve our Poisson problem
298  lines starting with '\#' are comment lines in Python) (available as \file{poisson.py}  \begin{python}
299  in the \ExampleDirectory):    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}  \begin{python}
375  from esys.escript import *  from esys.escript import *
376  from esys.escript.linearPDEs import Poisson  from esys.escript.linearPDEs import Poisson
377  from esys.finley import Rectangle  from esys.finley import Rectangle
378    import numpy
379    import matplotlib
380    import pylab
381  # generate domain:  # generate domain:
382  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
383  # define characteristic function of Gamma^D  # define characteristic function of Gamma^D
# Line 283  gammaD = whereZero(x[0])+whereZero(x[1]) Line 387  gammaD = whereZero(x[0])+whereZero(x[1])
387  mypde = Poisson(domain=mydomain)  mypde = Poisson(domain=mydomain)
388  mypde.setValue(f=1,q=gammaD)  mypde.setValue(f=1,q=gammaD)
389  u = mypde.getSolution()  u = mypde.getSolution()
390  # write u to an external file  # interpolate u to a matplotlib grid:
391  saveVTK("u.xml",sol=u)  x_grid = numpy.linspace(0.,1.,50)
392  \end{python}  y_grid = numpy.linspace(0.,1.,50)
393  The last statement writes the solution tagged with the name "sol" to the external file \file{u.xml} in  x=mydomain.getX()[0].toListOfTuples()
394  \VTK file format. \VTK is a software library  y=mydomain.getX()[1].toListOfTuples()
395  for the visualization of scientific, engineering and analytical data and is freely available  z=interpolate(u,mydomain.getX().getFunctionSpace())
396  from \url{http://www.vtk.org}. There are a variety of graphical user interfaces  z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )
397  for \VTK available, for instance \mayavi which can be downloaded from \url{http://mayavi.sourceforge.net/} but is also available on most  # interpolate u to a rectangular grid:
398  \LINUX distributions.  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      run-escript 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 two-dimensional 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}  \begin{figure}
414  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult.eps}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}
415  \caption{Visualization of the Poisson Equation Solution for $f=1$}  \caption{Visualization of the Poisson Equation Solution for $f=1$}
416  \label{fig:FirstSteps.3}  \label{fig:FirstSteps.3}
417  \end{figure}  \end{figure}
418    
419  You can edit the script file using your favourite text editor (or the Integrated DeveLopment Environment IDLE  \subsection{Visualization using \VTK}
 for Python, see \url{http://idlefork.sourceforge.net}). If the script file has the name \file{poisson.py} \index{scripts!\file{poisson.py}} you can run the  
 script from any shell using the command:  
 \begin{python}  
 python poisson.py  
 \end{python}  
 After the script has (hopefully successfully) been completed you will find the file \file{u.xml} in the current  
 directory. An easy way to visualize the results is the command  
 \begin{python}  
 mayavi -d u.xml -m SurfaceMap &  
 \end{python}  
 to show the results, see \fig{fig:FirstSteps.3}.  
   
 \subsection{Visualization with \pyvisi}  
   
 \subsection{Domain description using with \pycad}  
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      run-escript poisson\hackscore VTK.py
452      mayavi2 -d u.xml -m SurfaceMap
453    \end{verbatim}
454    The result is shown in Figure~\fig{fig:FirstSteps.3}.

Legend:
Removed from v.993  
changed lines
  Added in v.2923

  ViewVC Help
Powered by ViewVC 1.1.26