Binary file doc/HyperElasticity.png has changed

--- a/doc/doc.bib	Fri Dec 20 09:09:56 2013 +0100
+++ b/doc/doc.bib	Fri Dec 20 11:08:10 2013 +0100
@@ -76,6 +76,12 @@
  year = {2012},
 }
 
+@Article{Formaggia_matr,
+ author = {Formaggia, Luca},
+ title = {Advanced Programming for Scientific Computing. Lecture title: Storing techniques for sparse matrices},
+ year = {2012},
+}
+
 @book{logg2012automated,
   title={Automated solution of differential equations by the finite element method: The fenics book},
   author={Logg, Anders and Mardal, Kent-Andre and Wells, Garth},

Binary file doc/doc.pdf has changed

--- a/doc/doc.tex	Fri Dec 20 09:09:56 2013 +0100
+++ b/doc/doc.tex	Fri Dec 20 11:08:10 2013 +0100
@@ -1,4 +1,4 @@
-\documentclass[a4paper,10pt]{book}
+\documentclass[a4paper,11pt]{book}
 \usepackage[utf8x]{inputenc}
 \usepackage{graphicx}
 \usepackage{amsmath}
@@ -83,8 +83,12 @@
 \end{center}
 \end{titlepage}
 
+\frontmatter
+
 \tableofcontents
 
+\mainmatter
+
 \chapter{Introduction}
 Fem-fenics is an open source package (pkg) 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 
@@ -109,10 +113,15 @@
 \begin{center}
 \url{http://gedeone-gsoc.blogspot.com/}
 \end{center}
-Finally, the API is available as Appendix but also at the following address
+The API is available as Appendix \ref{app} but also at the following address
 \begin{center}
 \url{http://octave.sourceforge.net/fem-fenics/overview.html}
 \end{center}
+and if you would like to contribute to the project or give a look to the source code
+you can clone it from the following repository using mercurial
+\begin{center}
+\url{http://sourceforge.net/p/octave/fem-fenics/} .
+\end{center}
 
 \chapter{Introduction to Fem-fenics}\label{intr} 
 
@@ -394,7 +403,7 @@
   \item make the interface as simple as possible.
 \end{itemize}
 
-\section{General layout of a class}
+\section{General layout of a class}\label{class}
 Seven new classes are implemented for dealing with FEniCS objects and for using them inside Octave:
 \begin{itemize}
   \item \textbf{boundarycondition} stores and builds a dolfin::DirichletBC
@@ -870,18 +879,19 @@
 where the Finite Elements denoted with * are not yet fully supported inside FEniCS.
 
 \section{General layout of a function}
-There are three general kinds of functions in the code: functions which create an abstract problem
-(wrappers to UFL),
-functions which create the specific instance of a problem (wrapper to FEniCS) and functions
-which discretize the problem and generate the matrices.
-\subsection{Wrappers to UFL}
-As stated in section \ref{genlayout}, a problem is divided in two files:a \textit{.ufl} file 
+There are two general kinds of functions in the code: functions which create an abstract problem
+(wrappers to UFL) and
+functions which create the specific instance of a problem and discretize it (wrapper to DOLFIN).
+\section{Wrappers to UFL}
+As stated in section \ref{genlayout}, a problem is divided in two files: a \texttt{.ufl} file 
 where the abstract problem is described in Unified Form Language (UFL),
-and a script file \textit{.m} where a specific problem is implemented and solved.
-We suppose that they are called Poisson.ufl and Poisson.m .
+and a script file \texttt{.m} where a specific problem is implemented and solved.
+We suppose that they are called \texttt{Poisson.ufl} and \texttt{Poisson.m} .
 In order to use the information stored in the UFL file, i.e. the bilinear and the linear form, 
-they have to be imported inside Octave.
-When the UFL file is compiled using the ffc compiler, a header file \textit{Poisson.h} is generated.
+they have to be ``imported'' inside Octave. This is done using the 
+functions \texttt{import\_ufl\_BilinearForm, import\_ufl\_LinearForm, ...} .
+\subsection{Generation of code on the fly}
+When a UFL file is compiled using the ffc compiler, a header file \texttt{Poisson.h} is generated.
 In this header file, it is defined the Poisson class, which derives from dolfin::Form,
 and the constructor for the bilinear and linear form are set.
 This file is thus available only at compilation time, but it has to be included somehow
@@ -890,11 +900,12 @@
 change the values of the coefficient for a specific problem; 
 but, as we want to let the user free to write his own 
 variational problem, a different approach has been adopted.
-The ufl file is compiled at run time and generates a header file.
+The \texttt{ufl} file is compiled at run time and generates its header file.
 Then, a Poisson.cc file is written from a template which takes as input the name 
 of the header file and is compiled including the Poisson.h file;
-now the corresponding octave functions for the specific problem is available and will be used from
-BilinearForm, LinearForm, FunctionSpace, ... .
+now the corresponding Octave functions for the specific problem are available and
+will be later used from
+\texttt{BilinearForm, LinearForm, FunctionSpace, ...} .
 As an example it is presented the import\_ufl\_BilinearForm function.
 
 \begin{lstlisting}[language=Octave]
@@ -902,7 +913,7 @@
 
  ...
  
-  %the function which writes the var-prob.cc file
+  %the function which writes the var_prob.cc file (see below)
   generate_rhs (var_prob);
   
   %the function which writes the makefile
@@ -956,27 +967,29 @@
 \end{lstlisting}
 
 
-\subsection{Wrappers to DOLFIN}
+\section{Wrappers to DOLFIN}
 The general layout of a function is very simple and it is composed of 4 steps which we describe using an example:
 \begin{lstlisting}
-DEFUN_DLD (fem_fs, args, , "initialize a fs from a mesh")
+DEFUN_DLD (FunctionSpace, args, , "initialize a FunctionSpace 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
+          
+          // 2 extract the data stored in the Octave class as a DOLFIN object
           const dolfin::Mesh & mshd = msho.get_msh ();
-          // 3 build the new object using dolfin
+          
+          // 3 build a new object or extract the information needed 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
+          
+          // 4 convert the new object from DOLFIN to Octave and return it
           octave_value retval = new functionspace(g);
-           return retval;
+          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.
-\subsubsection{Sparse Matrices}
-\subsubsection{Polymorphism}
-\subsubsection{DirichletBC and Coefficient}
+
+\subsection{DirichletBC and Coefficient}
  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. 
  It has thus been derived from dolfin::Expression a class "expression" which has as private member 
@@ -1023,6 +1036,116 @@
 overloads the \verb$eval()$ method of the  dolfin::Expression class and it is evaluated at 
 run time using the octave function \verb$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.
+
+\subsection{Sparse Matrices}
+The \texttt{assemble} function discretize the continuos problem and
+return a sparse matrix. To deal with problems of big size, they are stored
+using a compressed technique \cite{Formaggia_matr} both in DOLFIN and in Octave.
+Unfortunately, DOLFIN uses row major orientation while Octave uses 
+column major orientation. They have thus to be converted efficiently from
+one type to the other.
+\begin{lstlisting}[language=C++]
+ 
+#include "form.h"
+#include "boundarycondition.h"
+
+DEFUN_DLD (assemble, args, nargout, " ")
+{
+  int nargin = args.length ();
+  octave_value_list retval;
+
+  if (! boundarycondition_type_loaded)
+    {
+      boundarycondition::register_type ();
+      boundarycondition_type_loaded = true;
+      mlock ();
+    }
+
+  if (! form_type_loaded)
+    {
+      form::register_type ();
+      form_type_loaded = true;
+      mlock ();
+    }
+
+  ...
+  
+  // Extract form object from the input
+  const form & frm = static_cast<const form&> (args(0).get_rep ());
+
+  const dolfin::Form & a = frm.get_form ();
+  a.check ();
+
+  ...
+  
+  // Assemble the Matrix in DOLFIN
+  dolfin::parameters["linear_algebra_backend"] = "uBLAS";
+  dolfin::Matrix A;
+  dolfin::assemble (A, a);
+
+  // Extract BC from input and apply BC
+  ...
+  const boundarycondition & bc
+    = static_cast<const boundarycondition&> (args(i).get_rep ());
+
+  const std::vector<boost::shared_ptr <const dolfin::DirichletBC> > 
+    & pbc = bc.get_bc ();
+
+  for (std::size_t j = 0; j < pbc.size (); ++j)
+    pbc[j]->apply(A);
+
+
+  // Get capacity of the dolfin sparse matrix
+  boost::tuples::tuple<const std::size_t*, 
+		       const std::size_t*, 
+	               const double*, int> 
+    aa = A.data ();
+
+  // Create the Matrix for Octave
+
+  ...
+  
+  for (std::size_t i = 0; i < nr; ++i)
+    {
+      A.getrow (i, cidx_tmp, data_tmp);
+      nz += cidx_tmp.size ();
+
+      for (octave_idx_type j = 0; 
+	  j < cidx_tmp.size (); ++j)
+	{
+	  orow [ii + j] = i;
+	  oc [ii + j] = cidx_tmp [j];
+	  ov [ii + j] = data_tmp [j];
+	}
+
+      ii = nz;
+    }
+
+  dims(0) = ii;
+  ridx.resize (dims);
+  cidx.resize (dims);
+  data.resize (dims);
+
+  SparseMatrix sm (data, ridx, cidx, nr, nc);
+  retval(0) = sm;
+
+  ...
+  
+  return retval;
+}
+\end{lstlisting}
+
+
+\subsection{Polymorphism}
+The objects which belong to the new classes presented in section \ref{class}
+have to overload some of the methods already available in Octave.
+For example, we want to be able to \texttt{plot} a \texttt{Mesh} or a \texttt{function}, to \texttt{save} it
+and to evaluate it at a specific point in the space (\texttt{feval}).
+As Octave is a dynamically typed language it could be a difficult task to achieve, but hopefully
+the Octave interpreter takes care of it and it is enough to put the polymorphic function
+in a folder named as the type. For example, in the \texttt{@function} folder
+inside the Fem-fenics directory we can find the \texttt{plot, save, feval} functions.
+
 \iffalse
 
 \paragraph{other function}
@@ -1055,20 +1178,6 @@
         This will be illustrated below in the HyperElasticity example.
 \fi
 
-\subsection{Wrapper to FEniCS}
-
-\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.
-\fi
-
-
 
 \chapter{More Advanced Examples}\label{exem}
 In this chapter more examples are provided.
@@ -1155,6 +1264,8 @@
 
 
 
+
+
 # Define essential boundary
 
 bc1 = DirichletBC (SubSpace (V, 1), @(x,y) [0; -sin(5.0*x)], 1);
@@ -1298,7 +1409,7 @@
 
 mesh = Mesh ('lshape.xml');
 
-# Define function spaces (P2-P1). From ufl file
+# Define function spaces (P2-P1). UFL file
 #  V = VectorElement("CG", triangle, 2)
 #  Q = FiniteElement("CG", triangle, 1)
 V = FunctionSpace ('VelocityUpdate', mesh);
@@ -1444,9 +1555,9 @@
 q = TestFunction(Q)
 
 # Set parameter values
+nu = 0.01
 dt = 0.01
 T = 3
-nu = 0.01
 
 # Define time-dependent pressure BC
 p_in = Expression("sin(3.0*t)", t=0.0)
@@ -1486,7 +1597,6 @@
 L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
 
 
-
 # Assemble matrices
 A1 = assemble(a1)
 A2 = assemble(a2)
@@ -1714,9 +1824,7 @@
 bcr = DirichletBC (V, Rotation, 2);
 bcs = {bcl, bcr};
 
-
 # Define source and boundary traction functions
-
 B = Constant ('B', [0.0; -0.5; 0.0]);
 T = Constant ('T', [0.1; 0.0; 0.0]);
 
@@ -1924,12 +2032,16 @@
 \end{parcolumns}
 \end{changemargin}
 
+\iffalse
 \section{Fictitious Domain}
 A penalization method to take into account obstacles in incompressible viscous flows
-
+\fi
 \newpage 
+
+\backmatter
+
 \appendix
-\chapter{API reference}
+\chapter{API reference}\label{app}
 
 \section{Import problem defined with ufl}
 \subsection*{import\_ufl\_BilinearForm}
@@ -2044,6 +2156,7 @@
     }
 
 \end{lstlisting} 
+\iffalse
 \paragraph {Warning} If in the Makefile.in you write something like
 \begin{lstlisting}[language=make]
     HAVE_DOLFIN_H = @HAVE_DOLFIN_H@  
@@ -2053,6 +2166,7 @@
     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.
+\fi
 
 \bibliographystyle{unsrt} 
 \bibliography{doc}