--- a/doc/doc.tex +++ b/doc/doc.tex @@ -79,19 +79,19 @@ \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 +Fem-fenics is an external package for Octave, which means that it can be installed only once that Octave has been +successfully installed on the PC. Furthermore, as Fem-fenics is based on Fenics, +it is also needed 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 +Once that Octave and Fenics are correctly installed, to install Fem-fenics 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 +That's all! For any problem during the installation don't hesitate to contact us. +To be sure that everything is working fine, load the fem-fenics pkg and run one of the examples provided within the package: \begin{verbatim} @@ -99,13 +99,13 @@ >> femfenics_examples() \end{verbatim} -For a description of the examples, you are referred to chapter \ref{exem}. +For a description of the examples, look at 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 +the form compiler FFC and the header file dolfin.h should also be available on the machine. +They are managed automatically by Fenics if it is installed as a binary package or with Dorsal. +If it has been done manually, please be sure that they are available before starting the installation of Fem-fenics. \section{General layout and first example} @@ -142,8 +142,8 @@ \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, +The only thing that has to be specified at this stage is the space of Finite Elements + used for the discretization of $H_{0, \Gamma_{D}}^{1}$. In this case, we choose the space of continuous lagrangian polynomial of degree one \verb|FiniteElement("Lagrange", triangle, 1)|, but many more possibilities are available. @@ -168,7 +168,7 @@ \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 +The parameters are setted as follow \begin{itemize} \item $\Omega = [0, 1]\times[0, 1]$ \item $\Gamma_{D} = {(0, y) \cup(1, y)} \ \subset \partial\Omega$ @@ -177,13 +177,12 @@ \item $g = \sin(5x)$ \end{itemize} -As a first thing we need to load into Octave the pkg that we have previously installed +As a first thing we need to load into Octave the pkgs previously installed \begin{verbatim} pkg load fem-fenics msh \end{verbatim} -We can thus import the ufl file inside Octave. 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: +The ufl file can thus be imported inside Octave. For every specific element defined inside the ufl file +there is a specific function which stores it for later use \begin{itemize} \item \verb$ufl_import_FunctionSpace ('Poisson')$ is a function which looks for the finite element space defined inside the file called Laplace.ufl; if everything is ok, it generates a function @@ -193,85 +192,86 @@ \item \verb$ufl_import_LinearForm ('Poisson')$ is a function which looks for the linear form. \end{itemize} + 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 (and only in this case), -we can just call \verb$import_ufl_Problem ('Poisson')$ which generates at once all -the three functions described above. +the \verb$import_ufl_Problem ('Poisson')$ can be used; it generates at once all +the three functions described above \begin{verbatim} ufl_import_Problem ('Poisson'); \end{verbatim} -To set the concrete elements which define our problem, +To set the concrete elements which define the 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 +It can be managed easily using the msh pkg. For a uniform structured mesh \begin{verbatim} x = y = linspace (0, 1, 33); 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()$: +fem-fenics mesh using the function \verb$Mesh ()$: \begin{verbatim} mesh = 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$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: +If instead of an Octave mesh we had a mesh stored in a different format, +the program dolfin-convert can try to convert it to the dolfin xml format. +Furthermore, \verb$Mesh ()$ accepts as arguments also a string with the name of +the dolfin xml file which contains the mesh. For example, for a mesh saved +in the gmsh format, do as follows: \begin{verbatim} Shell: dolfin-convert msh.gmsh msh.xml \end{verbatim} -and then inside our Octave script: +and then inside the 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: +because the finite element type and the polynomial degree have yet been specified in the ufl file: \begin{verbatim} V = FunctionSpace('Poisson', mesh); \end{verbatim} -We can now apply essential BC using \verb$DirichletBC ()$. 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. In our case, we apply homogenous boundary condition on the left and right side of the square -\begin{verbatim} - bc = DirichletBC(V, @(x, y) 0.0, [2;4]); -\end{verbatim} -The last thing that we have to do before solving the problem, is to set the coefficients specified -in the ufl file. It is important that they are called in the same way as in the ufl file. -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$Expression ()$ to which we have to pass as argument a string, +Essential BC can now be applied using \verb$DirichletBC ()$; this function receives as argument the functional space, +a function handle which specifies the value to set, and the label of the sides where the BC applies. +In this case, homogenous boundary conditions hold on the left and right side of the square +\begin{verbatim} + bc = DirichletBC(V, @(x, y) 0.0, [2; 4]); +\end{verbatim} +The last thing to do before solving the problem, is to set the coefficients specified +in the ufl file. It is important that they are called in the same way as in the ufl file: +the source term 'f' and the normal flux 'g'. +To set them, the function \verb$Expression ()$ can be used passing as argument a string, which specifies the name of the coefficient, and a function handle with the value required: \begin{verbatim} 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)); \end{verbatim} -Another possibility for dealing with the coefficients defined in the ufl file would be to use +Another possibility for dealing with the coefficients defined in the ufl file would have been to use the function \verb$Constant ()$ or \verb$Function ()$. -The coefficient can thus be used together with the FunctionSpace to set +The coefficients can thus be used together with the FunctionSpace to set the Bilinear and the Linear form \begin{verbatim} a = BilinearForm ('Poisson', V, V); L = LinearForm ('Poisson', V, f, g); \end{verbatim} -We can now obtain the discretized representation of our operator using the -functions \verb$assemble ()$ or \verb$assemble_system ()$, which also allow us -to specify the BC(s) that we want to apply. +The discretized representation of our operator is obtained using the +functions \verb$assemble ()$ or \verb$assemble_system ()$, which also allow +to specify the BC(s) to apply. Whenever possible, it is better to use the \verb$assemble_system ()$ function because it keeps the symmetry of the matrix while setting the entries for the BC: \begin{verbatim} [A, b] = assemble_system (a, L, 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 to solve the linear system. -For the moment we use the easisest possibility, i.e. the backslash command to solve the linear system: +Here A is a sparse matrix and b is a column vector. All the +functionalities available within Octave can now be exploited to solve the linear system. +The easisest possibility is the backslash command: \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 \verb$plot ()$, or save it \verb$save ()$ in the vtu +Once that the solution has been obtained, the vector is converted into a +fem-fenics function and plotted \verb$plot ()$ or saved \verb$save ()$ in the vtu format. \begin{verbatim} u = Function ('u', V, sol); @@ -280,7 +280,7 @@ \end{verbatim} The complete code for the Poisson problem is reported below, while -in figure \ref{Poissonfig} is presented the output. +in figure \ref{Poissonfig} the output is presented. \begin{figure} \begin{center} \includegraphics[height=7 cm,keepaspectratio=true]{./Fem-fenics_poisson.png} @@ -323,88 +323,236 @@ \chapter{Implementation}\label{impl} + +Fem-fenics aims to fill a gap in Octave: even if there are packages 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. + +The goal of the project is thus to provide a package which can be used to solve user defined problems +and which is able to exploit the functionality provided with Octave. + +\begin{figure} + \begin{center} + \includegraphics[height=10 cm,keepaspectratio=true]{./code_layout.png} + \caption{General layout of the package} + \label{Codelayout} + \end{center} +\end{figure} + +Instead of writing a library from scratch, an interface to one of the finite element library +which are already available has been created. +Among the many libraries taken into account, the one whch was best suited for our +purposes seemed to be the FEniCS project. It ``is a collection of free, open source, software +components with the common goal to enable automated solution of pde.'' +In particular, Dolfin is the C++/Python interface of FEniCS, providing a consistent Problem +Solving Environment for ODE and PDE. The idea has been to create wrappers in Octave for Dolfin, +in a similar way to what it has been done for Python. +This is a very natural choice, because Octave is mainly written in script language +and in C++. It is in fact possible to implement an Octave interpreter function in C++ through the +native oct-file interface or, coversely, to use Octave's Matrix/Array Classes in a C++ application. + +The works can be summirized as follows (fig. \ref{Codelayout}): + +the elements already available in Octave for the resolution of PDE (Mesh and Linear +Algebra) have been exploited, and wrappers to FEniCS functions added. +To allow exchanges between this programs, the necessary functions +for converting an Octave mesh/matrix into a FEniCS one and viceversa have been written. + 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. +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 + \item \textbf{coefficient} stores an expression object which is used for the + evaluation of user defined values + \item \textbf{expression} is needed for internal use only as explained below + \item \textbf{form} stores a general dolfin::Form and can be used either for a + dolfn::BilinearForm as well as for a dolfin::LinearForm + \item \textbf{function} for the dolfin::Function objects + \item \textbf{functionspace} stores the user defined FunctionSpace + \item \textbf{mesh} converts a PDE-tool like mesh structure in a dolfin::Mesh +\end{itemize} +The general layout of a class is quite simple and its main purpose is the storage +of the corresponding FEniCS objects. They are managed throughout +\verb$boost::shared_ptr< >$ to the corresponding FEniCS type. +All the classes also 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. +All these classes derive from \verb$octave_base_value$. For example, the form class implementation + follows, while classes which differs from the general layout are presented below. + +\begin{lstlisting} +#ifndef _FORM_OCTAVE_ +#define _FORM_OCTAVE_ + +#include <memory> +#include <vector> +#include <dolfin.h> +#include <octave/oct.h> + +class form : public octave_base_value +{ + + public: + + form () : octave_base_value () {} + + form (const dolfin::Form _frm) + : octave_base_value (), frm (new dolfin::Form (_frm)) {} + + form (boost::shared_ptr <const dolfin::Form> _frm) + : octave_base_value (), frm (_frm) {} + + void + print (std::ostream& os, bool pr_as_read_syntax = false) const + { + os << "Form " << ": is a form of rank " << frm->rank () + << " with " << frm->num_coefficients () + << " coefficients" << std::endl; + } + + ~form(void) {} + + bool is_defined (void) const { return true; } + + const dolfin::Form & get_form (void) const { return (*frm); } + + const boost::shared_ptr <const dolfin::Form> & + get_pform (void) const { return frm; } + + private: + + boost::shared_ptr <const dolfin::Form> frm; + + DECLARE_OCTAVE_ALLOCATOR; + DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA; + +}; + +static bool form_type_loaded = false; + +DEFINE_OCTAVE_ALLOCATOR (form); +DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (form, "form", "form"); +#endif + +\end{lstlisting} + +\subsection{mesh} +In addition to usual methods, the mesh class implemens functionalities which allow to deal with meshes +as they are currently available with the msh pkg, i.e. in the (p, e, t) format. +Dolfin instead use a xml format for storing informations relative to a mesh. +We have therefore added a constructor +\begin{verbatim} + mesh (Array<double>& p, Array<octave_idx_type>& e, + Array<octave_idx_type>& t); +\end{verbatim} +which accept as input a mesh in (p, e, t) format and convert it into a xml one. Before exploring +the code in more details,the main differences between the two storing formats are presented using +the very simple but rather instructive example of a unit square mesh with just two elements, fig. \ref{mesh}. + + +\paragraph{pet} +\begin{figure} + \begin{center} + \includegraphics[height=5 cm,keepaspectratio=true]{./mesh_1.png} + \caption{The (very) simple mesh for our example} + \label{mesh} + \end{center} +\end{figure} +A mesh is represented using the three matrices p, e, t, and, using msh, +we can easily obtain the mesh for our example typing +\begin{verbatim} + mesh = msh2m_structured_mesh ([0 1], [0 1], 1, [11 12 12 13]) +\end{verbatim} + +The matrix p stores information about the coordinates of the vertices +\begin{verbatim} + >> mesh.p +\end{verbatim} + +$ +\begin{array}{rrrrl} + \; 0 & 0 & 1 & 1 & \quad \text{x-coordinates} \\ + \; 0 & 1 & 0 & 1 & \quad \text{y-coordinates} \\ +\end{array} +$ +\newline +\newline +Thus the vertex in the $n^{th}$ column will be labelled as the vertex number $n$, and so on.. + +The matrix t stores information about the connectivity +\begin{verbatim} + >> mesh.t +\end{verbatim} + +$ +\begin{array}{rrl} + \; 1 & 1 & \quad \text{number of the first vertex of the element} \\ + \; 3 & 4 & \quad \text{number of the second vertex of the element} \\ + \; 4 & 2 & \quad \text{number of the third vertex of the element} \\ + \; 0 & 0 & \\ +\end{array} +$ +\newline +\newline +The first element is thus the one obtained connecting vertices 1-3-4 and so on.. + +The matrix e stores information related to every side edge, like the number of the vertices of the boundary elements, and the number of the geometrical border containing the edge. (This last information could be very useful later in this project, when we will have to find a way to apply boundary conditions to our problem) +\begin{verbatim} + >> mesh.e +\end{verbatim} + +$ +\begin{array}{rrrrl} + \; 1 & 3 & 2 & 1 & \quad \text{first vertex of the side edge} \\ + \; 3 & 4 & 4 & 2 & \quad \text{second vertex of the side edge} \\ + \; 0 & 0 & 0 & 0 & \text{} \\ + \; 0 & 0 & 0 & 0 & \text{} \\ + \; 11 & 12 & 12 & 13 & \quad \text{label of the geometrical border containing the edge} \\ + \; 0 & 0 & 0 & 0 & \text{}\\ + \; 1 & 1 & 1 & 1 & \text{} \\ +\end{array} +$ +\newline +\newline +The side edge between vertex 1-3 is labelled 11, between 3-4 12... + +\paragraph{dolfin xml} A mesh is an object of the dolfin::Mesh class which stores information only about +the coordinates of the vertices (like p) and the information about the connectivity (like t). +A mesh can thus be manipulated using the functions and the methods of the class. +The information about boundaries is not directly stored in the mesh. +The mesh used in th example is stored as + +\begin{verbatim} +<?xml version="1.0"?> +<dolfin xmlns:dolfin="http://fenicsproject.org"> + <mesh celltype="triangle" dim="2"> + <vertices size="4"> + <vertex index="0" x="0.000e+00" y="0.000e+00" /> + <vertex index="1" x="0.000e+00" y="1.000e+00" /> + <vertex index="2" x="1.000e+00" y="0.000e+00" /> + <vertex index="3" x="1.000e+00" y="1.000e+00" /> + </vertices> + <cells size="2"> + <triangle index="0" v0="0" v1="2" v2="3" /> + <triangle index="1" v0="0" v1="1" v2="3" /> + </cells> + </mesh> +</dolfin> +\end{verbatim} + +\subsection{Shared pointer} 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. - - + +\iffalse \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: @@ -607,8 +755,6 @@ \subsection{Sparse Matrices} -\subsection{Shared pointer} - \subsection{Polymorphism} \subsection{Code release} @@ -631,11 +777,8 @@ \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: + +To reach this goal, 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.