--- a +++ b/doc/doc.tex @@ -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}