Binary file doc/Fem-fenics_poisson.png has changed

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/doc.bib	Sun Nov 10 16:24:36 2013 +0000
@@ -0,0 +1,36 @@
+\begin{thebibliography}{9}
+
+
+@Unpublished{ufl,
+ title = {\url{http://fenicsproject.org/documentation/ufl/1.2.0/user/user_manual.html}},
+}
+
+@Unpublished{instfenics,
+ title = {\url{http://fenicsproject.org/download/}},
+}
+
+@Unpublished{instoctave,
+ title = {\url{http://www.gnu.org/software/octave/download.html}},
+}
+
+@Unpublished{msh,
+ title = {\url{http://octave.sourceforge.net/msh/index.html}},
+}
+
+@Unpublished{fpl,
+ title = {\url{http://octave.sourceforge.net/fpl/index.html}},
+}
+
+@Unpublished{secs1d,
+ title = {\url{http://octave.sourceforge.net/secs1d/index.html}},
+}
+
+@Unpublished{bim,
+ title = {\url{http://octave.sourceforge.net/bim/index.html}},
+}
+
+@Unpublished{api,
+ title = {\url{http://octave.sourceforge.net/fem-fenics/overview.html}},
+}
+
+\end{thebibliography} 
\ No newline at end of file

Binary file doc/doc.pdf has changed

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/doc.tex	Sun Nov 10 16:24:36 2013 +0000
@@ -0,0 +1,883 @@
+\documentclass[a4paper,10pt]{book}
+\usepackage[utf8x]{inputenc}
+\usepackage{graphicx}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{caption}
+\usepackage{subcaption}
+\usepackage{cclicenses}
+\usepackage{url}
+\usepackage{listings}
+\usepackage{hyperref}
+\usepackage{listings}
+\usepackage{color}
+\usepackage[english]{babel}
+\usepackage[T1]{fontenc}
+
+
+\definecolor{dkgreen}{rgb}{0,0.6,0}
+\definecolor{gray}{rgb}{0.5,0.5,0.5}
+\definecolor{mauve}{rgb}{0.58,0,0.82}
+
+\lstset{frame=tb,
+  language=Octave,
+  aboveskip=3mm,
+  belowskip=3mm,
+  showstringspaces=false,
+  columns=flexible,
+  basicstyle={\small\ttfamily},
+  numbers=left,
+  numberstyle=\tiny\color{gray},
+  keywordstyle=\color{blue},
+  commentstyle=\color{dkgreen},
+  stringstyle=\color{mauve},
+  breaklines=true,
+  breakatwhitespace=true
+  tabsize=3
+}
+
+\author{Marco Vassallo}
+\title{ \textbf{Fem-fenics} \\ \bigskip \textit{General purpose Finite Element library \\ for GNU-Octave}\\ 
+\bigskip
+       \textsc{ Work in progress \\(help and remarks are welcome)}}
+
+
+\begin{document}
+
+\maketitle
+\tableofcontents
+
+\chapter{Introduction}
+Fem-Fenics is an open-source package for the resolution of partial differential equations with Octave.
+The project has been developed during the Google Summer of Code 2013 with the help and the sustain of the GNU-Octave 
+community under the supervision of prof. De Falco.\\
+
+The report is structured as follows:
+\begin{itemize}
+  \item in chapter \ref{intr} we provide a simple reference guide for beginners
+  \item in chapter \ref{impl} is presented a detailed explanation of the relevant parts of the program. In this way, the 
+  interested reader can see what there is ``behind'' and expecially anyone interested in it can learn quickly how
+  it is possible to extend the code and contribute to the project.
+  \item in chapter \ref{exem} more examples are provided. For a lot of them, we present the octave script 
+  alongside with the code for Fenics (in C++ and/or Python) in order to provide the user with a quick reference
+  guide.
+\end{itemize}
+
+If you think that going inside the report could be boring, it is available a wiki at
+\begin{center}
+\url{http://wiki.octave.org/Fem-fenics}
+\end{center}
+while if you want to see how the project has grown during the time you can give a look at
+\begin{center}
+\url{http://gedeone-gsoc.blogspot.com/}
+\end{center}
+Finally, the API is available at the following address
+\begin{center}
+\url{http://octave.sourceforge.net/fem-fenics/overview.html}
+\end{center}
+
+\chapter{Introduction to Fem-fenics}\label{intr} 
+
+\section{Installation}
+Fem-fenics is an external package for Octave. It means that you can install it only once that you
+have successfully installed Octave on your PC. Furthermore, as Fem-fenics is based on Fenics,
+you also need a running version of the latter. They can be easily installed following the guidelines provided
+on the official Octave \cite{instoctave} and Fenics \cite{instfenics} websites.
+Once that you have got Octave and Fenics, you can just launch Octave (which now is provided with a new
+amazing GUI) and type
+
+\begin{verbatim}
+ >> pkg install fem-fenics -forge
+\end{verbatim}
+
+That's all! If you encounter any problem during the installation don't hesitate to contact us.
+To be sure that everything is working fine, you can load the fem-fenics pkg and run
+one of the examples provided within the package:
+
+\begin{verbatim}
+ >> pkg load fem-fenics
+ >> femfenics_examples()
+\end{verbatim}
+
+For a description of the examples, you are referred to chapter \ref{exem}.
+\\
+\\
+\textbf{NOTE} For completing the installation process successfully,
+the form compiler FFC and the header file dolfin.h should also be available on your machine.
+They are managed automatically by Fenics if you install it as a binary package or with Dorsal.
+If you have done it manually, please be sure that they are available before starting the
+installation of Fem-fenics.
+
+\section{General layout}
+
+A generic problem has to be solved in two steps:
+\begin{enumerate}
+\item a \textbf{.ufl file} where the abstract problem is described: this file has to be written in Unified Form Language (UFL),
+       which is a domain specific language for defining discrete variational forms and functionals in a notation 
+       close to pen-and-paper formulation. UFL is easy to learn, and the User manual provides explanations
+       and examples \cite{ufl}.
+\item a script file \textbf{.m} where the abstract problem is imported and a specific problem is implemented and solved:
+      this is the script file where the fem-fenics functions described in the following chapters are used.
+\end{enumerate}
+
+We provide immediately a simple example in order to familiarize the user with the code.
+
+\paragraph{The Poisson equation}
+In this example, we show how it is possible to solve the Poisson equation with mixed Boundary Conditions.
+If we indicate with $\Omega$ the domain and with $\Gamma = \Gamma_{N} \cup \Gamma_{D}$ the
+boundaries, the problem can be expressed as
+\begin{align*}
+  \Delta u &= f \qquad \text{on } \Omega \\
+  u &= 0 \qquad \text{on } \Gamma_{D} \\
+  \nabla u \cdot n &= g \qquad \text{on } \Gamma_{N}
+\end{align*}
+where $f, \, g$ are data which represent the source for and the flux
+of the scalar variable $u$.
+A possible variational formulation of the problem is: \\
+find $u \in H_{0, \Gamma_{D}}^{1} :$
+\begin{align*}
+  a(u, v) &= L(v) \qquad \forall v \in H_{0, \Gamma_{D}}^{1} \\
+  a(u, v) &= \int_{\Omega} \nabla u \cdot \nabla v \\
+  L(v)    &= \int_{\Omega} f v + \int_{\Gamma_{N}} g v \\
+\end{align*}
+
+The abstract problem can thus be written in the \verb|Poisson.ufl| file immediately.
+The only thing that we have to specify at this stage is the space of Finite Elements
+which we want to use for the discretization of $H_{0, \Gamma_{D}}^{1}$. In our case,
+we choose the space of continuous lagrangian polynomial of degree one 
+\verb|FiniteElement("Lagrange", triangle, 1)|, but many more
+possibilities are available.
+\subparagraph{Poisson.ufl}
+
+\begin{lstlisting}
+element = FiniteElement("Lagrange", triangle, 1)
+
+u = TrialFunction(element)
+v = TestFunction(element)
+
+f = Coefficient(element)
+g = Coefficient(element)
+
+a = inner(grad(u), grad(v))*dx
+L = f*v*dx + g*v*ds
+\end{lstlisting}
+
+It is always a good idea to check if the ufl code is correctly written before importing it into Octave:
+\begin{verbatim}
+ >> ffc -l dolfin Poisson.ufl
+\end{verbatim} 
+shouldn't produce any error.
+We can now implement and solve a specific instance of the Poisson problem in Octave.
+We choose to set the parameters as follow
+\begin{itemize}
+ \item $\Omega = [0, 1]\times[0, 1]$
+ \item $\Gamma_{D} = {(0, y) \cup(1, y)} \ \subset \partial\Omega$
+ \item $\Gamma_{N} = {(x, 0) \cup(x, 1)} \ \subset \partial\Omega$
+ \item $f = 10 \exp \dfrac{(x-0.5)^{2} + (y-0.5)^{2}}{0.02}$
+ \item $g = \sin(5x)$
+\end{itemize}
+\iffalse
+As a first thing we need to load into Octave the pkg that we have previously installed
+
+    pkg load fem-fenics
+
+We can thus initialize the environment and import the .ufl file inside Octave. 
+When we call the function \verb$fem_init_env ()$, the new types defined inside fem-fenics 
+are registered to the octave-parser. From the .ufl file, we have to generate the 
+corresponding functions for fem-fenics. There is a specific function which seeks 
+every specific element defined inside the .ufl file:
+
+    \verb$fem_create_fs ('name')$ is a function which looks for the finite element 
+    space defined inside the file called name.ufl; if everything is ok, it generates a function
+    \verb$fem_fs_name$
+    \verb$fem_create_rhs ('name')$ is a function which looks for the rhs of the 
+    equation, i.e. for the bilinear form defined inside name.ufl; if successful, 
+    it generates the function \verb$fem_rhs_name$
+    \verb$fem_create_lhs ('name')$ is a function which looks for the linear 
+    form defined inside name.ufl; if successful, it generates the function \verb$fem_lhs_name$
+
+In some cases one could be interested in using these functions separately but if, 
+as in our example, all the three elements are defined in the same .ufl file, 
+we can just call \verb$fem_create_all ('name')$ which generates at once all the three functions described above.
+\begin{verbatim}
+    fem_init_env ();
+    fem_create_all ('Poisson');
+\end{verbatim}
+To set the concrete elements which define our problem, 
+the first things to do is to create a mesh. 
+It can be managed easily using the msh pkg. In our case, we create a uniform mesh
+\begin{verbatim}
+    x = y = linspace (0, 1, 32);
+    msho = msh2m_structured_mesh (x, y, 1, 1:4); 
+\end{verbatim}
+Once that the mesh is available, we can thus initialize the 
+fem-fenics mesh using the function \verb$fem_init_mesh()$:
+\begin{verbatim}
+    mshd = fem_init_mesh (msho);
+\end{verbatim}
+If instead of an Octave mesh you have a mesh stored in a different format,
+you can try to convert it to the dolfin xml format using the program dolfin-convert.
+In fact, \verb$fem_init_mesh()$ accepts as arguments also a string with the name of 
+the dolfin xml file which contains your mesh. For example, if you have a mesh saved
+in the gmsh format, you can do as follows:
+\begin{verbatim}
+Shell:
+    dolfin-convert msh.gmsh msh.xml
+\end{verbatim}
+and then inside our Octave script:
+\begin{verbatim}
+    mshd = fem_init_mesh ('msh.xml');
+\end{verbatim}
+To initialize the functional space, we have to specify as argument only the fem-fenics mesh,
+because the finite element type and the polynomial degree have been specified in the .ufl file:
+\begin{verbatim}
+    V  = fem_fs_Poisson (mshd);
+\end{verbatim}
+We can now apply essential BC \verb$using fem_bc ()$. This function receives as argument the functional space,
+a function handle which specifies the value that we want to set, and the label of the sides where we want
+to apply the BC.
+\begin{verbatim}
+    bc = fem_bc (V, @(x,y) 0, [2, 4]);
+\end{verbatim}
+The last thing that we have to do before solving the problem, is to set the coefficients. 
+In our case they are the source term 'f' and the normal flux 'g'. 
+To set them, we can for example use the function \verb$fem_coeff ()$ to which we have to pass as argument a string,
+which specifies the name of the coefficient as stated in the .ufl file, and a function handle with the value:
+\begin{verbatim}
+  f = fem_coeff ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
+  g = fem_coeff ('g', @(x,y) sin (5.0 * x));
+\end{verbatim}
+Another possibility for dealing with the coefficients defined in the .ufl file would be to use 
+the function \verb$fem_func ()$, but I will describe it in a future post.
+
+We can now obtain the discretized representation of our operator using the 
+functions \verb$fem_r(l)hs_Poisson ()$. We have to pass as arguments two things: 
+the BC(s) that we want to apply and the coefficient.
+In our case, 'f' and 'g' are both related to the lhs of the equation and 
+so they have to be passed only to \verb$fem_lhs_Poisson ()$. On the other side, 
+the BC(s) has to be specified both for the rhs and for the lhs because 
+they are imposed setting directly to the right values the corresponding entry:
+\begin{verbatim}
+  A = fem_rhs_Poisson (V, bc);
+  b = fem_lhs_Poisson (V, f, g, bc);
+\end{verbatim}
+Here A is a sparse matrix and b is a column vector. We can thus use all the 
+functionalities available within Octave: if we are interested in the spectrum 
+of the matrix we can use the eig() function, but for the moment we use the backslash 
+command to solve the linear system:
+\begin{verbatim}
+  u = A \ b;
+\end{verbatim}
+Once that the solution has been obtained, we can transform the vector into a 
+fem-fenics function and plot it using \verb$fem_plot ()$, or save it with \verb$fem_save ()$.
+\begin{verbatim}
+  func = fem_func ('u', V, u)
+  fem_plot (func);
+\end{verbatim}
+At the moment \verb$fem_plot ()$ takes as input only arguments of type fem-function, 
+but we will extend it later so that it could also plot variables of type fem-mesh (and maybe even more).
+The result that we obtain should be something like this
+\fi
+
+A possible implementation for the Poisson problem is reported below, while
+in figure \ref{Poissonfig} is presented the output.
+\begin{figure}
+ \begin{center}
+  \includegraphics[height=7 cm,keepaspectratio=true]{./Fem-fenics_poisson.png}
+   \caption{The result for the Poisson equation}
+   \label{Poissonfig}
+  \end{center}
+\end{figure}
+
+\subparagraph{Poisson.m}
+\begin{lstlisting}
+#load the pkg and import the ufl problem
+pkg load fem-fenics msh
+import_ufl_Problem ('Poisson')
+ 
+# Create the mesh and define function space
+x = y = linspace (0, 1, 33);
+mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
+V = FunctionSpace('Poisson', mesh);
+
+# Define boundary condition and source term
+bc = DirichletBC(V, @(x, y) 0.0, [2;4]);
+f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
+g = Expression ('g', @(x,y) sin (5.0 * x));
+
+#Create the Bilinear and the Linear form
+a = BilinearForm ('Poisson', V, V);
+L = LinearForm ('Poisson', V, f, g);
+ 
+#Extract the matrix and compute the solution
+[A, b] = assemble_system (a, L, bc);
+sol = A \ b;
+u = Function ('u', V, sol);
+ 
+# Save solution in VTK format and plot it
+save (u, 'poisson')
+plot (u);
+
+\end{lstlisting}
+
+
+
+\chapter{Implementation}\label{impl} 
+Two main ideas have guided us throughout the realization of the pkg:
+
+\begin{itemize}
+  \item keep the syntax as close as possible to the original one in Fenics
+  \item make the interface as simple as possible.
+\end{itemize}
+\iffalse
+Fem-fenics aims to fill a gap in Octave: even if there are pkgs for the creation of mesh \cite{msh},
+for the postprocessing of data \cite{fpl} and for the resolution of some specific pde \cite{secs1d} \cite{bim},
+no general purpose finite element library is available.
+Our goal is thus to provide a pkg which can be used to solve user defined problem and which is able to
+exploit the functionality provided from Octave. Instead of writing a library from scratch, we decided to
+write an interface for one of the finite element library which are already available.
+We considered a lot of possibilities: Dune Elmer Fenics FreeFem++ Hermes LifeV
+
+``The FEniCS Project is a collection of free, open source, software
+components with the common goal to enable automated solution of pde.''
+
+Dolfin is a C++/Python interface of FEniCS, providing a consistent Problem
+Solving Environment for ODE and PDE.
+Why DOLFIN?
+C++ library
+Provides data structures and algorithms for computational meshes
+and finite element assembly (... and also a namespace)
+A lot of examples and broadly supported
+The code is well written/organized
+
+
+Our idea is to create wrappers in Octave for Dolfin, which is the problem-solving C++ environment of Fenics ,
+in a similar way to what it has been done for Python. 
+
+The general structure of the code is summarized in figure . The start point is a.
+\fi
+\section{General layout of a class}
+All these classes derive from \verb$octave_base_value$. 
+\iffalse
+By now, the general layout of a class is very simple and looks like:
+As a general rule, now all the private members are handled using a \verb$boost::shared_ptr< >$ to the corresponding DOLFIN type.  
+This is done for two main reasons:
+    in our program we have to refer in several places to resources which are built dynamically 
+    (e.g. a FunctionSpace contain a reference to the Mesh) and we want that it is destroyed only
+    when the last references is destroyed.
+    they are widely used inside DOLFIN and it is thus easier to deal with it if we use the same types.
+All the classes implement a constructor which takes as argument the corresponding dolfin type and a default 
+constructor which is necessary to register a type in the Octave interpreter.
+
+\paragraph{mesh}  The mesh class derives publicly from \verb$octave_base_value$ and dolfin::Mesh and it provides 
+the minimum of functionalities needed for our TDD.
+In addition to usual methods, we have implemented functionalities which allow us to deal with meshes 
+currently available with the msh pkg. We have thus added the constructor
+\newline
+ \verb$mesh (Array<double>& p, Array<octave_idx_type>& e, Array<octave_idx_type>& t);$
+\newline
+and the method
+
+  \verb$octave_scalar_map get_pet (void) const;$
+
+We have thus implement the functions which allow us to initialize an object of this class from Octave
+\newline
+  \verb$  DEFUN_DLD (fem_init_mesh, args, , "...")$
+\newline
+and to get it back in the (p, e, t) format
+\newline
+  \verb$  DEFUN_DLD (fem_get_mesh, args, ,"... ")$
+\newline
+We can thus now run in Octave a script like:
+\begin{verbatim} 
+  msh1 = msh2m_structured_mesh(1:6,1:4,0,3:6);
+  msh2 = fem_init_mesh (msh1);
+  msh3 = fem_get_mesh (msh2);
+\end{verbatim}
+We can thus check that the msh1 and msh3 variables represent the same mesh and with the command whos we get
+\begin{verbatim}
+   Name       Size                     Bytes      Class
+--------------------------------------------------------------------------
+    msh1        1x1                       2240      struct
+    msh2        0x0                          0      mesh
+    msh3        1x1                       2240      struct
+\end{verbatim}
+ which confirms that msh2 is a variable of our mesh class.
+
+
+\paragraph{functionspace}
+A dolfin::FunctionSpace is defined by specifying a mesh and the type of the finite element which we want to use. 
+The mesh is handled by our class, while the FE are specified inside the .ufl file. Possible choices are:
+
+    Argyris                                  ARG
+    Arnold–Winther                    AW
+    Brezzi–Douglas–Marini        BDM
+    Crouzeix–Raviart                  CR
+    Discontinuous Lagrange      DG
+    Hermite                                HER
+    Lagrange                              CG
+    Mardal–Tai–Winther             MTW
+    Morley                                  MOR
+    Nédélec 1st kind H (curl)     N1curl
+    Nédélec 2nd kind H (curl)    N2curl
+    Raviart–Thomas                   RT
+
+Once that the .ufl file has been compiled, in the output file it is defined a namespace where we can find 
+the corresponding FunctionSpace and which only needs a mesh to be initialized.
+
+\paragraph{boundary condition}
+This class is used for dealing with essential BC. We need three elements for the definition of an object of
+type dolfin::DirichletBC :
+
+    a dolfin::FunctionSpace where the BC is defined, which is handled by our class functionspace described above
+    a dolfin::Function or a dolfin::Expression which represents the value that we want to assign.
+    This is done deriving a new class from dolfin::Expression and overloading the eval function.
+    a dolfin::subdomain which specifies where we want to apply the BC. 
+\paragraph{coefficient}
+This class is used to store a std::string where we save the name of the coefficient, 
+and an object of type expression, which is a class derived from dolfin::Expression. 
+The class expression would be explained in more detail in the next post.
+\fi
+\section{General layout of a function}
+\iffalse
+The functions which we describe below are used to set and get the values for the variable that are commonly 
+used in a FEM problem.
+
+The main functions which are available right now are:
+
+    \verb$fem_init_mesh$: this function can take as input either a string which specify 
+    the name of the file to read or a mesh produced with the msh pkg
+    \verb$fem_bc$: take as input the functional space, a function handle which specifies 
+    the value that we want to set on the boundary, and the label of the boundaries where we want to apply the bc
+    coefficient: take as input the name of one of the coefficients defined in 
+    the .ufl file and a function handle with the value that we want to set
+
+The following functions work fine but at the moment are directly compiled including 
+the output file obtained from the .ufl file (which should be avoided as discussed below) :
+
+    \verb$fem_fs$: take as input the mesh where we want to define the functional space
+    \verb$fem_rhs$: take as input the function space where the bilinear form is defined, 
+    a list of coefficients and a list of BCs that we want to apply. It returns a 
+    sparse matrix which contains the discretization of the bilinear operator
+    \verb$fem_lhs$: take the same input as the function \verb$fem_rhs$. If you have essential BC 
+    in the problem, you need to pass them also to this function in such a way that 
+    the DOF corresponding to essential nodes are set to the right value
+
+The general layout of a function is very simple and is composed of 4 steps which we describe using an example:
+\begin{lstlisting}
+DEFUN_DLD (fem_fs, args, , "initialize a fs from a mesh")
+{
+          // 1 read data
+          const mesh & msho = static_cast<const mesh&> (args(0).get_rep ());
+          // 2 convert the data from octave to dolfin
+          const dolfin::Mesh & mshd = msho.get_msh ();
+          // 3 build the new object using dolfin
+          boost::shared_ptr <const dolfin::FunctionSpace> g (new Laplace::FunctionSpace (mshd));
+          // 4 convert the new object from dolfin to Octave and return it
+          octave_value retval = new functionspace(g);
+           return retval;
+}
+\end{lstlisting}
+All the functions presented above follow this general structure, and thus here we present 
+in detail only functions which present some differences.
+
+\paragraph{fem bc and fem coeff}
+ These two functions take as input a function handle which cannot be directly evaluated by
+ a dolfin function to set, respectively, the value on the boundary or the value of the coefficient. 
+ We have thus derived from dolfin::Expression a class "expression" which has as private member 
+ an octave function handle and which  overloads the function eval(). In this way, an object of 
+ the class expression can be initialized throughout a function handle and can be used inside dolfin because 
+ "it is"  a dolfin::Expression.
+\begin{lstlisting}
+class expression : public dolfin::Expression
+{
+
+ public:
+  expression (octave_fcn_handle & _f) :
+             dolfin::Expression (), f (new octave_fcn_handle (_f)) {}
+
+  void eval (dolfin::Array<double>& values, const dolfin::Array<double>& x) const
+  {
+    octave_value_list b;
+    b.resize (x.size ());
+    for (int i = 0; i < x.size (); ++i)
+      b(i) = x[i];
+    octave_value_list res = feval (f->function_value (), b);
+    values[0] = res(0).double_value ();
+  }
+
+ private:
+  octave_fcn_handle * f;
+}; 
+\end{lstlisting}
+
+\paragraph{fem lhs and fem rhs}
+ For these functions, we need to set the possible coefficient of the (bi)linear form  and to apply the BC.
+
+For the coefficient: 
+\begin{lstlisting}
+    // read the coefficient 
+    const coefficient & cf = static_cast <const coefficient&> (args (i).get_rep ());
+    // get the correct position of the coefficient in the declaration of the problem 
+    std::size_t n = a.coefficient_number (cf.get_str ());
+    // extract the information.. we use again the expression class descripted above 
+    const boost::shared_ptr<const expression> & pexp = cf.get_expr (); 
+    //set the coefficient for the (bi)linear form 
+    a.set_coefficient (n, pexp);
+\end{lstlisting}
+For the BC
+At the moment we are simply using the method apply() of the class DirichletBC, 
+but it does not preserve the symmetry of the matrix. We have thus planned to insert 
+later a function which instead of the method apply() calls the function \verb$assemble_system()$ , 
+which preserves the symmetry but which builds together the lhs and the rhs.
+The sparse matrix is then built from the dolfin::Matrix following the general guidelines of Octave.
+
+
+\paragraph{other function}
+
+
+    SubSpace allows to extract a subspace from a vectorial one. 
+    For example, if our space is P2 x P0 we can extract the one or 
+    the other and then apply BC only where it is necessary.
+    \verb$fem_eval$ takes as input a Function and a coordinate and returns a 
+    vector representing the value of the function at this point.
+    for dealing with form of rank 0, i.e. with functional, we have now 
+    added the functions \verb$fem_create_functional$ to create it from a .ufl file. 
+    We have thus extended the function assemble which returns the corresponding double value.
+    \verb$plot_2d$ and \verb$plot_3d$: these functions allow us to plot a function specifying 
+    a mesh and the value of the function at every node of the mesh. 
+    This is something which could be useful also outside of fem-fenics.
+
+
+\section{Implementation Details}
+The relevant implementation details which the user should know are:
+
+    all the objects are managed using \verb$boost::shared_ptr <>$. 
+    It means that the same resource can be shared by more objects and useless copies 
+    should be avoided. For example, if we have two different functional spaces in the same problem, 
+    like with Navier-Stokes for the velocity and the pressure, the mesh is shared between 
+    them and no one has its own copy.
+    The BC are imposed directly to the mesh setting to zero all the off diagonal elements 
+    in the corresponding line. This means that we could loose the symmetry of the matrix, if any. 
+    The simmetry is lost because we are using the method apply() of the class dolfin::DirichletBC. 
+    We have thus planned to insert later a function which instead of the method apply() calls the 
+    function \verb$assemble_system()$ , which preserves the symmetry of the system but which builds 
+    together the lhs and the rhs.
+    The coefficient of the variational problem can be specified using either a fem-coefficient 
+    or a fem-function. They are different objects which behave in different ways: a fem-coefficient 
+    object overloads the eval() method of the  dolfin::Expression class and it is evaluated at 
+    run time using the octave function feval (). A fem-function instead doesn't need to be evaluated 
+    because it is assembled copying element-by-element the values contained in the input vector.
+    
+    The relevant implementation details which the user should know are:
+
+    all the objects are managed using \verb$boost::shared_ptr <>$. 
+    It means that the same resource can be shared by more objects and useless copies should be avoided. 
+    For example, if we have two different functional spaces in the same problem, like with Navier-Stokes f
+    or the velocity and the pressure, the mesh is shared between them and no one has its own copy. 
+
+    The essential BC are imposed directly to the matrix with the command assemble(), 
+    which sets to zero all the off diagonal elements in the corresponding line, sets to 1 
+    the diagonal element and sets to the exact value the rhs. This means that we could loose
+    the symmetry of the matrix, if any. To avoid this problem and preserve the symmetry of 
+    the system it is available the \verb$assemble_system()$ command which builds at once the lhs and the rhs. 
+
+    The coefficient of the variational problem can be specified using either 
+    an Expression(), a Constant() or a Function(). They are different objects 
+    which behave in different ways: an Expession or a Constant object overloads 
+    the eval() method of the dolfin::Expression class and it is evaluated at run 
+    time using the octave function feval (). A Function instead doesn't need to 
+    be evaluated because it is assembled copying element-by-element the values 
+    contained in the input vector. 
+    
+    
+     We have split the construction of the form into two steps:
+
+        We set all the coefficients of the form using the function which we create on the fly. 
+        They will be named \verb$ProblemName_BilinearForm$ or \verb$ProblemName_LinearForm$.
+        Then we apply specific BC to the form using the assemble() function and we get back the matrix. 
+        If we are assembling the whole system and we want to keep the symmetry of the matrix (if any), 
+        we can instead use the command \verb$assemble_system$ (). Finally, if we are solving a non-linear problem 
+        and we need to apply essential BC, we should provide to the function also the vector with the 
+        tentative solution in order to modify the entries corresponding to the boundary values. 
+        This will be illustrated below in the HyperElasticity example.
+\fi
+
+\subsection{Mesh generation and conversion}
+
+\subsection{Sparse Matrices}
+
+\subsection{Shared pointer}
+
+\subsection{Polymorphism}
+
+\subsection{Code release}
+
+\subsection{Code on the fly}
+
+\iffalse
+For the creation on the fly of the code from the header file, 
+Juan Pablo provided me a python code which makes a great job. I have spent some 
+time adapting it to my problem, but when I finally got a working code, we realized 
+that it was probably enough to use the functions available inside Octave because 
+the problem was rather simple. The pyton code is however available here, while the 
+final solution adopted can be found there.
+For the problem related with \verb$fem_init_env ()$, we are still working on it. 
+\fi
+
+\subsection{Autoconf}
+ In this section we want to discuss how we can write a config.ac and a Makefile.in files which:
+\begin{itemize}
+    \item check if a program is available and stop if it is not
+    \item check if a header file is available and issue a warning if not, but go ahead with the compilation
+\end{itemize}
+We want to speak about it because, even if it is not strictly related to the fem-fenics library, 
+I hope it could be helpful for someone else because some solutions which could seem right at a 
+first sight are definitely wrong.
+
+As stated above, if we want to generate automatically a Makefile, we need two components:
+
+\paragraph{configure.ac} Is a file which checks whether the program/header is available or not 
+and sets consequently the values of some variables.
+\begin{lstlisting}
+    # Checks if the program mkoctfile is available and sets the variable HAVE_MKOCTFILE consequently
+    AC_CHECK_PROG([HAVE_MKOCTFILE], [mkoctfile], [yes], [no])
+    # if mkoctfile is not available, it issues an error and stops the compilation
+    if [test $HAVE_MKOCTFILE = "no"]; then 
+      AC_MSG_ERROR([mkoctfile required to install $PACKAGE_NAME])
+    fi
+
+    #Checks if the header dolfin.h is available; if it is available, the value of the ac_dolfin_cpp_flags is substituted with -DHAVE_DOLFIN_H, otherwise it is left empty and a warning message is printed
+    AC_CHECK_HEADER([dolfin.h],
+      [AC_SUBST(ac_dolfin_cpp_flags,-DHAVE_DOLFIN_H)  AC_SUBST(ac_dolfin_ld_flags,-ldolfin)],
+      [AC_MSG_WARN([dolfin headers could not be found, some functionalities will be disabled, don't worry your package will still be working, though.])] ).
+
+    # It generates the Makefile, using the template described below
+    AC_CONFIG_FILES([Makefile])
+\end{lstlisting} 
+\paragraph{Makefile.ac} This file is a template for the Makefile, which will be automatically generated when the configure.ac 
+file is executed. The values of the variable \verb$ac_dolfin_cpp_flags$ and \verb$ac_dolfin_ld_flags$ are substituted with the 
+results obtained above:
+\begin{lstlisting}
+    CPPFLAGS += @ac_dolfin_cpp_flags@
+    LDFLAGS += @ac_dolfin_ld_flags@
+\end{lstlisting} 
+In this way, if dolfin.h is available, CPPFLAGS contains also the flag  \verb$-DHAVE_DOLFIN_H.$
+
+\paragraph {program.cc}  Our .cc program, should thus include the header dolfin.h only if 
+\verb$-DHAVE_DOLFIN_H$ is defined at compilation time.
+For example
+
+\begin{lstlisting}
+    #ifdef HAVE_DOLFIN_H
+    #include <dolfin.h> 
+    #endif
+    int main ()
+    {  
+
+    #ifndef HAVE_DOLFIN_H
+        error("program: the program was built without support for dolfin");
+    #else 
+      /* Body of your function */
+    #endif
+     return 0;
+    }
+
+\end{lstlisting} 
+\paragraph {Warning} If in the Makefile.in you write something like
+\begin{lstlisting} 
+    HAVE_DOLFIN_H = @HAVE_DOLFIN_H@  
+    ifdef HAVE_DOLFIN_H   
+      CPPFLAGS += -DHAVE_DOLFIN_H  
+      LIBS += -ldolfin
+    endif
+ \end{lstlisting} 
+ it doesn't work because the variable \verb$HAVE_DOLFIN_H$ seems to be always defined, even if the header is not available.
+
+\chapter{More Advanced Examples}\label{exem}
+\iffalse
+ With the following examples, we can see directly in action the new features and understand how they work.
+
+    Navier-Stokes: we learn how to deal with a vector-field problem and how we can save the solution using the 
+    \verb$fem_save$ () function. We also use the fem pkg to generate a mesh using gmesh.
+    Mixed-Poisson: we solve the Poisson problem presented in the previous posts using a mixed formulation, 
+    and we see how we can extract a scalar field from a vector one.
+    HyperElasticity: we exploit the fsolve () command to solve a non-linear problem. In particular, 
+    we see how to use the assemble() function to apply BC also in this situation.
+    Advection-Diffusion: we solve a time dependent problem using the lsode () command and save 
+    the solution using the pkg flp.
+
+For each problem, we refer the reader to the complete desciption on FEniCS or bim web-page, 
+while here we highlight only the  implementation detail relevant for our pkg.
+\fi
+\section{Navier-Stokes equation with Chorin-Temam projection algorithm}
+
+\paragraph{TentativeVelocity.ufl}
+\begin{lstlisting}
+# Copyright (C) 2010 Anders Logg
+# Define function spaces (P2-P1)
+V = VectorElement("CG", triangle, 2)
+Q = FiniteElement("CG", triangle, 1)
+
+# Define trial and test functions
+u = TrialFunction(V)
+v = TestFunction(V)
+
+# Define coefficients
+k  = Constant(triangle)
+u0 = Coefficient(V)
+f  = Coefficient(V)
+nu = 0.01
+
+# Define bilinear and linear forms
+eq = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \
+    nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
+a  = lhs(eq)
+L  = rhs(eq) 
+ 
+\end{lstlisting}
+
+\paragraph{PressureUpdate.ufl}
+\begin{lstlisting}
+ # Copyright (C) 2010 Anders Logg
+ # Define function spaces (P2-P1)
+V = VectorElement("CG", triangle, 2)
+Q = FiniteElement("CG", triangle, 1)
+
+# Define trial and test functions
+p = TrialFunction(Q)
+q = TestFunction(Q)
+
+# Define coefficients
+k  = Constant(triangle)
+u1 = Coefficient(V)
+
+# Define bilinear and linear forms
+a = inner(grad(p), grad(q))*dx
+L = -(1/k)*div(u1)*q*dx 
+ 
+\end{lstlisting}
+
+\paragraph{VelocityUpdate.ufl}
+\begin{lstlisting}
+ # Copyright (C) 2010 Anders Logg
+# Define function spaces (P2-P1)
+V = VectorElement("CG", triangle, 2)
+Q = FiniteElement("CG", triangle, 1)
+
+# Define trial and test functions
+u = TrialFunction(V)
+v = TestFunction(V)
+
+# Define coefficients
+k  = Constant(triangle)
+u1 = Coefficient(V)
+p1 = Coefficient(Q)
+
+# Define bilinear and linear forms
+a = inner(u, v)*dx
+L = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
+\end{lstlisting}
+
+\paragraph{NS.m}
+\begin{lstlisting}
+pkg load fem-fenics msh
+import_ufl_Problem ("TentativeVelocity");
+import_ufl_Problem ("VelocityUpdate");
+import_ufl_Problem ("PressureUpdate");
+ 
+# We can either load the mesh from the file as in Dolfin but 
+# we can also use the msh pkg to generate the L-shape domain
+L-shape-domain;
+mesh = Mesh (msho);
+ 
+# Define function spaces (P2-P1).
+V = FunctionSpace ('VelocityUpdate', mesh);
+Q = FunctionSpace ('PressureUpdate', mesh);
+ 
+# Set parameter values and define coefficients
+dt = 0.01;
+T = 3.;
+k = Constant ('k', dt);
+f = Constant ('f', [0; 0]);
+u0 = Expression ('u0', @(x,y) [0; 0]);
+
+# Define boundary conditions
+noslip = DirichletBC (V, @(x,y) [0; 0], [3, 4]);
+outflow = DirichletBC (Q, @(x,y) 0, 2);
+
+# Assemble matrices
+a1 = BilinearForm ('TentativeVelocity', V, V, k);
+a2 = BilinearForm ('PressureUpdate', Q, Q);
+a3 = BilinearForm ('VelocityUpdate', V, V);
+A1 = assemble (a1, noslip);
+A3 = assemble (a3, noslip);
+
+# Time-stepping
+t = dt; i = 0;
+while t < T
+ 
+  # Update pressure boundary condition
+  inflow = DirichletBC (Q, @(x,y) sin(3.0*t), 1);
+ 
+  # Compute tentative velocity step
+  L1 = LinearForm ('TentativeVelocity', V, k, u0, f);
+  b1 = assemble (L1, noslip);
+  utmp = A1 \ b1;
+  u1 = Function ('u1', V, utmp);
+ 
+  # Pressure correction
+  L2 = LinearForm ('PressureUpdate', Q, u1, k);
+  [A2, b2] = assemble_system (a2, L2, inflow, outflow);
+  ptmp = A2 \ b2;
+  p1 = Function ('p1', Q, ptmp);
+ 
+  # Velocity correction
+  L3 = LinearForm ('VelocityUpdate', V, k, u1, p1);
+  b3 = assemble (L3, noslip);
+  ut = A3 \ b3;
+  u1 = Function ('u0', V, ut);
+  
+  # Save to file
+  save (p1, sprintf ("p_%3.3d", ++i));
+  save (u1, sprintf ("u_%3.3d", i));
+ 
+  # Move to next time step
+  u0 = u1;
+  t += dt
+ 
+end
+\end{lstlisting}
+
+\paragraph{L-shape-domain.m}
+\begin{lstlisting}
+name = [tmpnam ".geo"];
+fid = fopen (name, "w");
+fputs (fid,"Point (1)  = {0, 0, 0, 0.1};\n");
+fputs (fid,"Point (2)  = {1, 0, 0, 0.1};\n");
+fputs (fid,"Point (3)  = {1, 0.5, 0, 0.1};\n");
+fputs (fid,"Point (4)  = {0.5, 0.5, 0, 0.1};\n");
+fputs (fid,"Point (5) = {0.5, 1, 0, 0.1};\n");
+fputs (fid,"Point (6) = {0, 1, 0,0.1};\n");
+ 
+fputs (fid,"Line (1)  = {5, 6};\n");
+fputs (fid,"Line (2) = {2, 3};\n");
+ 
+fputs (fid,"Line(3) = {6,1,2};\n");
+fputs (fid,"Line(4) = {5,4,3};\n");
+fputs (fid,"Line Loop(7) = {3,2,-4,1};\n");
+fputs (fid,"Plane Surface(8) = {7};\n");
+fclose (fid);
+msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),...
+                   "scale", 1,"clscale", .2);
+unlink (canonicalize_file_name (name));
+\end{lstlisting}
+
+\section{A penalization method to take into account obstacles in incompressible viscous flows}
+
+
+\newpage 
+
+\bibliographystyle{unsrt} 
+\bibliography{doc}
+\end{document}