doc/user/firststep.tex

% $Id$

\section{The First Steps}
\label{FirstSteps} 

\begin{figure}
\centerline{\includegraphics[width=\figwidth]{FirstStepDomain}}
\caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}
\label{fig:FirstSteps.1}
\end{figure}

In this chapter we will give an introduction how to use \escript to solve 
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}
is sufficient. It is helpful if the reader has some basic knowledge of PDEs \index{partial differential equation}.

The PDE \index{partial differential equation} we wish to solve is the Poisson equation \index{Poisson equation} 
\begin{equation}
-\Delta u =f 
\label{eq:FirstSteps.1}
\end{equation}
for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$
is the unit square 
\begin{equation}
\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 \}
\label{eq:FirstSteps.1b}
\end{equation}
The domain is shown in \fig{fig:FirstSteps.1}.

$\Delta$ denotes the Laplace operator\index{Laplace operator} which is defined by
\begin{equation}
\Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}
\label{eq:FirstSteps.1.1}
\end{equation}
where, for any function $w$ and any direction $i$, $u\hackscore{,i}$
denotes the partial derivative \index{partial derivative} of $u$ with respect to $i$.  
\footnote{Some readers
may be more familiar with the Laplace operator\index{Laplace operator} being written
as $\nabla^2$, and written in the form
\begin{equation*}
\nabla^2 u = \nabla^t \cdot \nabla u =  \frac{\partial^2 u}{\partial x\hackscore 0^2} 
+ \frac{\partial^2 u}{\partial  x\hackscore 1^2}
\end{equation*}
and \eqn{eq:FirstSteps.1} as
\begin{equation*}
-\nabla^2 u = f
\end{equation*}
}
Basically, in the subindex of a function, any index to the left of the comma denotes a spatial derivative with respect 
to the index. To get a more compact form we will write $w\hackscore{,ij}=(w\hackscore {,i})\hackscore{,j}$
which leads to
\begin{equation}
\Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii}
\label{eq:FirstSteps.1.1b}
\end{equation}
In some cases, and we will see examples for this in the next chapter,
the usage of the nested $\sum$ symbols blows up the formulas and therefore
it is convenient to use the Einstein summation convention \index{summation convention}. This 
drops the $\sum$ sign and assumes that a summation over a repeated index is performed 
("repeated index means summation"). For instance we write
\begin{eqnarray}
x\hackscore{i}y\hackscore{i}=\sum\hackscore{i=0}^2 x\hackscore{i}y\hackscore{i}   \\
x\hackscore{i}u\hackscore{,i}=\sum\hackscore{i=0}^2 x\hackscore{i}u\hackscore{,i}   \\
u\hackscore{,ii}=\sum\hackscore{i=0}^2 u\hackscore{,ii} \\
x\hackscore{ij}u\hackscore{i,j}=\sum\hackscore{j=0}^2\sum\hackscore{i=0}^2 x\hackscore{ij}u\hackscore{i,j}   \\
\label{eq:FirstSteps.1.1c}
\end{eqnarray}
With the summation convention we can write the Poisson equation \index{Poisson equation} as
\begin{equation}
- u\hackscore{,ii} =1 
\label{eq:FirstSteps.1.sum}
\end{equation}
On the boundary of the domain $\Omega$ the normal derivative $n\hackscore{i} u\hackscore{,i}$
of the solution $u$ shall be zero, ie. $u$ shall fulfill
the homogeneous Neumann boundary condition\index{Neumann
boundary condition!homogeneous}
\begin{equation}
n\hackscore{i} u\hackscore{,i}= 0 \;.
\label{eq:FirstSteps.2}
\end{equation}
$n=(n\hackscore{i})$ denotes the outer normal field
of the domain, see \fig{fig:FirstSteps.1}. Remember that we 
are applying the Einstein summation convention \index{summation convention}, i.e
$n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} +
n\hackscore{1} u\hackscore{,1}$. 
\footnote{Some readers may familiar with the notation
\begin{equation*}
\frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i}
\end{equation*}
for the normal derivative.}
The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the
set $\Gamma^N$ which is the top and right edge of the domain:
\begin{equation}
\Gamma^N=\{(x\hackscore 0;x\hackscore 1) \in \Omega | x\hackscore{0}=1 \mbox{ or } x\hackscore{1}=1  \}
\label{eq:FirstSteps.2b}
\end{equation}
On the bottom and the left edge of the domain which is defined
as 
\begin{equation}
\Gamma^D=\{(x\hackscore 0;x\hackscore 1) \in \Omega | x\hackscore{0}=0 \mbox{ or } x\hackscore{1}=0  \}
\label{eq:FirstSteps.2c}
\end{equation}
the solution shall be identically zero:
\begin{equation}
u=0 \; .
\label{eq:FirstSteps.2d}
\end{equation}
This kind of boundary condition is called a homogeneous Dirichlet boundary condition
\index{Dirichlet boundary condition!homogeneous}. The partial differential equation in \eqn{eq:FirstSteps.1.sum} together
with the Neumann boundary condition \eqn{eq:FirstSteps.2} and 
Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so
called boundary value
problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for 
the unknown
function $u$. 


\begin{figure}
\centerline{\includegraphics[width=\figwidth]{FirstStepMesh}}
\caption{Mesh of $4 \time 4$ elements on a rectangular domain.  Here
each element is a quadrilateral and described by four nodes, namely
the corner points. The solution is interpolated by a bi-linear
polynomial.}
\label{fig:FirstSteps.2}
\end{figure}

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

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

The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and
the right hand side $f$ of the PDE to constant $1$: 
\begin{python}
from esys.escript.linearPDEs import Poisson
mypde = Poisson(mydomain)
mypde.setValue(f=1)
\end{python}
We have not specified any boundary condition but the 
\Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann
boundary condition!homogeneous} defined by \eqn{eq:FirstSteps.2}. With this boundary 
condition the BVP\index{boundary value problem!BVP} we have defined has no unique solution. In fact, with any solution $u$
and any constant $C$ the function $u+C$ becomes a solution as well. We have to add 
a Dirichlet boundary condition \index{Dirichlet boundary condition}. This is done 
by defining a characteristic function \index{characteristic function} 
which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$ where Dirichlet boundary condition is set
and $0$ elsewhere. In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c},
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 
an object \var{x} which represents locations in the domain one uses
\begin{python}
x=mydomain.getX() \;.
\end{python}
In fact \var{x} is a \Data object which we will learn more about in Chapter~\ref{X}. At this stage we only have to know
that \var{x} has a 

In the first statement, the method \method{getX} of the \Domain \var{mydomain} 
gives access to locations 
in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which is
discussed in Chpater\ref{X} in more details. What we need to know here is that 
\var{x} has \Rank (=number of dimensions) and a \Shape (=tuple of dimensions) which can be checked by 
calling the \method{getRank} and \method{getShape} methods:
\begin{python}
print "rank ",x.getRank(),", shape ",x.getShape()
\end{python}
will print something like
\begin{python}
rank 1, shape (2,)
\end{python}
The \Data object also maintains type information which is represented by the 
\FunctionSpace of the object. For instance
\begin{python}
print x.getFunctionSpace()
\end{python}
will print 
\begin{python}
Function space type: Finley_Nodes on FinleyMesh 
\end{python}
which tells us that the coordinates are stored on the nodes of a \finley mesh. 
To get the  $x\hackscore{0}$ coordinates of the locations we use the
statement 
\begin{python}
x0=x[0]
\end{python}
Object \var{x0} 
is again a \Data object now with \Rank $0$ and 
\Shape $()$. It inherits the \FunctionSpace from \var{x}:
\begin{python}
print x0.getRank(),x0.getShape(),x0.getFunctionSpace()
\end{python}
will print
\begin{python}
0 () Function space type: Finley_Nodes on FinleyMesh 
\end{python}
We can now construct the function \var{gammaD} by
\begin{python}
gammaD=whereZero(x[0])+whereZero(x[1])
\end{python}
where
\code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (allmost) equal to zero 
and $0$ elsewhere. 
Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is 
equal to zero and $0$ elsewhere.
The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])} 
gives a function on the domain \var{mydomain} which is exactly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero.
Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from 
\begin{python}
print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace()
\end{python}
one gets 
\begin{python}
0 () Function space type: Finley_Nodes on FinleyMesh 
\end{python}
The additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the 
characteristic function \index{characteristic function} of the locations
of the domain where homogeneous Dirichlet boundary condition \index{Dirichlet boundary condition!homogeneous}
are set. The complete definition of our example is now: 
\begin{python}
from esys.linearPDEs import Poisson
x = mydomain.getX()
gammaD = whereZero(x[0])+whereZero(x[1])
mypde = Poisson(domain=mydomain)
mypde = setValue(f=1,q=gammaD)
\end{python}
The first statement imports the \Poisson class definition form the \linearPDEs module \escript package.
To get the solution of the Poisson equation defined by \var{mypde} we just have to call its
\method{getSolution}. 

Now we can write the script to solve our test problem (Remember that
lines starting with '\#' are comment lines in Python) (available as \file{poisson.py} 
in the \ExampleDirectory):
\begin{python}
from esys.escript import *
from esys.escript.linearPDEs import Poisson
from esys.finley import Rectangle
# generate domain:
mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
# define characteristic function of Gamma^D
x = mydomain.getX()
gammaD = whereZero(x[0])+whereZero(x[1])
# define PDE and get its solution u
mypde = Poisson(domain=mydomain)
mypde.setValue(f=1,q=gammaD)
u = mypde.getSolution()
# write u to an external file
saveVTK("u.xml",sol=u)
\end{python}
The last statement writes the solution tagged with the name "sol" to the external file \file{u.xml} in 
\VTK file format. \VTK is a software library
for the visualization of scientific, engineering and analytical data and is freely available
from \url{http://www.vtk.org}. There are a variaty of graphical user interfaces 
for \VTK available, for instance \mayavi which can be downloaded from \url{http://mayavi.sourceforge.net/} but is also available on most
\LINUX distributions.

\begin{figure}
\centerline{\includegraphics[width=\figwidth]{FirstStepResult.eps}}
\caption{Visualization of the Possion Equation Solution for $f=1$}
\label{fig:FirstSteps.3}
\end{figure}

You can edit the script file using your favourite text editor (or the Integrated DeveLopment Environment IDLE
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}.
1	jgs	102	% $Id$
2
3	jgs	121	\section{The First Steps}
4	jgs	102	\label{FirstSteps}
5
6			\begin{figure}
7			\centerline{\includegraphics[width=\figwidth]{FirstStepDomain}}
8			\caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}
9			\label{fig:FirstSteps.1}
10			\end{figure}
11
12	jgs	107	In this chapter we will give an introduction how to use \escript to solve
13			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}
14			is sufficient. It is helpful if the reader has some basic knowledge of PDEs \index{partial differential equation}.
15	jgs	102
16	jgs	107	The PDE \index{partial differential equation} we wish to solve is the Poisson equation \index{Poisson equation}
17	jgs	102	\begin{equation}
18			-\Delta u =f
19			\label{eq:FirstSteps.1}
20			\end{equation}
21	jgs	107	for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$
22	jgs	102	is the unit square
23			\begin{equation}
24			\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 \}
25			\label{eq:FirstSteps.1b}
26			\end{equation}
27	jgs	107	The domain is shown in \fig{fig:FirstSteps.1}.
28	jgs	102
29			$\Delta$ denotes the Laplace operator\index{Laplace operator} which is defined by
30			\begin{equation}
31			\Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}
32			\label{eq:FirstSteps.1.1}
33			\end{equation}
34	jgs	107	where, for any function $w$ and any direction $i$, $u\hackscore{,i}$
35			denotes the partial derivative \index{partial derivative} of $u$ with respect to $i$.
36	jgs	102	\footnote{Some readers
37			may be more familiar with the Laplace operator\index{Laplace operator} being written
38			as $\nabla^2$, and written in the form
39			\begin{equation*}
40	jgs	110	\nabla^2 u = \nabla^t \cdot \nabla u = \frac{\partial^2 u}{\partial x\hackscore 0^2}
41	jgs	102	+ \frac{\partial^2 u}{\partial x\hackscore 1^2}
42			\end{equation*}
43			and \eqn{eq:FirstSteps.1} as
44			\begin{equation*}
45			-\nabla^2 u = f
46			\end{equation*}
47			}
48	jgs	107	Basically, in the subindex of a function, any index to the left of the comma denotes a spatial derivative with respect
49	jgs	102	to the index. To get a more compact form we will write $w\hackscore{,ij}=(w\hackscore {,i})\hackscore{,j}$
50			which leads to
51			\begin{equation}
52			\Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii}
53			\label{eq:FirstSteps.1.1b}
54			\end{equation}
55			In some cases, and we will see examples for this in the next chapter,
56			the usage of the nested $\sum$ symbols blows up the formulas and therefore
57	jgs	107	it is convenient to use the Einstein summation convention \index{summation convention}. This
58			drops the $\sum$ sign and assumes that a summation over a repeated index is performed
59	jgs	102	("repeated index means summation"). For instance we write
60			\begin{eqnarray}
61			x\hackscore{i}y\hackscore{i}=\sum\hackscore{i=0}^2 x\hackscore{i}y\hackscore{i} \\
62			x\hackscore{i}u\hackscore{,i}=\sum\hackscore{i=0}^2 x\hackscore{i}u\hackscore{,i} \\
63			u\hackscore{,ii}=\sum\hackscore{i=0}^2 u\hackscore{,ii} \\
64	jgs	107	x\hackscore{ij}u\hackscore{i,j}=\sum\hackscore{j=0}^2\sum\hackscore{i=0}^2 x\hackscore{ij}u\hackscore{i,j} \\
65	jgs	102	\label{eq:FirstSteps.1.1c}
66			\end{eqnarray}
67			With the summation convention we can write the Poisson equation \index{Poisson equation} as
68			\begin{equation}
69			- u\hackscore{,ii} =1
70			\label{eq:FirstSteps.1.sum}
71			\end{equation}
72			On the boundary of the domain $\Omega$ the normal derivative $n\hackscore{i} u\hackscore{,i}$
73			of the solution $u$ shall be zero, ie. $u$ shall fulfill
74			the homogeneous Neumann boundary condition\index{Neumann
75			boundary condition!homogeneous}
76			\begin{equation}
77			n\hackscore{i} u\hackscore{,i}= 0 \;.
78			\label{eq:FirstSteps.2}
79			\end{equation}
80			$n=(n\hackscore{i})$ denotes the outer normal field
81			of the domain, see \fig{fig:FirstSteps.1}. Remember that we
82			are applying the Einstein summation convention \index{summation convention}, i.e
83	jgs	107	$n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} +
84			n\hackscore{1} u\hackscore{,1}$.
85	jgs	102	\footnote{Some readers may familiar with the notation
86			\begin{equation*}
87			\frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i}
88			\end{equation*}
89			for the normal derivative.}
90			The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the
91			set $\Gamma^N$ which is the top and right edge of the domain:
92			\begin{equation}
93			\Gamma^N=\{(x\hackscore 0;x\hackscore 1) \in \Omega \| x\hackscore{0}=1 \mbox{ or } x\hackscore{1}=1 \}
94			\label{eq:FirstSteps.2b}
95			\end{equation}
96	jgs	107	On the bottom and the left edge of the domain which is defined
97	jgs	102	as
98			\begin{equation}
99			\Gamma^D=\{(x\hackscore 0;x\hackscore 1) \in \Omega \| x\hackscore{0}=0 \mbox{ or } x\hackscore{1}=0 \}
100			\label{eq:FirstSteps.2c}
101			\end{equation}
102			the solution shall be identically zero:
103			\begin{equation}
104			u=0 \; .
105			\label{eq:FirstSteps.2d}
106			\end{equation}
107	jgs	107	This kind of boundary condition is called a homogeneous Dirichlet boundary condition
108	jgs	102	\index{Dirichlet boundary condition!homogeneous}. The partial differential equation in \eqn{eq:FirstSteps.1.sum} together
109			with the Neumann boundary condition \eqn{eq:FirstSteps.2} and
110			Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so
111			called boundary value
112	jgs	107	problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for
113			the unknown
114	jgs	102	function $u$.
115
116
117			\begin{figure}
118			\centerline{\includegraphics[width=\figwidth]{FirstStepMesh}}
119			\caption{Mesh of $4 \time 4$ elements on a rectangular domain. Here
120			each element is a quadrilateral and described by four nodes, namely
121			the corner points. The solution is interpolated by a bi-linear
122			polynomial.}
123			\label{fig:FirstSteps.2}
124			\end{figure}
125
126			In general the BVP\index{boundary value problem!BVP} cannot be solved analytically and numerical
127			methods have to be used construct an approximation of the solution
128			$u$. Here we will use the finite element method\index{finite element
129			method} (FEM\index{finite element
130			method!FEM}). The basic idea is to fill the domain with a
131	jgs	107	set of points called nodes. The solution is approximated by its
132	jgs	102	values on the nodes\index{finite element
133	jgs	107	method!nodes}. Moreover, the domain is subdivided into small,
134			sub-domain called elements \index{finite element
135	jgs	102	method!element}. On each element the solution is
136			represented by a polynomial of a certain degree through its values at
137			the nodes located in the element. The nodes and its connection through
138			elements is called a mesh\index{finite element
139	jgs	107	method!mesh}. \fig{fig:FirstSteps.2} shows an
140	jgs	102	example of a FEM mesh with four elements in the $x_0$ and four elements
141			in the $x_1$ direction over the unit square.
142			For more details we refer the reader to the literature, for instance
143			\Ref{Zienc,NumHand}.
144
145			\escript provides the class \Poisson to define a Poisson equation \index{Poisson equation}.
146			(We will discuss a more general form of a PDE \index{partial differential equation!PDE}
147			that can be defined through the \LinearPDE class later). The instantiation of
148			a \Poisson class object requires the specification of the domain $\Omega$. In \escript
149			the \Domain class objects are used to describe the geometry of a domain but it also
150			contains information about the discretization methods and the actual solver which is used
151			to solve the PDE. Here we are using the FEM library \finley \index{finite element
152			method}. The following statements create the \Domain object \var{mydomain} from the
153			\finley method \method{Rectangle}
154			\begin{python}
155	jgs	107	from esys.finley import Rectangle
156			mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
157	jgs	102	\end{python}
158			In this case the domain is a rectangle with the lower, left corner at point $(0,0)$ and
159			the right, upper corner at $(\var{l0},\var{l1})=(1,1)$.
160	jgs	107	The arguments \var{n0} and \var{n1} define the number of elements in $x\hackscore{0}$ and
161	jgs	102	$x\hackscore{1}$-direction respectively. For more details on \method{Rectangle} and
162			other \Domain generators within the \finley module,
163			see \Chap{CHAPTER ON FINLEY}.
164
165	jgs	107	The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and
166	jgs	102	the right hand side $f$ of the PDE to constant $1$:
167			\begin{python}
168	gross	565	from esys.escript.linearPDEs import Poisson
169	jgs	107	mypde = Poisson(mydomain)
170			mypde.setValue(f=1)
171	jgs	102	\end{python}
172			We have not specified any boundary condition but the
173			\Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann
174			boundary condition!homogeneous} defined by \eqn{eq:FirstSteps.2}. With this boundary
175			condition the BVP\index{boundary value problem!BVP} we have defined has no unique solution. In fact, with any solution $u$
176			and any constant $C$ the function $u+C$ becomes a solution as well. We have to add
177			a Dirichlet boundary condition \index{Dirichlet boundary condition}. This is done
178	jgs	107	by defining a characteristic function \index{characteristic function}
179			which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$ where Dirichlet boundary condition is set
180	jgs	102	and $0$ elsewhere. In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c},
181	gross	565	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
182			an object \var{x} which represents locations in the domain one uses
183	jgs	102	\begin{python}
184	gross	565	x=mydomain.getX() \;.
185	jgs	102	\end{python}
186	gross	565	In fact \var{x} is a \Data object which we will learn more about in Chapter~\ref{X}. At this stage we only have to know
187			that \var{x} has a
188
189	jgs	107	In the first statement, the method \method{getX} of the \Domain \var{mydomain}
190			gives access to locations
191	gross	565	in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which is
192			discussed in Chpater\ref{X} in more details. What we need to know here is that
193			\var{x} has \Rank (=number of dimensions) and a \Shape (=tuple of dimensions) which can be checked by
194			calling the \method{getRank} and \method{getShape} methods:
195			\begin{python}
196			print "rank ",x.getRank(),", shape ",x.getShape()
197			\end{python}
198			will print something like
199			\begin{python}
200			rank 1, shape (2,)
201			\end{python}
202			The \Data object also maintains type information which is represented by the
203			\FunctionSpace of the object. For instance
204			\begin{python}
205			print x.getFunctionSpace()
206			\end{python}
207			will print
208			\begin{python}
209			Function space type: Finley_Nodes on FinleyMesh
210			\end{python}
211			which tells us that the coordinates are stored on the nodes of a \finley mesh.
212			To get the $x\hackscore{0}$ coordinates of the locations we use the
213			statement
214			\begin{python}
215			x0=x[0]
216			\end{python}
217			Object \var{x0}
218			is again a \Data object now with \Rank $0$ and
219			\Shape $()$. It inherits the \FunctionSpace from \var{x}:
220			\begin{python}
221			print x0.getRank(),x0.getShape(),x0.getFunctionSpace()
222			\end{python}
223			will print
224			\begin{python}
225			0 () Function space type: Finley_Nodes on FinleyMesh
226			\end{python}
227			We can now construct the function \var{gammaD} by
228			\begin{python}
229			gammaD=whereZero(x[0])+whereZero(x[1])
230			\end{python}
231			where
232			\code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (allmost) equal to zero
233	jgs	107	and $0$ elsewhere.
234	gross	565	Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is
235	jgs	107	equal to zero and $0$ elsewhere.
236	gross	565	The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])}
237			gives a function on the domain \var{mydomain} which is exactly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero.
238			Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from
239			\begin{python}
240			print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace()
241			\end{python}
242			one gets
243			\begin{python}
244			0 () Function space type: Finley_Nodes on FinleyMesh
245			\end{python}
246	jgs	107	The additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the
247	jgs	102	characteristic function \index{characteristic function} of the locations
248			of the domain where homogeneous Dirichlet boundary condition \index{Dirichlet boundary condition!homogeneous}
249			are set. The complete definition of our example is now:
250			\begin{python}
251	jgs	107	from esys.linearPDEs import Poisson
252	jgs	102	x = mydomain.getX()
253	gross	565	gammaD = whereZero(x[0])+whereZero(x[1])
254	jgs	107	mypde = Poisson(domain=mydomain)
255			mypde = setValue(f=1,q=gammaD)
256	jgs	102	\end{python}
257	gross	565	The first statement imports the \Poisson class definition form the \linearPDEs module \escript package.
258	jgs	107	To get the solution of the Poisson equation defined by \var{mypde} we just have to call its
259	jgs	102	\method{getSolution}.
260
261			Now we can write the script to solve our test problem (Remember that
262	gross	569	lines starting with '\#' are comment lines in Python) (available as \file{poisson.py}
263	jgs	102	in the \ExampleDirectory):
264			\begin{python}
265	gross	565	from esys.escript import *
266	gross	569	from esys.escript.linearPDEs import Poisson
267	jgs	107	from esys.finley import Rectangle
268	jgs	102	# generate domain:
269	jgs	107	mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
270	jgs	102	# define characteristic function of Gamma^D
271			x = mydomain.getX()
272	gross	565	gammaD = whereZero(x[0])+whereZero(x[1])
273	jgs	102	# define PDE and get its solution u
274	gross	565	mypde = Poisson(domain=mydomain)
275			mypde.setValue(f=1,q=gammaD)
276	jgs	102	u = mypde.getSolution()
277			# write u to an external file
278	gross	565	saveVTK("u.xml",sol=u)
279	jgs	102	\end{python}
280	gross	565	The last statement writes the solution tagged with the name "sol" to the external file \file{u.xml} in
281			\VTK file format. \VTK is a software library
282	jgs	102	for the visualization of scientific, engineering and analytical data and is freely available
283	gross	565	from \url{http://www.vtk.org}. There are a variaty of graphical user interfaces
284			for \VTK available, for instance \mayavi which can be downloaded from \url{http://mayavi.sourceforge.net/} but is also available on most
285			\LINUX distributions.
286	jgs	102
287			\begin{figure}
288			\centerline{\includegraphics[width=\figwidth]{FirstStepResult.eps}}
289	gross	565	\caption{Visualization of the Possion Equation Solution for $f=1$}
290	jgs	102	\label{fig:FirstSteps.3}
291			\end{figure}
292
293	jgs	107	You can edit the script file using your favourite text editor (or the Integrated DeveLopment Environment IDLE
294	gross	569	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
295	jgs	102	script from any shell using the command:
296	gross	565	\begin{python}
297	gross	569	python poisson.py
298	gross	565	\end{python}
299	gross	569	After the script has (hopefully successfully) been completed you will find the file \file{u.xml} in the current
300	jgs	102	directory. An easy way to visualize the results is the command
301	gross	565	\begin{python}
302			mayavi -d u.xml -m SurfaceMap &
303			\end{python}
304			to show the results, see \fig{fig:FirstSteps.3}.
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision