This document is under active development and discussion!

If you find errors or omissions in this document, please don’t hesitate to submit an issue or open a pull request with a fix. We also encourage you to ask questions and discuss any aspects of the project on the Feel++ Gitter forum. New contributors are always welcome!

Preface

1. Discussion Forum

We’re always happy to help out with Feel++ or any other questions you might have. You can ask a question or signal an issue at the Gitter Feel++ salon.

Join the Feel++ chat

Join%20Chat

This book is available on Github. We use Gitter to discuss the changes in the book.

Join the Feel++ book chat

Join%20Chat

2. Conventions used in this book

The following typographical conventions are used in the book

Italic indicates new terms

typewriter is used on program listings as well as when referring to programming elements, e.g. functions, variables, statements, data types, environment variables or keywords.

$ typewriter or > typewriter displays commands that the user types literally without the $ or >.

this is a general note.
this is a general warning.
be cautious

3. Mathematical Notations

3.1. Geometry and Meshes

  • \(d=1,2,3\) geometrical dimension

  • \(\Omega \subset \mathbb{R}^d\)

  • \(K\) a cell or element of a mesh

  • \(h\) characteristic mesh size

  • \(k_{\mathrm{geo}}\) polynomial order of the geometrical transformation

  • \(\delta=(h,k_{\mathrm{geo}})\) discretization parameter pair for the geometrical transformation, default value \(k_{\mathrm{geo}}=1\) (straight cells or elements)

  • \(\varphi^K_\delta: \hat{K} \rightarrow K\), geometrical transformation

  • \(\mathcal{T}_{\delta}\) a triangulation, \(\mathcal{T}_\delta = \{ K\; | \; K=\varphi^K_\delta (\hat{K}) \} \)

  • \(\Omega_h \equiv \cup_K {K}\)

3.2. Spaces

  • \(P^k_{c,h} = \{ v_h \in C^0(\bar{\Omega}); \forall K \in \mathcal{T}_h,\ v_h \circ T_K \in \mathbb{P}^k\}\) Space of continuous piecewise polynomial of total degree \(\leq k\).

Introduction to Feel++

Discuss and Contribute
Use Issue 858 to drive development of this section. Your contributions make a difference. No contribution is too small.

4. What is Feel++?

Feel++ is a unified C++ implementation of Galerkin methods (finite and spectral element methods) in 1D, 2D and 3D to solve partial differential equations.

Feel++ is

  1. a versatile mathematical kernel solving easily problems using different techniques thus allowing testing and comparing methods, e.g. cG versus dG.

  2. a small and manageable library which nevertheless encompasses a wide range of numerical methods and techniques and in particular reduced order methods such as the reduced basis method.

  3. a software that follows closely the mathematical abstractions associated with partial differential equations (PDE) and in particular the finite element mathematical framework and variational formulations.

  4. a library that offers solving strategies that scales up to thousands and even tens of thousands of cores.

  5. a library entirely in C++ allowing to create C++ complex and typically non-linear multi-physics applications currently in industry, physics and health-care.

Quick Starts

5. Installation Quick Start

Using Feel++ inside Introduction is the recommended and fastest way to use Feel++. The Docker chapter is dedicated to Docker and using Feel++ in Docker.

We strongly encourage you to follow these steps if you begin with Feel++ in particular as an end-user.

People who would like to develop with and in Feel++ should read through the remaining sections of this chapter.

6. Usage Start

Start the Docker container feelpp/feelpp-base or feelpp/feelpp-toolboxes as follows

> docker run -it -v $HOME/feel:/feel feelpp/feelpp-toolboxes
these steps are explained in the chapter on Feel++ Containers.

Then run e.g. the Quickstart Laplacian that solves the Laplacian problem in Quickstart Laplacian sequential or in Quickstart Laplacian on 4 cores in parallel.

Quickstart Laplacian sequential
> feelpp_qs_laplacian_2d --config-file Testcases/quickstart/laplacian/feelpp2d/feelpp2d.cfg

The results are stored in Docker in

/feel/qs_laplacian/feelpp2d/np_1/exports/ensightgold/qs_laplacian/

and on your computer

$HOME/feel/qs_laplacian/feelpp2d/np_1/exports/ensightgold/qs_laplacian/

The mesh and solutions can be visualized using e.g. Parariew or Visit.

ParaView (recommended)

is an open-source, multi-platform data analysis and visualization application.

Visit

is a distributed, parallel visualization and graphical analysis tool for data defined on two- and three-dimensional (2D and 3D) meshes

Quickstart Laplacian on 4 cores in parallel
> mpirun -np 4 feelpp_qs_laplacian_2d --config-file Testcases/quickstart/laplacian/feelpp2d/feelpp2d.cfg

The results are stored in a simular place as above: just replace np_1 by np_4 in the paths above. The results should look like

ufeelpp2d

meshfeelpp2d

Solution

Mesh

7. Syntax Start

Here are some excerpts from Quickstart Laplacian that solves the Laplacian problem. It shows some of the features of Feel++ and in particular the domain specific language for Galerkin methods.

First we load the mesh, define the function space define some expressions

Laplacian problem in an arbitrary geometry, loading mesh and defining spaces
    tic();
    auto mesh = loadMesh(_mesh=new Mesh<Simplex<FEELPP_DIM,1>>);
    toc("loadMesh");

    tic();
    auto Vh = Pch<2>( mesh );
    auto u = Vh->element("u");
    auto mu = expr(soption(_name="functions.mu")); // diffusion term
    auto f = expr( soption(_name="functions.f"), "f" );
    auto r_1 = expr( soption(_name="functions.a"), "a" ); // Robin left hand side expression
    auto r_2 = expr( soption(_name="functions.b"), "b" ); // Robin right hand side expression
    auto n = expr( soption(_name="functions.c"), "c" ); // Neumann expression
    auto g = expr( soption(_name="functions.g"), "g" );
    auto v = Vh->element( g, "g" );
    toc("Vh");

Second we define the linear and bilinear forms to solve the problem

Laplacian problem in an arbitrary geometry, defining forms and solving
    tic();
    auto l = form1( _test=Vh );
    l = integrate(_range=elements(mesh),
                  _expr=f*id(v));
    l+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_2*id(v));
    l+=integrate(_range=markedfaces(mesh,"Neumann"), _expr=n*id(v));
    toc("l");

    tic();
    auto a = form2( _trial=Vh, _test=Vh);
    a = integrate(_range=elements(mesh),
                  _expr=mu*gradt(u)*trans(grad(v)) );
    a+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_1*idt(u)*id(v));
    a+=on(_range=markedfaces(mesh,"Dirichlet"), _rhs=l, _element=u, _expr=g );
    //! if no markers Robin Neumann or Dirichlet are present in the mesh then
    //! impose Dirichlet boundary conditions over the entire boundary
    if ( !mesh->hasAnyMarker({"Robin", "Neumann","Dirichlet"}) )
        a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=g );
    toc("a");

    tic();
    //! solve the linear system, find u s.t. a(u,v)=l(v) for all v
    if ( !boption( "no-solve" ) )
        a.solve(_rhs=l,_solution=u);
    toc("a.solve");

More explanations are available in Learning by examples.

Installing Feel++

8. Getting Started

This section describes the available ways to to download, compile and install Feel++.

8.1. Docker

Using Feel++ inside Introduction is the recommended and fastest way to use Feel++. The Docker is dedicated to Docker and Feel++ Containers is dedicated to Feel++ in Docker.

We strongly encourage you to follow these steps if you begin with Feel++ in particular as an end-user.

People who would like to develop with and in Feel++ should read through the remaining sections of this chapter.

8.2. System requirements

8.2.1. Compilers

Feel++ uses C++14 compilers such as GCC6 and Clang. Currently it is not mandatory to have a C++14 stantard library but it will be soon.

There used to be a major compatibility issue between llvm/clang and GCC compilers since GCC5 released the ABI tag which makes it impossible to compile Feel++ using llvm/clang with GCC5 or GCC6 standard libraries for a time. Please see the following table to understand the working C++ compiler / C++ standard library combinations.

Table 1. Table C++ compilers and standard libraries combinations
Compiler Standard Library

clang (3.6, 3.7, 3.8)

libstdc++ 4.9

clang

libc++ (corresponding clang version)

clang (3.8(requires patches), 3.9)

libstdc++ 6

GCC 6

libstdc++ 6

GCC 6.2.1 seems to be problematic on debian/testing — the tests in the testsuite fail. — GCC 6.3.1 or GCC 6.2.0 don’t have any problems.

8.2.2. Required tools and libraries

Other than C++14 compilers, Feel++ requires only a few tools and libraries, namely CMake, Boost C++ libraries and an MPI implementation such as open-mpi or mpich. The table below provides information regarding the minimum and maximum version supported. A — means it has not necessarily been tested with the latest version but we do not expect any issues. Note that for MPI, an implementation with MPI-IO support would be best.

Table 2. Table required tools to compile Feel++
Name Minimum Version Maximum Version Notes

CMake

3.0

 — 

MPI

 — 

 — 

openmpi or mpich

Boost

1.55

1.63

Here is a list of libraries that we recommend to use jointly with Feel++.

Table 3. Table optional external libraries
Library Minimum Version Maximum Version Notes

HDF5

1.8.6

1.8.16

Enables high performance I/O; Enables MED Support; Be careful on Debian/sid a more recent version of HDF5 breaks MED support

PETSc

3.2

3.7

Last is best; a requirement for parallel and high performance computing

SLEPc

3.2

3.7

last is best; a requirement for eigenvalue problem; depends on PETSc

Gmsh

2.8.7

2.16

last is best; a requirement if you want to be able to read many file formats; HDF5 version in Debian/sid currently breaks MED format support.

Superlu

superlu and superlu_dist

Suitesparse

umfpack (colamd,amd)

OpenTURNS

2.0

Uncertainty quantification

Here is a list of tools that we recommend to use jointly with Feel++.

Table 4. Table of recommended tools
Tool License Notes

Computer Aided Design

Gmsh

Open Source

Mesh Generation

Gmsh

Open Source

MeshGems

Commercial

Post-Processing

Paraview

Open Source

Ensight

Commercial

Octave

Open Source

Gmsh

Open Source

Note that all these packages are available under Debian GNU/Linux and Ubuntu. Once you have installed those dependencies, you can go to Compiling.

8.2.5. Suggested tools

Here is a list of tools that we suggest to use jointly with Feel++.

Table 5. Table of suggested tools
Tool License Notes

Computer Aided Design (CAD)

Freecad

Open Source

Salome

Open Source

HDF5 version in Debian/sid currently breaks MED format support.

Modeling, Compilation and Simulation Environment

Open Modelica

Open Source

Debugging and Profiling

Google perftools

Open Source

Valgrind

Open Source

8.3. Feel++ on Linux

We now turn to the installation of the Feel++ dependencies on Linux. Feel++ is currently support on Ubuntu (16.04, 16.10) and Debian (Sid, Testing).

8.3.1. Ubuntu

Ubuntu 16.10 Yaketti Yak

Here is the suggested installation of the Feel++ dependencies on Ubuntu 16.10

$ sudo apt-get -qq update
$ sudo apt-get install automake autoconf libtool libboost-all-dev\
  bash-completion emacs24 gmsh libgmsh-dev libopenturns-dev \
  libbz2-dev libhdf5-openmpi-dev libeigen3-dev libcgal-dev \
  libopenblas-dev libcln-dev libcppunit-dev libopenmpi-dev \
  libann-dev libglpk-dev libpetsc3.7-dev libslepc3.7-dev \
  liblapack-dev libmpfr-dev paraview python-dev libhwloc-dev \
  libvtk6-dev libpcre3-dev python-h5py python-urllib3 xterm tmux \
  screen python-numpy python-vtk6 python-six python-ply wget \
  bison sudo xauth cmake flex gcc-6 g++-6 clang-3.9 \
  clang++-3.9 git ipython openmpi-bin pkg-config
Ubuntu 16.04

Here is the suggested installation of the Feel++ dependencies on Ubuntu LTS 16.04

$ sudo apt-get install autoconf automake bash-completion bison\
 clang++-3.8 clang-3.8 cmake emacs24 flex g++-6 gcc-6 git gmsh\
  ipython libann-dev libbz2-dev libcgal-dev libcln-dev \
  libcppunit-dev libeigen3-dev libglpk-dev libgmsh-dev \
  libhdf5-openmpi-dev libhwloc-dev liblapack-dev libmpfr-dev\
   libopenblas-dev libopenmpi-dev libopenturns-dev libpcre3-dev \
   libpetsc3.6.2-dev libproj-dev libslepc3.6.1-dev libtool \
   libvtk6-dev openmpi-bin paraview pkg-config python-dev \
   python-h5py python-numpy python-ply python-six \
   python-urllib3 python-vtk6 screen sudo tmux wget xauth xterm
We are unfortunately stung by the ABI change in GCC 6 when using clang. You need to recompile the Boost C++ libraries to be able to use clang, see the section in the Annexes on Compiling Boost.

8.3.2. Debian

Debian Sid and Testing

At the time of writing there is little difference between Sid and Testing, here is the recommend dependencies installation command line:

$ apt-get -y install \
    autoconf automake bash-completion bison cmake emacs24 \
    flex git gmsh ipython libann-dev libboost-all-dev \
    libbz2-dev libcgal-dev libcln-dev libcppunit-dev \
    libeigen3-dev libglpk-dev libgmsh-dev \
    libhdf5-openmpi-dev libhwloc-dev liblapack-dev \
    libmpfr-dev libopenblas-dev libopenmpi-dev \
    libopenturns-dev libpcre3-dev libtool libvtk6-dev \
    openmpi-bin paraview petsc-dev pkg-config python-dev \
    python-h5py python-numpy python-ply python-six \
    python-urllib3 python-vtk6 screen slepc-dev sudo \
    tmux wget xauth xterm zsh
Older distributions

Unfortunately the older distributions have the ABI GCC issue with clang, e.g. Debian/jessie, or they are too old to support a simple installation procedure.

8.4. Mac OS X

Feel++ is supported on Mac OSX, starting from OS X 10.9 Mavericks to OS X 10.12 Sierra using Homebrew or MacPorts.

8.4.1. First step

Xcode is required on Mac OSX to install Feel++.

The easiest way to do so is to go through the Apple Store application and to search for Xcode. Xcode will provide the programming environment, e.g clang, for the next steps.

8.4.2. Homebrew

Introduction to HomeBrew

Homebrew is a free/open source software introduced to simplify the installation of other free/open source software on MacOS X. Homebrew is distributed under the BSD 2 Clause (NetBSD) license. For more information, visit their website.

Installation

To install the latest version of Homebrew, simply visit their website and follow the instructions. Each new package Homebrew installs is built into an intermediate place called the Cellar (usually /usr/local/Cellar) and then the packages are symlinked into /usr/local (default).

Key commands

Homebrew base command is brew. Here is a list of base available commands:

  • brew doctor: Check if the system has any problem with the current installation of Homebrew;

  • brew install mypackage: This command installs the package mypackage;

  • brew install [--devel|--HEAD] mypackage: These options respectively installs either the development version or the HEAD version of the package mypackage, if such versions are specified in the Formula file;

  • brew uninstall mypackage: This command allows to uninstall the package mypackage.

Formulas

A Formula is a Ruby script format specific to Homebrew. It allows to describe the installation process of a package. Feel[]+ uses specific Formulae that you can get in the Feel+ github repository: feelpp/homebrew-feelpp.

Installation

This section is aimed at users that do not have Homebrew already installed.
In order to build Feel++ from Homebrew, you have to do the following steps:

First install Homebrew

pass:[\( /usr/bin/ruby -e "\)](curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

then check your Homebrew installation and fix warnings/errors if necessary

$ brew doctor

Install Homebrew-science tap to get the scientific software recommended or suggested for Feel++.

$ brew tap homebrew/homebrew-science

you should see something like

==> Tapping homebrew/science
Cloning into '/usr/local/Homebrew/Library/Taps/homebrew/homebrew-science'...
remote: Counting objects: 661, done.
remote: Compressing objects: 100% (656/656), done.
remote: Total 661 (delta 0), reused 65 (delta 0), pack-reused 0
Receiving objects: 100% (661/661), 591.93 KiB | 0 bytes/s, done.
Tapped 644 formulae (680 files, 1.9M)

Next you install Feel++ tap with

brew tap feelpp/homebrew-feelpp

you should read something like

==> Tapping feelpp/feelpp
Cloning into '/usr/local/Homebrew/Library/Taps/feelpp/homebrew-feelpp'...
remote: Counting objects: 5, done.
remote: Compressing objects: 100% (5/5), done.
remote: Total 5 (delta 0), reused 4 (delta 0), pack-reused 0
Unpacking objects: 100% (5/5), done.
Tapped 1 formula (30 files, 60.7K)

The final step is to either install Feel++

$ brew install feelpp

or just Feel++ dependencies if you plan to build Feel++ from sources yourself

$ brew install --only-dependencies feelpp

Note If you encounter problems, you can fix them using brew doctor. A frequent issue is to force open-mpi with brew link --overwrite open-mpi

Advanced usage

If Homebrew is already installed on your system, you might want to customize your installation for the correct dependencies to be met for Feel++.

Feel++ Dependencies

You can browse Feel++ dependencies using the following command:

$ brew deps feelpp | column

you get the list of formulas Feel++ depends on for its installation

ann                fftw                libtool                slepc
arpack                gcc                metis                suite-sparse
autoconf        glpk                mumps                sundials
automake        gmp                netcdf                superlu
boost                gmsh                open-mpi        superlu_dist
cln                hdf5                parmetis        szip
cmake                hwloc                petsc                tbb
eigen                hypre                scalapack        veclibfort
Customizing builds

If you want to customize the compilation process for a dependency (Set debug mode, Remove checking steps, Remove the link with certain libraries, etc.), you can access to the building options with the info flag. For exemple, with open-mpi:

$ brew info open-mpi

You get various information about the open-mpi formula

open-mpi: stable 2.0.1 (bottled), HEAD
High performance message passing library
https://www.open-mpi.org/
Conflicts with: lcdf-typetools, mpich
/usr/local/Cellar/open-mpi/2.0.1 (688 files, 8.6M) *
  Built from source on 2016-09-26 at 10:36:46 with: --c++11 --with-mpi-thread-multiple
From: https://github.com/Homebrew/homebrew-core/blob/master/Formula/open-mpi.rb
==> Dependencies
Required: libevent ✔
==> Requirements
Recommended: fortran ✔
Optional: java ✔
==> Options
--c++11
        Build using C++11 mode
--with-cxx-bindings
        Enable C++ MPI bindings (deprecated as of MPI-3.0)
--with-java
        Build with java support
--with-mpi-thread-multiple
        Enable MPI_THREAD_MULTIPLE
--without-fortran
        Build without fortran support
--HEAD
        Install HEAD version

Then, you then just have to pass the needed flags, when installing the dependency.

Important: boost has to be installed with mpi and c++11 support and mumps needs to be installed with the following scotch5 support.

8.4.3. MacPorts

Introduction

MacPorts is an open-source community projet which aims to design an easy-to-use system for compiling, installing and upgrading open-source software on Mac OS X operating system. It is distributed under BSD License and facilitate the access to thousands of ports (software) without installing or compiling open-source software. MacPorts provides a single software tree which includes the latest stable releases of approximately 17700 ports targeting the current Mac OS X release (10.9). If you want more information, please visit their website.

MacPorts Installation

To install the latest version of MacPorts, please go to Installing MacPorts page and follow the instructions. The simplest way is to install it with the Mac OS X Installer using the pkg file provided on their website. It is recommended that you install X11 (X Window System) which is normally used to display X11 applications.
If you have installed with the package installer (MacPorts-2.x.x.pkg) that means MacPorts will be installed in /opt/local. From now on, we will suppose that macports has been installed in /opt/local which is the default MacPorts location. Note that from now on, all tools installed by MacPorts will be installed in /opt/local/bin or /opt/local/sbin for example (that’s here you’ll find gcc4.7 or later e.g /opt/local/bin/g++-mp-4.7 once being installed).

Key commands

In your command-line, the software MacPorts is called by the command port. Here is a list of key commands for using MacPorts, if you want more informations please go to MacPorts Commands.

  • sudo port -v selfupdate: This action should be used regularly to update the local tree with the global MacPorts ports. The option -v enables verbose which generates verbose messages.

  • port info mypackage: This action is used to get information about a port. (description, license, maintainer, etc.)

  • sudo port install mypackage: This action install the port mypackage.

  • sudo port uninstall mypackage: This action uninstall the port mypackage.

  • port installed: This action displays all ports installed and their versions, variants and activation status. You can also use the -v option to also display the platform and CPU architecture(s) for which the ports were built, and any variants which were explicitly negated.

  • sudo port upgrade mypackage: This action updgrades installed ports and their dependencies when a Portfile in the repository has been updated. To avoid the upgrade of a port’s dependencies, use the option -n.

Portfile

A Portfile is a TCL script which usually contains simple keyword values and TCL expressions. Each package/port has a corresponding Portfile but it’s only a part of a port description. Feel[]+ provides some mandatory Portfiles for its compilation which are either not available in MacPorts or are buggy but Feel+ also provides some Portfiles which are already available in MacPorts such as gmsh or petsc. They usually provide either some fixes to ensure Feel++ works properly or new version not yet available in MacPorts. These Portfiles are installed in ports/macosx/macports.

Installation

To be able to install Feel++, add the following line in /opt/local/etc/macports/source.conf at the top of the file before any other sources:

file:///<path to feel top directory>/ports/macosx/macports

Once it’s done, type in a command-line:

 $ cd <your path to feel top directory>/ports/macosx/macports
 $ sudo portindex -f

You should have an output like this:

Reading port index in pass:[\(<\)]your path to feel top directorypass:[\(>\)]/ports/macosx/macports
Adding port science/feel++
Adding port science/gmsh
Adding port science/petsc

Total number of ports parsed:   3
Ports successfully parsed:      3
Ports failed:                   0
Up-to-date ports skipped:       0

Your are now able to type

$ sudo port install feel++

It might take some time (possibly an entire day) to compile all the requirements for Feel++ to compile properly. If you have several cores on your MacBook Pro, iMac or MacBook, we suggest that you configure macports to use all or some of them.

To do that uncomment the following line in the file /opt/local/etc/macports/macports.conf

buildmakejobs        0 pass:[\(\#\)] all the cores

At the end of the sudo port install feel++, you have all dependencies installed. To build all the Makefile, \cmake is automatically launched but can have some libraries may not be found but they are not mandatory for build Feel++, only the features related to the missing libraries will be missing.

Missing ports

cmake can build Makefiles even if some packages are missing (latex2html, VTK …​). It’s not necessary to install them but you can complete the installation with MacPorts, cmake will find them by itself once they have been installed.

8.5. Building Feel++

Once the steps to install on Linux or MacOS X has been followed, we explain, in this section, how to download and build Feel++ from source.

8.5.1. For the impatient

First retrieve the source

$ git clone https://github.com/feelpp/feelpp.git

Create a build directory

$ mkdir build
$ cd build

Configure Feel++

$ CXX=clang++ ../feelpp/configure -r

Compile the Feel++ library

$ make feelpp
you can speed up the make process by passing the option -j<N> where N is the number of concurrent make sub-processes. It compiles N files at a time and respect dependencies. For example -j4 compiles 4 C++ files at a time.
Be aware that Feel++ consumes memory. The Feel++ library compile with 2Go of RAM. But to be more comfortable, 4Go or more would be best. The more, the better.

Compile your first Feel++ applications

$ make quickstart

Execute your first Feel++ application in sequential

$ cd quickstart
$ ./feelpp_qs_laplacian_2d --config-file qs_laplacian_2d.cfg

Execute your first Feel++ application using 4 mpi processes

$ mpirun -np 4 feelpp_qs_laplacian_2d --config-file qs_laplacian_2d.cfg

8.5.2. Downloading sources

Using Tarballs

Feel is distributed as tarballs following each major release. The tarballs are available on the link:https://github.com/feelpp/feelpp/releases[Feel Releases] web page.

Download the latest tarball, then uncompress it with:

$ tar -xzf feelpp-X.YY.0.tar.gz
$ cd feelpp-X.YY.0

You can now move to the section Using cmake.

Using Git

Alternatively, you can download the sources of Feel++ directly from the Git repository.

$ git clone  https://github.com/feelpp/feelpp.git

You should read something like

Cloning into 'feelpp'...
remote: Counting objects: 129304, done.
remote: Compressing objects: 100% (18/18), done.
remote: Total 129304 (delta 6), reused 0 (delta 0), pack-reused 129283
Receiving objects: 100% (129304/129304), 150.52 MiB | 1.69 MiB/s, done.
Resolving deltas: 100% (94184/94184), done.
Checking out files: 100% (7237/7237), done.
$ cd feelpp

The first level directory tree is as follows

$ tree -L 1 -d | column
.                        ├── databases                ├── research
├── applications        ├── doc                        ├── testsuite
├── benchmarks                ├── feel                └── tools
├── cmake                ├── ports                14 directories
├── contrib                ├── projects
├── data                ├── quickstart

8.5.3. Configuring Feel++

For now on, we assume that clang++ has been installed in /usr/bin. Yor mileage may vary depending on your installation of course.

It is not allowed to build the library in the top source directory.

It is recommended to have a directory (e.g. FEEL) in which you have both the sources and build directories, as follows

$ ls FEEL
feelpp/ # Sources
feel.opt/ # Build directory

feelpp is the top directory where the source have been downloaded, using git or tarballs.

Using cmake

The configuration step with cmake is as follows

$ cd FEEL/feel.opt
$ cmake ../feelpp -DCMAKE_CXX_COMPILER=/usr/bin/clang++-3.6 -DCMAKE_C_COMPILER=/usr/bin/clang-3.6 -DCMAKE_BUILD_TYPE=RelWithDebInfo

CMake supports different build type that you can set with -DCMAKE_BUILD_TYPE (case insensitive) : * None * Debug : typically -g * Release : typically -O3 -DNDEBUG * MinSizeRel : typically -Os * RelWithDebInfo : typically -g -O2 -DNDEBUG

Using configure

Alternatively you can use the configure script which calls cmake. configure --help will provide the following help.

Listing Configure help
Options:
 -b, --build                         build type: Debug, Release, RelWithDebInfo
 -d, --debug                         debug mode
-rd, --relwithdebinfo                relwithdebinfo mode
 -r, --release                       release mode
     --std=c++xx                     c++ standard: c++14, c++1z (default: c++14)
     --stdlib=libxx                  c++ standard library: stdc++(GCC), c++(CLANG) (default: stdc++)
     --max-order=x                   maximum polynomial order to instantiate(default: 3)
     --cxxflags                      override cxxflags
     --cmakeflags                    add extra cmake flags
     --prefix=PATH                   define install path
 -v, --verbose                       enable verbose output
 -h, --help                          help page
     --<package>-dir=PACKAGE_PATH    define <package> install directory
     --disable-<package>             disable <package>
     --generator=GENERATOR           cmake generator

We display below a set of possible configurations:

Compile using Release build type, default c compiler and libstdc

Listing compiling using default compilers
$ ../feelpp/configure -r

Compile using Release build type, clang compiler and libstdc

Listing compiling using clang++
$ CXX=clang++ ../feelpp/configure -r

Compile using Debug build type, clang compiler and libc

Listing compiling using clang/libc in Debug mode
CXX=clang++ ../feelpp/configure -d -stdlib=c++

8.5.4. Compiling Feel++

Once cmake or configure have done their work successfully, you are ready to compile Feel++

$ make

You can speed up the compilation process, if you have a multicore processor by specifying the number of parallel jobs make will be allowed to spawn using the -j flag:

Listing build Feel++ library using 4 concurrent jobs
$ make -j4 feelpp
From now on, all commands should be typed in build directory (e.g feel.opt) or its subdirectories.

8.5.5. Running the Feel++ Testsuite

If you encounter issues with Feel++, you can run the testsuite and send the resulting report. Feel++ has more than 300 tests running daily on our servers. Most of the tests are run both in sequential and in parallel.

The testsuite is in the testsuite directory.

$ cd testsuite

The following command will compile 10 tests at a time

$ make -j10
Listing: Running the Feel++ testsuite
$ ctest -j4 -R .

It will run 4 tests at a time thanks to the option -j4.

9. Docker

Docker is the recommended way if you are beginning using Feel++.

This chapter explains step by step how to get the Feel++ Container System(FCS), how to execute a precompiled application, how to parameter and run models.

9.1. Introduction

Container based technologies are revolutionizing development, deployment and execution of softwares. Containers encapsulate a software and allow to run seamlessly on different platforms — clusters, workstations, laptops — The developer doesn’t have to worry about specific environments and users spend less time in configuring and installing the software. Containers appear to be lightweight virtual machines (VMs) — they are started in a fraction of a second — but they, in fact, have important differences.

One of the differences is the isolation process. The VMs share only the hypervisor, the OS and hardware whereas containers may share, between each other, large parts of filesystems rather than having copies. Another difference is that, unlike in VMs, processes in a container are similar to native processes and they do not incur the overhead due to the VM hypervisor. The figure below illustrates these fundamental differences. We see in particular that the applications 2 and 3 are sharing lib 2 without redundancy.

Figure 1. Figure : VMs vs Containers

Docker is a container technology providing:

  1. an engine to start and stop containers,

  2. a user friendly interface from the creation to the distribution of containers and

  3. a hub  — cloud service for container distribution — that provides publicly a huge number of containers to download and avoid duplicating work.

9.2. Installation

This section covers briefly the installation of Docker. It should be a relatively simple smooth process to install Docker.

9.2.1. Channels

Docker offers two channels: the stable and beta channels.

stable channel

is fully baked and tested software providing a reliable platform to work with. Releases are not frequent.

beta channel

offers cutting edge features and experimental versions of the Docker Engine. This is a continuation of the initial Beta program of Docker to experiment with the latest features in development. It incurs far more instabilities than the stable channel but releases are done frequently — possibly several releases per month.

In the latter we shall consider only installing and using the stable channel.

9.2.2. Installing Docker

At the time of writing this section, Docker is available on Linux, Mac and Windows.

Mac and Windows

The support for Mac and Windows as Host OS was recently released and Docker Inc provides installation processes Docker For Mac and Docker for Windows which are the recommended way of installing Docker on these platforms.

Linux

Most Linux distributions have their own packages but they tend to lag behind the stable releases of Docker which could be a serious issue considering the development speed of Docker.

To follow Docker releases, it is probably best to use the packages distributed by Docker.

Installing Binaries

The last possibility is to use Docker Binaries to install Docker. This should be used at the last resort if packages are provided neither by your distribution nor by Docker Inc.

Tested with Docker 1.12

At the time of writing this book, the Docker version we used is Docker 1.12. All commands have been tested with this version.

9.2.3. Running without sudo

On Linux, Docker is a priviledged binary, you need to prefix all your commands with sudo, e.g. on Ubuntu. You need first to belong to the docker group with the following command on Ubuntu

$ sudo usermod -aG docker

It creates the docker group if it doesn’t already exist and adds the current user to the docker group. Then you need to log out and log in again. Similar process is available on other distributions. You need also to restart the docker service

$ sudo service docker restart
From now on, we omit the sudo command when using Docker for the sake of brevity.
Adding a user to the docker group has security implications. On a shared machine, you should consider reading the Docker security page.

9.2.4. Checking Docker

We now check your installation by running docker version To make sure everything is installed correctly and working, try running the docker version command. You should see something like the following on Linux or Mac.

Listing : Output of docker version on Linux
> docker version
Client:
 Version:      1.12.1
 API version:  1.24
 Go version:   go1.6.3
 Git commit:   23cf638
 Built:        Mon, 10 Oct 2016 21:38:17 +1300
 OS/Arch:      linux/amd64

Server:
 Version:      1.12.1
 API version:  1.24
 Go version:   go1.6.3
 Git commit:   23cf638
 Built:        Mon, 10 Oct 2016 21:38:17 +1300
 OS/Arch:      linux/amd64
Listing : Output of docker version on Mac
> docker version
Client:
Version: 1.12.6
API version: 1.24
Go version: go1.6.4
Git commit: 78d1802
Built: Wed Jan 11 00:23:16 2017
OS/Arch: darwin/amd64

Server:
Version: 1.12.6
API version: 1.24
Go version: go1.6.4
Git commit: 78d1802
Built: Wed Jan 11 00:23:16 2017
OS/Arch: linux/amd64

If so, you are ready for the next step. If instead you get something like

Listing : Bad response output of docker version on Linux
> docker version
Client:
 Version:      1.12.1
 API version:  1.24
 Go version:   go1.6.3
 Git commit:   23cf638
 Built:        Mon, 10 Oct 2016 21:38:17 +1300
 OS/Arch:      linux/amd64
Cannot connect to the Docker daemon. Is the docker daemon running on this host?
Listing : Bad response output of docker version on Mac
> docker version
Client:
 Version:      1.12.6
 API version:  1.24
 Go version:   go1.6.4
 Git commit:   78d1802
 Built:        Wed Jan 11 00:23:16 2017
 OS/Arch:      darwin/amd64
Error response from daemon: Bad response from Docker engine

then it means that the Docker daemon is not running or that the client cannot access it.

To investigate the problem you can try running the daemon manually — e.g. sudo docker daemon. This should give you some informations of what might have gone wrong with your installation.

10. Feel++ Containers

Feel++ leverages the power of Docker and provides a stack of container images.

10.1. First steps

To test Docker is installed properly, try

$ docker run feelpp/feelpp-env  echo 'Hello World!'

We have called the docker run command which takes care of executing containers. We passed the argument feelpp/feelpp-env which is a Feel++ Ubuntu 16.10 container with the required programming and execution environment for Feel++.

feelpp/ in feelpp/feelpp-env provides the organization name (or namespace) of the image and feelpp-env is the image name. Note also that Docker specifies a more complete name feelpp/feelpp-env:latest including the tag name :latest. We will see later how we defined the latest tag at the Feel++ organization. See Feel++ Container System for more details.

This may take a while depending on your internet connection but eventually you should see something like

Unable to find image 'feelpp/feelpp-env:latest' locally (1)
latest: Pulling from feelpp/feelpp-env
8e21f82d32cf: Pull complete
[...]
0a8dee947f9b: Pull complete
Digest: sha256:457539dbd781594eccd4ddf26a7aefdf08a2fff9dbeb1f601a22d9e7e3761fbc
Status: Downloaded newer image for feelpp/feelpp-env:latest
Hello World!
1 The first line tells us that there is no local copy of this Feel++ image. Docker checks automatically online on the Docker Hub if an image is available.

Once the image is downloaded, Docker launches the container and executes the command we provided echo 'Hello World!' from inside the container. The result of the command is showed on the last line of the output log above.

If you run the command again, you won’t see the download part and the command will be executed very fast.

We can ask Docker to give us a shell using the following command

$ docker run -it feelpp/feelpp-env

It provides a shell prompt from inside the container which is very similar to what you obtain when login with ssh on a remote machine. The flags -i and -t tell Docker to provide an interactive session (-i) with a TTY attached (-t).

10.1.1. Feel++ Container System

The Feel++ Container System (FCS) is organized in layers and provides a set of images.

10.1.2. Naming

The naming convention of the FCS allows the user to know where they come from and where they are stored on the Docker Hub. The name of the images is built as follows

feelpp/feelp-<component>[:tag]

where

  • feelpp/ is the namespace of the image and organization name

  • feelpp-<component> the image name and Feel++ component

  • [:tag] an optional tag for the image, by default set to :latest

Feel++ images(components) are defined as layers in the FCS in the table below.

Table 6. Table of the current components of the FCS
Component Description Built From

feelpp-env

Execution and Programming environment

<OS>

feelpp-libs

Feel++ libraries and tools

feelpp-env

feelpp-base

Feel++ base applications

feelpp-libs

feelpp-toolboxes

Feel++ toolboxes

feelpp-toolboxes

| Note: feelpp-env depends on an operating system image <OS>, the recommended and default <OS> is Ubuntu 16.10. In the future, we will build upon the next Ubuntu LTS or Debian Stable releases.

10.1.3. Tags

By default, the :latest tag is assumed in the name of the images, for example when running

$ docker run -it feelpp/feelpp-base

it is in fact feelpp/feelpp-base:latest which is being launched. The following table displays how the different images depend from one another.

Image Built from

feelpp-env:latest

Ubuntu 16.10

feelpp-libs:latest

feelpp-env:latest

feelpp-base:latest

feelpp-libs:latest

feelpp-toolboxes:latest

feelpp-base:latest

10.1.4. Host OS

As we said before the default Host OS is Ubuntu 16.10. However Docker shines in continuous integration. It provides a large set of operating system to build upon and allows to check the software in various contexts. The FCS takes advantage of Docker to build feelpp-libs for several operating systems provided by feelpp-env and with different compilers any time a commit in the Feel++ repository is done.

Table 7. Table providing the list of supported Host OS
Operating system version feelpp-env Tags Compilers

Ubuntu

16.10

ubuntu-16.10, latest

GCC 6.x, Clang 3.9

Ubuntu

16.04

ubuntu-16.04

GCC 6.x, Clang 3.8

Debian

sid

debian-sid

GCC 6.x, Clang 3.9,4.0

Debian

testing

debian-testing

GCC 6.x, Clang 3.9

If you are interested in testing Feel++ in these systems, you can run these flavors.

10.1.5. Containers

feelpp-env

feelpp-env provides the Host OS and Feel++ dependencies.

$ docker run feelpp/feelpp-env
feelpp-libs

feelpp-libs builds from feelpp-env and provides:

  1. the Feel++ libraries

  2. the Feel++ mesh partitioner application

$ docker run feelpp/feelpp-libs
feelpp-base

feelpp-base builds from feelpp-libs and provides two basic applications:

  1. feelpp_qs_laplacian_*: 2D and 3D laplacian problem

  2. feelpp_qs_stokes_*: 2D stokes problem

$ docker run feelpp/feelpp-base
feelpp-toolboxes

feelpp-toolboxes builds from feelpp-base and provides

$ docker run feelpp/feelpp-toolboxes

10.2. Running Feel++ Applications

To run Feel++ applications in docker, you need first to create a directory where you will store the Feel++ simulation files. For example, type

> mkdir $HOME/feel

and then type the following docker command

> docker run -it -v $HOME/feel:/feel feelpp/feelpp-libs

The previous command will execute the latest feelpp/feelpp-libs docker image in interactive mode in a terminal (-ti) and mount $HOME/feel in the directory /feel of the docker image.

Running the command df inside the Docker container launched by the previous command

feelpp@4e7b485faf8e:~$ df

will get you this kind of output

Filesystem     1K-blocks      Used Available Use% Mounted on
none           982046716 505681144 426457452  55% /
tmpfs          132020292         0 132020292   0% /dev
tmpfs          132020292         0 132020292   0% /sys/fs/cgroup
/dev/sda2      982046716 505681144 426457452  55% /feel
shm                65536         0     65536   0% /dev/shm

You see on the last but one line the directory $HOME/feel mounted on /feel in the Docker image.

Note that mouting a host sub-directory on /feel is mandatory. If you don’t, the Feel++ applications will exit due to lack of permissions. If you prefer running inside the docker environment you can type unset FEELPP_REPOSITORY and then all results from Feel++ applications will be store in $HOME/feel. But then you will have to use `rsync or ssh to copy your results out of the docker image if needed.

Learning Feel++

11. The Laplacian

11.1. Problem statement

We are interested in this section in the conforming finite element approximation of the following problem:

Laplacian problem

Look for \(u\) such that

\[\left\{\begin{split} -\Delta u &= f \text{ in } \Omega\\ u &= g \text{ on } \partial \Omega_D\\ \frac{\partial u}{\partial n} &=h \text{ on } \partial \Omega_N\\ \frac{\partial u}{\partial n} + u &=l \text{ on } \partial \Omega_R \end{split}\right.\]
\(\partial \Omega_D\), \(\partial \Omega_N\) and \(\partial \Omega_R\) can be empty sets. In the case \(\partial \Omega_D =\partial \Omega_R = \emptyset\), then the solution is known up to a constant.

In the implementation presented later, \(\partial \Omega_D =\partial \Omega_N = \partial \Omega_R = \emptyset\), then we set Dirichlet boundary conditions all over the boundary. The problem then reads like a standard laplacian with inhomogeneous Dirichlet boundary conditions:

Laplacian Problem with inhomogeneous Dirichlet conditions

Look for \(u\) such that

Inhomogeneous Dirichlet Laplacian problem
\[-\Delta u = f\ \text{ in } \Omega,\quad u = g \text{ on } \partial \Omega\]

11.2. Variational formulation

We assume that \(f, h, l \in L^2(\Omega)\). The weak formulation of the problem then reads:

Laplacian problem variational formulation

Look for \(u \in H^1_{g,\Gamma_D}(\Omega)\) such that

Variational formulation
\[\displaystyle\int_\Omega \nabla u \cdot \nabla v +\int_{\Gamma_R} u v = \displaystyle \int_\Omega f\ v+ \int_{\Gamma_N} g\ v + \int_{\Gamma_R} l\ v,\quad \forall v \in H^1_{0,\Gamma_D}(\Omega)\]

11.3. Conforming Approximation

We now turn to the finite element approximation using Lagrange finite element. We assume \(\Omega\) to be a segment in 1D, a polygon in 2D or a polyhedron in 3D. We denote \(V_\delta \subset H^1(\Omega)\) an approximation space such that \(V_{g,\delta} \equiv P^k_{c,\delta}\cap H^1_{g,\Gamma_D}(\Omega)\).

The weak formulation reads:

Laplacian problem weak formulation

Look for \(u_\delta \in V_\delta \) such that

\[\displaystyle\int_{\Omega_\delta} \nabla u_{\delta} \cdot \nabla v_\delta +\int_{\Gamma_{R,\delta}} u_\delta\ v_\delta = \displaystyle \int_{\Omega_\delta} f\ v_\delta+ \int_{\Gamma_{N,\delta}} g\ v_\delta + \int_{\Gamma_{R,\delta}} l\ v_\delta,\quad \forall v_\delta \in V_{0,\delta}\]
from now on, we omit \(\delta\) to lighten the notations. Be careful that it appears both the geometrical and approximation level.

11.4. Feel++ Implementation

In Feel++, \(V_{g,\delta}\) is not built but rather \(P^k_{c,\delta}\).

The Dirichlet boundary conditions can be treated using different techniques and we use from now on the elimination technique.

We start with the mesh

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<FEELPP_DIM,1>>);
the keyword auto enables type inference, for more details see Wikipedia C++11 page.

Next the discretization setting by first defining Vh=Pch<k>(mesh) \(\equiv P^k_{c,h}\), then elements of Vh and expressions f, n and g given by command line options or configuration file.

    auto Vh = Pch<2>( mesh );
    auto u = Vh->element("u");
    auto mu = doption(_name="mu");
    auto f = expr( soption(_name="functions.f"), "f" );
    auto r_1 = expr( soption(_name="functions.a"), "a" ); // Robin left hand side expression
    auto r_2 = expr( soption(_name="functions.b"), "b" ); // Robin right hand side expression
    auto n = expr( soption(_name="functions.c"), "c" ); // Neumann expression
    auto g = expr( soption(_name="functions.g"), "g" );
    auto v = Vh->element( g, "g" );

at the following line

    auto v = Vh->element( g, "g" );

v is set to the expression g, which means more precisely that v is the interpolant of g in Vh.

the variational formulation is implemented below, we define the bilinear form a and linear form l and we set strongly the Dirichlet boundary conditions with the keyword on using elimination. If we don’t find Dirichlet, Neumann or Robin in the list of physical markers in the mesh data structure then we impose Dirichlet boundary conditions all over the boundary.

    auto l = form1( _test=Vh );
    l = integrate(_range=elements(mesh),
                  _expr=f*id(v));
    l+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_2*id(v));
    l+=integrate(_range=markedfaces(mesh,"Neumann"), _expr=n*id(v));
    toc("l");

    tic();
    auto a = form2( _trial=Vh, _test=Vh);
    a = integrate(_range=elements(mesh),
                  _expr=mu*gradt(u)*trans(grad(v)) );
    a+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_1*idt(u)*id(v));
    a+=on(_range=markedfaces(mesh,"Dirichlet"), _rhs=l, _element=u, _expr=g );
    //! if no markers Robin Neumann or Dirichlet are present in the mesh then
    //! impose Dirichlet boundary conditions over the entire boundary
    if ( !mesh->hasAnyMarker({"Robin", "Neumann","Dirichlet"}) )
        a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=g );

We have the following correspondance:

Element sets Domain

elements(mesh)

\(\Omega\)

boundaryfaces(mesh)

\(\partial \Omega\)

markedfaces(mesh,"Dirichlet")

\(\Gamma_D\)

markedfaces(mesh,"Neumann")

\(\Gamma_R\)

markedfaces(mesh,"Robin")

\(\Gamma_R\)

next we solve the algebraic problem

Listing: solve algebraic system
    //! solve the linear system, find u s.t. a(u,v)=l(v) for all v
    if ( !boption( "no-solve" ) )
        a.solve(_rhs=l,_solution=u);

next we compute the \(L^2\) norm of \(u_\delta-g\), it could serve as an \(L^2\) error if \(g\) was manufactured to be the exact solution of the Laplacian problem.

    cout << "||u_h-g||_L2=" << normL2(_range=elements(mesh), _expr=idv(u)-g) << std::endl;

and finally we export the results, by default it is in the ensight gold format and the files can be read with Paraview and Ensight. We save both \(u\) and \(g\).

Listing: export Laplacian results
    auto e = exporter( _mesh=mesh );
    e->addRegions();
    e->add( "u", u );
    e->add( "g", v );
    e->save();

11.5. Testcases

The Feel++ Implementation comes with testcases in 2D and 3D.

11.5.1. circle

circle is a 2D testcase where \(\Omega\) is a disk whose boundary has been split such that \(\partial \Omega=\partial \Omega_D \cup \partial \Omega_N \cup \partial \Omega_R\).

Here are some results we can observe after use the following command

cd Testcases/quickstart/circle
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_2d --config-file circle.cfg

This give us some data such as solution of our problem or the mesh used in the application.

ucircle

meshCircle

Solution \(u_\delta\)

Mesh

11.5.2. feelpp2d and feelpp3d

This testcase solves the Laplacian problem in \(\Omega\) an quadrangle or hexadra containing the letters of Feel++

feelpp2d

After running the following command

cd Testcases/quickstart/feelpp2d
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_2d --config-file feelpp2d.cfg

we obtain the result \(u_\delta\) and also the mesh

ufeelpp2d

/images/Laplacian/TestCases/Feelpp2d/meshfeelpp2d.png[]

Solution \(u_\delta\)

Mesh

feelpp3d

We can launch this application with the current line

cd Testcases/quickstart/feelpp3d
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_3d --config-file feelpp3d.cfg

When it’s finish, we can extract some informations

ufeelpp3d

meshfeelpp3d

Solution \(u_\delta\)

Mesh

Programming Feel++

12. Step by step into Feel++ Programming

Christophe Prud’homme <@prudhomm> v1.0, 2017/03/07

12.1. Output Directories

merge with 02-SettingUpEnvironment.adoc and rewrite pending!

Feel++ generates various files that are spread over various directories. For this tutorial, it would be beneficial to check the content of these files to familiarize yourself with Feel++.

12.1.1. Environment variables

Some of Feel++ behavior can be driven by environment variables such as

  • FEELPP_REPOSITORY

  • FEELPP_WORKDIR

Both variables should point to the same place. They define the root directory where the simulation results will be stored. By default they are set to $HOME/feel.

If you want to change the root directory where the results are stored, define e.g. FEELPP_REPOSITORY. For example in the docker image feel/apps:latest, it points to /feel.

For running,

docker run -it -v $HOME/feel:/feel feelpp/apps:latest

should get you this output

# Feel++ applications

This image provides the Feel++ applications in
 . advection
 . fluid mechanics
 . solid mechanics
 . fluid-structure interaction
 . thermodynamics
The testcases are in $HOME/Testcases/ and the results are in
$HOME/feel

Follow the steps described in

http://book.feelpp.org/content/Applications/Models/readme.html

to use the models.

feelpp@50381de2bd23:~$

and here is the result of echo $FEELPP_REPOSITORY in the docker image

feelpp@50381de2bd23:~pass:[\( echo \)]FEELPP_REPOSITORY
/feel

12.1.2. Results

  • Results $FEELPP_REPOSITORY/feel/<your_app_name>/np_1

12.1.3. Log files

Feel++ uses Google Glog.

  • Log files : $FEELPP_REPOSITORY/feel/<your_app_name>/np_1/logs

  • Mesh : $FEELPP_REPOSITORY/feel/<your_app_name>/np_1

  • Config files : $FEELPP_REPOSITORY/<your_build_folder>/doc/manual/tutorial

12.2. Setting up the Feel++ Environment

merge with 01-OutputDirectories.adoc and rewrite pending!

12.2.1. Minimal Example

Let’s begin with our first program using the Feel[]+ framework. To start, you include the Feel+ headers.

We use the C[]+ `namespace` to avoid `Feel::` prefix before Feel+ objects.

We initialize the environment variables through the Feel++ Environment class, which can be found here.

#include <feel/feel.hpp>

int main( int argc, char* argv[] )
{
    using namespace Feel;

    Environment env( _argc=argc, _argv=argv,
                     _about=about( _name="env",
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org") );
    std::cout << "proc " << Environment::rank()
              <<" of "<< Environment::numberOfProcessors()
              << std::endl;

}

and the config file

myapp-solver-type=cg
# myapp-pc-type=ilu

12.2.2. Adding options

We pass command line options using the Boost Program Options, library using the prefix po:: which is a Feel[]+ alias for the Boost::program_options namespace. To add a new Feel+ option, we must create a new Feel[]+ `options_description`. You must add the default Feel+ options and the new one that we choose here as a double value. Note that the default value will be assigned if not specified by the user.

12.2.3. Compilation execution and logs

To compile a tutorial, just use the GNU make command.

make feelpp_tut_<appname>

where <appname> is the name of the application you wish to compile (here, myapp). Go to the execution directory as specified in the program, and execute it.You can list the log files created :

ls /tmp/<your login>/feelpp/feelpp_tut_myapp/

If you open one of these log, you should be able to see your value and the processor number used to compute. You can run your application on several processors using MPI :

mpirun -np 2 feelpp_tut_myapp

Note that there will be one log for each processor in that case.

12.2.4. Config files

A config file can be parsed to the program to profile your options. The default config paths are,

  • current dir

  • $HOME/Feelpp/config/

  • $INSTALL_PREFIX/share/Feelpp/config/

then you have to write inside one of these folders a file called <app_name>.cfg or feelpp_<app_name>.cfg. For example, our myapp.cfg would look like :

value=0.53

Note that you can specify the config file through the option --config-file=<path>

It’s also possible to give several configuration files with the option --config-files <path1> <path2> <path3>

./feelpp_tut_myapp --config-files ex1.cfg ex2.cfg ex3.cfg

In the case where some options are duplicated in the files, the priority is given at the end :

  • ex3.cfg can overwrite options in ex2.cfg and ex1.cfg

  • ex2.cfg can overwrite options in ex1.cfg

All files in --config-files can overwrite options given by --config-file. And all options in the command line can overwrite all options given in cfg files.

12.2.5. Initializing PETSc, SLEPc and other third party libraries

PETSc is a suite of data structures and routines for the scalable (parallel) solution of scientific applications modeled by partial differential equations. It employs the MPI standard for parallelism.

Feel++ supports the PETSc framework, the Environment takes care of initializing the associated PETSc environment.

12.3. Loading a Mesh

The next step is to load a mesh.

The loadMesh function has a _name option set by default as the default value of the --gmsh.filename option that point either to a .geo, either to a .msh, or a .h5 file. Meshes in general are more detail into this section.

auto mesh=loadMesh( _mesh=new Mesh<Simplex<2>> );

12.3.1. Exporting the Mesh for visualisation

See this section for more details about exporting and visualizing meshes.

12.3.2. Implementation

An implementation reads as follows:

#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelvf/vf.hpp>

int main( int argc, char** argv )
{
    using namespace Feel;
    using namespace Feel::vf;
    // initialize Feel++ Environment
    Environment env( _argc=argc, _argv=argv,
                     _about=about( _name="mymesh" ,
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org" ) );
    // tag::mesh[]
    // create a mesh with GMSH using Feel++ geometry tool
    auto mesh = loadMesh(_mesh=new  Mesh<Simplex<2>>);
    // end::mesh[]

    LOG(INFO) << "mesh " << soption(_name="gmsh.filename") << " loaded";

    LOG(INFO) << "volume =" << integrate( _range=elements( mesh ), _expr=cst( 1. ) ).evaluate();
    LOG(INFO) << "surface = " << integrate( _range=boundaryfaces( mesh ), _expr=cst( 1. ) ).evaluate();
}

and the associated config file

[gmsh]
hsize=1e-1

12.4. Defining and using expressions

The next step is to construct a function space over the mesh.

12.4.1. Step by step explanations

  • We start by loading a Mesh in 2D.

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
  • then we define some expression through the command line of config file: g is a scalar field and f is a vector field, here is an example how to enter them :

./feelpp_tut_myexpression --a=3 --functions.g="a*x*y:x:y:a" --functions.f="{sin(pi*x),cos(pi*y)}:x:y"

You can print back the expression to the screen to check that everything is ok. You want to use as expression a*x+b*y, you have to define a and b as option (either in your code, either in the library).

  • then we compute the gradient of g and f.

    auto grad_g=grad<2>(g);
    auto grad_f=grad(f);
    std::cout << "grad(g)=" << grad_g << std::endl;
    std::cout << "grad(f)=" << grad_f << std::endl;
template argument are given to grad to specify the shape of the gradient: in the case of \(\nabla g\), it is \(1\times2\) and \(2\times 2\) for \(\nabla f\) since we are in 2D.
  • then we compute the laplacian of g and f.

    auto laplacian_g=laplacian(g);
    std::cout << "laplacian(g)=" << laplacian_g << std::endl;

    auto laplacian_f=laplacian(f);
    std::cout << "laplacian(f)=" << laplacian_f << std::endl;
  • then we compute the divergence of f.

    auto div_f=div(f);
    std::cout << "div(f)=" << div_f << std::endl;
  • and the curl of f

    auto curl_f=curl(f);
    std::cout << "curl(f)=" << curl_f << std::endl;
  • Finally we evaluate these expressions at one point given by the option x and y.

12.4.2. Implementation

An implementation reads as follows:

#include <feel/feelcore/environment.hpp>
#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feelvf/ginac.hpp>
using namespace Feel;

inline
po::options_description
makeOptions()
{
    po::options_description EXPRoptions( "DAR options" );
    EXPRoptions.add_options()
    ( "a", po::value<double>()->default_value( 1 ), "a parameter" )
    ( "b", po::value<double>()->default_value( 2 ), "a parameter" )
    ;
    return EXPRoptions;
}


int main(int argc, char**argv )
{
    Environment env( _argc=argc, _argv=argv,
                     _desc=makeOptions(),
                     _about=about(_name="myexpression",
                                  _author="Feel++ Consortium",
                                  _email="feelpp-devel@feelpp.org"));

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);

    auto g = expr(soption(_name="functions.g"));
    std::cout << "g=" << g << std::endl;

    auto f = expr<2,1>(soption(_name="functions.f"));
    std::cout << "f=" << f << std::endl;

    double aVal = doption("a")+doption("b");
    std::map<std::string,double> myMap; myMap["aVal"]=aVal;
    auto i = expr(soption("functions.i"),myMap);
    std::cout << "i=" << i << std::endl;

    auto grad_g=grad<2>(g);
    auto grad_f=grad(f);
    std::cout << "grad(g)=" << grad_g << std::endl;
    std::cout << "grad(f)=" << grad_f << std::endl;

    auto laplacian_g=laplacian(g);
    std::cout << "laplacian(g)=" << laplacian_g << std::endl;

    auto laplacian_f=laplacian(f);
    std::cout << "laplacian(f)=" << laplacian_f << std::endl;


    auto div_f=div(f);
    std::cout << "div(f)=" << div_f << std::endl;

    auto curl_f=curl(f);
    std::cout << "curl(f)=" << curl_f << std::endl;

    std::cout << "Evaluation  at  (" << doption("x") << "," << doption("y") << "):" << std::endl;
    std::cout << "           g(x,y)=" << g.evaluate() << std::endl;
    std::cout << "           f(x,y)=" << f.evaluate() << std::endl;
    std::cout << "           i(x,y)=" << i.evaluate() << std::endl;
    std::cout << "Gradient:\n";
    std::cout << "     grad(g)(x,y)=" << grad_g.evaluate() << std::endl;
    std::cout << "     grad(f)(x,y)=" << grad_f.evaluate() << std::endl;
    std::cout << "Divergence:\n";
    std::cout << "      div(f)(x,y)=" << div_f.evaluate() << std::endl;
    std::cout << "Curl:\n";
    std::cout << "     curl(f)(x,y)=" << curl_f.evaluate() << std::endl;
    std::cout << "Laplacian:\n";
    std::cout << "laplacian(g)(x,y)=" << laplacian_g.evaluate() << std::endl;
    std::cout << "laplacian(f)(x,y)=" << laplacian_f.evaluate() << std::endl;
}

and the associated config file

a=12
b=-1
[functions]
g=(a-x)*x+(a/b)*y^3:x:y:a:b
f={1,1}:x:y
i=(x-aVal)*y:x:y:aVal

12.4.3. Execution

$ ./feelpp_tut_myexpression

or

$ ./feelpp_tut_myexpression --a=3 --functions.g="<your_function>" --functions.f="<your_function>"

We start with the following function \(\)g=1\(\) and \(\)f=(1,1)\(\).

 $./feelpp_tut_myexpression --functions.g=1:x:y --functions.f="{1,1}:x:y"

and get something like this

g=1
f={1,1}
i=(x-aVal)*y
grad(g)=[[0,0]]
grad(f)=[[0,0],[0,0]]
laplacian(g)=[[0]]
laplacian(f)=[[0],[0]]
div(f)=[[0]]
curl(f)=[[0]]
Evaluation  at  (0,0):
           g(x,y)=1
           f(x,y)=1
1
           i(x,y)=-0
Gradient:
     grad(g)(x,y)=0 0
     grad(f)(x,y)=0 0
0 0
Divergence:
      div(f)(x,y)=0
Curl:
     curl(f)(x,y)=0
Laplacian:
laplacian(g)(x,y)=0
laplacian(f)(x,y)=0
0

The symbolic calculus system worked as expected.

12.5. Evaluating function

Once you have created an element, you may want to give it a value, that can depends on a lot of parameters ( mainly spaces, but others may apply ).

To do so, Feel++ relies on expressions. We may use various kind of expressions :

12.5.1. Built-in

Let’s begin with the evaluation of the expression \(sin(\pi x)\) on a unit circle.

First at all, we define the unit circle and its function space :

  // circle - geometrical order: 2
  auto mesh = unitCircle<2>();

  // circle - geometrical order: 1
  auto meshp1 = unitCircle<1>();

  auto Xh = Pch<2>( mesh );

Then the expression we would like to evaluate :

  auto myExpr = sin(pi*Px());

Px() refers to the variable \(\)x\(\) of our space.

With this,we can project it on our function space :

  auto v = project( _space=Xh, _range=elements(mesh),
      _expr=myExpr);

The expression will be evaluated on each point of our mesh.

In order to visualize the result, we create an exporter, named exhi, and add to it the projection.

  auto v = project( _space=Xh, _range=elements(mesh),
      _expr=myExpr);
Code
#include <feel/feel.hpp>

int main(int argc, char**argv )
{
  using namespace Feel;
  Environment env( _argc=argc, _argv=argv,
      _about=about(_name="myexporter",
        _author="Christophe Prud'homme",
        _email="christophe.prudhomme@feelpp.org"));

  // circle - geometrical order: 2
  auto mesh = unitCircle<2>();

  // circle - geometrical order: 1
  auto meshp1 = unitCircle<1>();

  auto Xh = Pch<2>( mesh );

  auto myExpr = sin(pi*Px());

  auto v = project( _space=Xh, _range=elements(mesh),
      _expr=myExpr);

  auto exhi = exporter( _mesh=mesh, _name="exhi" );
  auto exlo = exporter( _mesh=meshp1, _name="exlo" );
  auto exhilo = exporter( _mesh=lagrangeP1(_space=Xh)->mesh(),_name="exhilo");

  int max = 10; double dt = 0.1;
  double time = 0;
  for (int i = 0; i<max; i++)
  {
    exhilo->step( time )->add( "vhilo", v );
    exlo->step( time )->add( "vlo", v );
    exhi->step( time )->add( "vhi", v );
    time += dt;
    exhi->save();
    exlo->save();
    exhilo->save();
  }
}

The list of the Feel++ Keyword is here.

12.5.2. Hard Coded

In this second method, we will use Functor :

struct Functor
{
    static const size_type context = vm::JACOBIAN|vm::POINT;
    typedef double value_type;
    typedef Feel::uint16_type uint16_type;
    static const uint16_type rank = 0;
    static const uint16_type imorder = 1;
    static const bool imIsPoly = true;
    myFunctor( int i ): coord(i){}
    double operator()( uint16_type a, uint16_type b, ublas::vector<double> const& x, ublas::vector<double> const& n ) const
        {
            return xtag::coord[];
        }
    int coord = 0;
};

We create a unit square meshed by triangles and we define the associated function space :

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    auto Xh= Pch<2>(mesh);

From this space, we can define two elements, here one equals to the variable \(x\) and the other to the variable \(y\), obtain from Functor class.

    auto u = Xh.element(idf(Functor(0))); // Will contain x
    auto v = Xh.element(idf(Functor(1))); // Will contain y

The data exportation is the final step to visualize our expression \(\)x\(\) and \(\)y\(\) on the defined mesh.

[import:"exporter",unindent:"true"](../codes/myfunctor.cpp)

Code
#include <feel/feelcore/environment.hpp>
#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feeldiscr/pch.hpp>
#include <feel/feeldiscr/pchv.hpp>
#include <feel/feelvf/projectors.hpp>
#include <feel/feeldiscr/operatorlagrangep1.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelvf/function.hpp>
using namespace Feel;

namespace Feel {
struct Functor
{
    static const size_type context = vm::JACOBIAN|vm::POINT;
    typedef double value_type;
    typedef Feel::uint16_type uint16_type;
    static const uint16_type rank = 0;
    static const uint16_type imorder = 1;
    static const bool imIsPoly = true;
    myFunctor( int i ): coord(i){}
    double operator()( uint16_type a, uint16_type b, ublas::vector<double> const& x, ublas::vector<double> const& n ) const
        {
            return xtag::coord[];
        }
    int coord = 0;
};
}

int main(int argc, char**argv )
{
    Environment env( _argc=argc, _argv=argv,
                     _about=about(_name="myfunctor",
                                  _author="Feel++ Consortium",
                                  _email="feelpp-devel@feelpp.org"));


    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    auto Xh= Pch<2>(mesh);


    auto u = Xh.element(idf(Functor(0))); // Will contain x
    auto v = Xh.element(idf(Functor(1))); // Will contain y

    auto ex1 = exporter(_mesh=mesh,_name="ex1");
    ex1->add("u",u);
    ex1->add("v",v);
    ex1->save();
}

The list of the Feel++ Keyword is here.

12.6. Visualizing functions over a mesh

The next step is to visualize function over the mesh. The source code that generate our output is available in myexporter.cpp, presented in the previous section.

You can visualize data via :

The results files are in

  • $HOME/feel/myexporter/np_1

or

  • $FEELPP_WORKDIR/feel/myexporter/np_1.

We discriminate output directories based on the name of the simulation (the parameter _name in the environment), the number of process (mpirun -np …​) and the type of the chosen exporter

--exporter.format={ensight|ensightgold|gmsh|...}

12.7. Spaces and elements

You’ve learned how to discretize the space you want to compute on. You now have to learn how to define and use function spaces and elements of functions spaces. For advanced informations on this subject, you can look there.

12.7.1. Constructing a function space

  • Loading a Mesh in 2D

auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
  • For basic function spaces, we have predetermined constructors:

auto Xh = Pch<2>( mesh );
  • Defining an element

auto u = Xh->element( "u" );
auto w = Xh->element( "w" );

One can also use :

  • Pdh<ORDER>(mesh) : Polynomial Discontinuous

  • Pvh<ORDER>(mesh) : Polynomial Continuous Vectorial

  • Pdhv<ORDER>(mesh) : Polynomial Discontinuous Vectorial

  • Pchm<ORDER>(mesh) : Polynomial Continuous Matrix

  • Ned1h<ORDER>(mesh) : Nedelec function spaces

12.7.2. Implementation

The implementation reads are follows

#include <feel/feel.hpp>
using namespace Feel;

int main( int argc, char** argv )
{

    //Initialize Feel++ Environment
    Environment env( _argc=argc, _argv=argv,
                     _about=about( _name="myfunctionspace",
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org" )  );

    // create the mesh
    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);

    // function space \( X_h \) using order 2 Lagrange basis functions
    auto Xh = Pch<2>( mesh );

    auto g = expr<4>( soption(_name="functions.g"));
    auto gradg = grad<3>(g);

    // elements of \( u,w \in X_h \)
    auto u = Xh->element( "u" );
    auto w = Xh->element( "w" );

                // build the interpolant of u
    u.on( _range=elements( mesh ), _expr=g );

                // build the interpolant of the interpolation error
    w.on( _range=elements( mesh ), _expr=idv( u )-g );

    // compute L2 norms \(||\cdot||_{L^2}\)
    double L2g = normL2( elements( mesh ), g );
    double H1g = normL2( elements( mesh ), _expr=g,_grad_expr=gradg );
    double L2uerror = normL2( elements( mesh ), ( idv( u )-g ) );
    double H1uerror = normH1( elements( mesh ), _expr=( idv( u )-g ),
                              _grad_expr=( gradv( u )-gradg ) );
    if ( Environment::isMasterRank() )
    {
        std::cout << "||u-g||_0 = " << L2uerror/L2g << std::endl;
        std::cout << "||u-g||_1 = " << H1uerror/H1g << std::endl;
    }

    // export for post-processing
    auto e = exporter( _mesh=mesh );

                 // save interpolant
    e->add( "g", u );

                 // save interpolant of interpolation error
    e->add( "u-g", w );

    e->save();

}

and the associated config file

[gmsh]
# default is hypercube
shape=hypercube
# you can try ellipsoid
#shape=ellipsoid

[functions]
g=x*y*(y*2.2+4.3)*(2.1*x-1.3):x:y

12.8. Computing integrals over mesh

The next step is to compute integrals over the mesh ( See this for detailed methods ).

12.8.1. Step by step explanations

  • We start by loading a Mesh in 2D

  • then we define the Feel[]+ expression that we are going to integrate using the soption function that retrieves the command line option string `functions.g`. We then transform this string into a Feel+ expression using expr().

  • then we compute two integrals over the domain and its boundary respectively

    • \(\int_\Omega g\)

    • \(\int_{\partial \Omega} g\)

  • and we print the results to the screen.

Only the rank 0 process (thanks to Environment) isMasterRank() prints to the screen as the result is the same over all mpi processes if the application was run in parallel. Note also that the code actually prints the expression passed by the user through the command line option functions.g.

12.8.2. Some results

We start with the following function \(g=1\). Recall that by default the domain is the unit square, hence the \(\int_\Omega g\) and \(\int_{\partial\Omega} g\) should be equal to 1 and 4 respectively.

./feelpp_tut_myintegrals --functions.g=1
int_Omega 1 = 1
int_{boundary of Omega} 1 = 4

Now we try \(g=x\). We need to tell Feel++ what are the symbols associated with the expression: here the symbol is x and it works as follows

./feelpp_tut_myintegrals --functions.g=x:x
int_Omega x = 0.5
int_{boundary of Omega} x = 2
remember that there is a separator : between the expression and each symbol composing it.

Now we try \(g=x y\). We need to tell Feel++ what are the symbols associated with the expression: here the symbol is x and y and it works as follows

./feelpp_tut_myintegrals --functions.g=x*y:x:y
int_Omega x*y = 0.25
int_{boundary of Omega} x*y = 1

More complicated functions are of course doable, such as \(g=\sin( x y ).\)

./feelpp_tut_myintegrals --functions.g="sin(x*y):x:y"
int_Omega sin(x*y) = 0.239812
int_{boundary of Omega} sin(x*y) = 0.919395

Here is the last example in parallel over 4 processors which returns, of course, the exact same results as in sequential

mpirun -np 4 ./feelpp_doc_myintegrals --functions.g="sin(x*y):x:y"
int_Omega sin(x*y) = 0.239812
int_{boundary of Omega} sin(x*y) = 0.919395

Finally we can change the type of domain and compute the area and perimeter of the unit disk as follows

./feelpp_doc_myintegrals --functions.g="1:x:y" --gmsh.domain.shape=ellipsoid --gmsh.hsize=0.05
int_Omega 1 = 0.784137
int_{boundary of Omega} 1 = 3.14033

Note that we don’t get the exact results due to the fact that [stem]:[\Omega_h = \cup_{K \in \mathcal{T}_h} K] which we use for the numerical integration is different from the exact domain \(\Omega = \{ (x,y)\in \mathbb{R}^2 | x^2+y^2 < 1\}\).

12.8.3. Implementation

To compile just type

$ ./feelpp_tut_myintegrals

The complete code reads as follows :

#include <feel/feel.hpp>


using namespace Feel;

int
main( int argc, char** argv )
{
    // Initialize Feel++ Environment
    Environment env( _argc=argc, _argv=argv,
                     _about=about( _name="myintegrals" ,
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org" ) );

    /// [mesh]
    // create the mesh (specify the dimension of geometric entity)
    auto mesh = loadMesh( _mesh=new Mesh<Simplex<2>> );
    /// [mesh]

    /// [expression]
    // our function to integrate
    auto g = expr( soption(_name="functions.g") );
    /// [expression]

    /// [integrals]
    // compute integral of g (global contribution): \(\int_{\partial \Omega} g\)
    auto intf_1 = integrate( _range = elements( mesh ),
                                 _expr = g ).evaluate();

    // compute integral g on boundary: \( \int_{\partial \Omega} g \)
    auto intf_2 = integrate( _range = boundaryfaces( mesh ),
                             _expr = g ).evaluate();

    // compute integral of grad f (global contribution): \( \int_{\Omega} \nabla g \)
    auto grad_g = grad<2>(g);
    auto intgrad_f = integrate( _range = elements( mesh ),
                                _expr = grad_g ).evaluate();

    // only the process with rank 0 prints to the screen to avoid clutter
    if ( Environment::isMasterRank() )
        std::cout << "int_Omega " << g << " = " << intf_1  << std::endl
                  << "int_{boundary of Omega} " << g << " = " << intf_2 << std::endl
                  << "int_Omega grad " << g << " = "
                  << "int_Omega  " << grad_g << " = "
                  << intgrad_f  << std::endl;
    /// [integrals]
}
/// [all]
[functions]
g=x*y*(y*2.2+4.3)*(2.1*x-1.3):x:y

12.9. Solve a partial differential equation

With all the previously notions we approach, the definition of a partial differential equation and boundary conditions are our next step. More details on these aspects can be retrieve at this page.

Variational formulation

This example refers to a laplacian problem, define by

Strong formulation
\[-\Delta u=1 \text{ in } \Omega=[0,1]^2, \quad u=1 \text{ on } \partial \Omega\]

After turning the Strong formulation into its weak form, we have

\[\int_\Omega \nabla u \cdot \nabla v = \int_\Omega v,\quad u=1 \text{ on } \partial \Omega\]

where \(u\) is the unknown and \(v\) a test function. The left side is known as the bilinear form \(a\) and the right side is the linear form \(l\). This equation can be so written :

\[\int_\Omega a(u,v) = \int_\Omega l(v), \quad u=1 \text{ on } \partial \Omega\]
Implementation

The steps to implement this problem are

  • Loading a 2D mesh, creating the function space \(V_h\), composed of piecewise polynomial functions of order 2, and its associated elements

 auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
 auto Vh = Pch<2>( mesh );
 auto u = Vh->element();
 auto v = Vh->element();
  • Define the linear form \(l\) with test function space \(V_h\)

auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh),
              _expr=id(v));
  • Define the bilinear form \(a\) with \(V_h\) as test and trial function spaces

auto a = form2( _trial=Vh, _test=Vh);
a = integrate(_range=elements(mesh),
              _expr=gradt(u)*trans(grad(v)) );

form1 and form2 are used to define respectively the left and right side of our partial differential equation.

  • Add Dirichlet boundary condition on \(\)u\(\)

a+=on(_range=boundaryfaces(mesh),
      _rhs=l, _element=u, _expr=cst(0.) );

We impose, in this case, \(u=0\) on \(\partial\Omega\), with the keyword on.

  • Solving the problem

a.solve(_rhs=l,_solution=u);
  • Exporting the solution

auto e = exporter( _mesh=mesh );
e->add( "u", u );
e->save();

The complete code reads as follows :

#include <feel/feel.hpp>

int main(int argc, char**argv )
{
     using namespace Feel;

     Environment env( _argc=argc, _argv=argv,
                      _about=about(_name="qs_laplacian",
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org"));

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    auto Vh = Pch<2>( mesh );
    auto u = Vh->element();
    auto v = Vh->element();

    auto l = form1( _test=Vh );
    l = integrate(_range=elements(mesh),
                  _expr=id(v));

    auto a = form2( _trial=Vh, _test=Vh);
    a = integrate(_range=elements(mesh),
                  _expr=gradt(u)*trans(grad(v)) );
    a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=cst(0.) );
    a.solve(_rhs=l,_solution=u);

    auto e = exporter( _mesh=mesh );
    e->add( "u", u );
    e->save();
}

and the corresponding config file

[gmsh]
hsize=1e-1

12.10. Using a backend

12.10.1. Introduction

After the discretization process, one may have to solve a (non) linear system. Feel++ interfaces with PETSc/SLEPc and Eigen3. Consider this system

\[A x = b\]

We call Backend an object that manages the solution strategy to solve it. Some explanation are available at Solver and Preconditioner.

Feel++ provides a default backend that is mostly hidden to the final user. In many examples, you do not have to take care of the backend. You change the backend behavior via the command line or config files. For example

./feelpp_doc_mybackend --backend.pc-type=id

will use the identity matrix as a right preconditionner for the default backend. The size of the preconditionner will be defined from the size of the A matrix.

If you try to solve a different system \(A_1 y= c\) (in size) with the same backend or the default without rebuilding it, it will fail.

backend(_rebuild=true)->solve(_matrix=A1,_rhs=c,_sol=y);

Each of that options can be retrieved via the --help-lib argument in the command line.

Non default Backend

You may need to manage more than one backend in an application: you have different systems to solve and you want to keep some already computed objects such as preconditioners.

  • The default backend is in fact an unnamed backend: in order to distinguish between backend you have to name them. For example

        po::options_description app_options( "MyBackend options" );
        app_options.add(backend_options("myBackend"));

        Environment env(_argc=argc, _argv=argv,
                        _desc = app_options,
                        _about=about(_name="mybackend",
                                     _author="Feel++ Consortium",
                                     _email="feelpp-devel@feelpp.org"));
  • After that, you create the backend object:

        // create a backend
        boost::shared_ptr<Backend<double>> myBackend(backend(_name="myBackend"));
the backend’s name has to match the name you gave at the options step.
  • Then, you load meshes, creates spaces etc. At solve time, or you solve with the default backend:

        // solve a(u,v)=l(v)
        if ( Environment::isMasterRank() )
                std::cout << "With default backend\n";
        a.solve(_rhs=l,_solution=u1); // Compute with default backend

One of the important backend option is to be able to monitor the residuals and iteration count

./feelpp_tut_mybackend --pc-type=id --ksp-monitor=true --myBackend.ksp-monitor=true
  • Finally you can create a named backend:

        // solve a(u,v)=l(v)
        if ( Environment::isMasterRank() )
                std::cout << "With named backend\n";
        a.solveb(_rhs=l,_solution=u2, _backend=myBackend); // Compute with myBackend

12.10.2. Implementation

The implementation reads as follows:

#include <feel/feel.hpp>
using namespace Feel;

int main(int argc, char**argv )
{
        po::options_description app_options( "MyBackend options" );
        app_options.add(backend_options("myBackend"));

        Environment env(_argc=argc, _argv=argv,
                        _desc = app_options,
                        _about=about(_name="mybackend",
                                     _author="Feel++ Consortium",
                                     _email="feelpp-devel@feelpp.org"));

        // create a backend
        boost::shared_ptr<Backend<double>> myBackend(backend(_name="myBackend"));

        // create the mesh
        auto mesh = loadMesh(_mesh=new Mesh<Simplex< 2 > > );

        // function space
        auto Vh = Pch<2>( mesh );

        // element in Vh
        auto u  = Vh->element();
        auto u1 = Vh->element();
        auto u2 = Vh->element();

        // left hand side
        auto a = form2( _trial=Vh, _test=Vh );
        a = integrate(_range=elements(mesh),
                      _expr=expr(soption("functions.alpha"))*trace(gradt(u)*trans(grad(u))) );

        // right hand side
        auto l = form1( _test=Vh );
        l = integrate(_range=elements(mesh),
                      _expr=expr(soption("functions.f"))*id(u));

        // BC
        a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
              _expr=expr(soption("functions.g")));

        // solve a(u,v)=l(v)
        if ( Environment::isMasterRank() )
                std::cout << "With default backend\n";
        a.solve(_rhs=l,_solution=u1); // Compute with default backend
        // solve a(u,v)=l(v)
        if ( Environment::isMasterRank() )
                std::cout << "With named backend\n";
        a.solveb(_rhs=l,_solution=u2, _backend=myBackend); // Compute with myBackend

        // save results
        auto e = exporter( _mesh=mesh );
        e->step(0) -> add( "u", u1 );
        e->step(1) -> add( "u", u2 );
        e->save();
}

and the associated config file

ksp-monitor=true
ksp-rtol=1e-5
backend.verbose=true
pc-type=id

[functions]
alpha=1
f=x+y:x:y
g=sin(x):x

[myBackend]
backend.verbose=true
ksp-monitor=true
ksp-rtol=1e-5
pc-type=lu

12.11. Defining a Model

12.11.1. Introduction

It is well known an equation can describe a huge range of physical problems. Each of theses problems will have a particular environment, but the equation to solve will be the same. To make our program applicable to theses range of problem, we have defined a model. Models definitions can be retrieve in this section.

12.11.2. What is a model

A model is defined by :

  • a Name

  • a Description

  • a Model

  • Parameters

  • Materials

  • Boundary Conditions

  • Post Processing

Parameters

A parameter is a non physical property for a model.

Materials

To retrieve the materials properties, we use :

ModelMaterials materials = model.materials();
BoundaryConditions

Thanks to GiNaC, we handle boundary conditions (Dirichlet, Neumann, Robin) as expression. You have to indicate in the json file the quantity to handle (velocity, pressure…​) and the associated expression.

map_scalar_field<2> bc_u { model.boundaryConditions().getScalarFields<2>("heat","dirichlet") };

We can apply theses boundary condition this way

  for(auto it : bc_u){
    if(boption("myVerbose") && Environment::isMasterRank() )
      std::cout << "[BC] - Applying " << it.second << " on " << it.first << std::endl;
    a+=on(_range=markedfaces(mesh,it.first), _rhs=l, _element=u, _expr=it.second );
  }
Code
  for(auto it : materials)
  {
    auto mat = material(it);
    if(boption("myVerbose") && Environment::isMasterRank() )
      std::cout << "[Materials] - Laoding data for " << it.second.name() << " that apply on marker " << it.first  << " with diffusion coef ["
#if MODEL_DIM == 3
        << "[" << it.second.k11() << "," << it.second.k12() << "," << it.second.k13() << "],"
        << "[" << it.second.k12() << "," << it.second.k22() << "," << it.second.k23() << "],"
        << "[" << it.second.k13() << "," << it.second.k23() << "," << it.second.k33() << "]]"
#else
        << "[" << it.second.k11() << "," << it.second.k12() << "],"
        << "[" << it.second.k12() << "," << it.second.k22() << "]]"
#endif
        << std::endl;
    k11.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k11()));
    k12.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k12()));
    k22.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k22()));
#if MODEL_DIM == 3
    k13 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k13());
    k23 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k23());
    k33 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k33());
#endif
  }
#if MODEL_DIM == 2
  a += integrate(_range=elements(mesh),_expr=inner(mat<MODEL_DIM,MODEL_DIM>(idv(k11), idv(k12), idv(k12), idv(k22) )*trans(gradt(u)),trans(grad(v))) );
#else
  a += integrate(_range=elements(mesh),_expr=inner(mat<MODEL_DIM,MODEL_DIM>(idv(k11), idv(k12), idv(k13), idv(k12), idv(k22), idv(k23), idv(k31), idv(k32), idv(k33))*trans(gradt(u)),trans(grad(v))) );
#endif
PostProcessing

TODO: explanation pending.

Example

We have set up an example : an anisotropic laplacian.

#include <feel/feel.hpp>
#include <feel/feelmodels/modelproperties.hpp>
int main(int argc, char**argv )
{
  using namespace Feel;
  po::options_description laplacianoptions( "Laplacian options" );
  laplacianoptions.add_options()
    ("myVerbose", po::value< bool >()-> default_value( true ), "Display information during execution")
    ;
  Environment env( _argc=argc, _argv=argv,
      _desc=laplacianoptions,
      _about=about(_name="aniso_laplacian",
        _author="Feel++ Consortium",
        _email="feelpp-devel@feelpp.org"));
  ModelProperties model; // Will load --mod-file
  map_scalar_field<2> bc_u { model.boundaryConditions().getScalarFields<2>("heat","dirichlet") };
  ModelMaterials materials = model.materials();
  if(boption("myVerbose") && Environment::isMasterRank() )
    std::cout << "Model " << Environment::expand( soption("mod-file")) << " loaded." << std::endl;
  auto f = expr( soption(_name="functions.f"), "f" );
  auto mesh = loadMesh(_mesh=new Mesh<Simplex<MODEL_DIM>>);
  auto Vh = Pch<2>( mesh );
  auto u = Vh->element();
  auto v = Vh->element();
  auto k11 = Vh->element();
  auto k12 = Vh->element();
  auto k22 = Vh->element();
#if MODEL_DIM == 3
  auto k13 = Vh->element();
  auto k23 = Vh->element();
  auto k33 = Vh->element();
#endif
  auto a = form2( _trial=Vh, _test=Vh);
  auto l = form1( _test=Vh );
  l = integrate(_range=elements(mesh),_expr=f*id(v));
  for(auto it : materials)
  {
    auto mat = material(it);
    if(boption("myVerbose") && Environment::isMasterRank() )
      std::cout << "[Materials] - Laoding data for " << it.second.name() << " that apply on marker " << it.first  << " with diffusion coef ["
#if MODEL_DIM == 3
        << "[" << it.second.k11() << "," << it.second.k12() << "," << it.second.k13() << "],"
        << "[" << it.second.k12() << "," << it.second.k22() << "," << it.second.k23() << "],"
        << "[" << it.second.k13() << "," << it.second.k23() << "," << it.second.k33() << "]]"
#else
        << "[" << it.second.k11() << "," << it.second.k12() << "],"
        << "[" << it.second.k12() << "," << it.second.k22() << "]]"
#endif
        << std::endl;
    k11.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k11()));
    k12.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k12()));
    k22.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k22()));
#if MODEL_DIM == 3
    k13 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k13());
    k23 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k23());
    k33 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k33());
#endif
  }
#if MODEL_DIM == 2
  a += integrate(_range=elements(mesh),_expr=inner(mat<MODEL_DIM,MODEL_DIM>(idv(k11), idv(k12), idv(k12), idv(k22) )*trans(gradt(u)),trans(grad(v))) );
#else
  a += integrate(_range=elements(mesh),_expr=inner(mat<MODEL_DIM,MODEL_DIM>(idv(k11), idv(k12), idv(k13), idv(k12), idv(k22), idv(k23), idv(k31), idv(k32), idv(k33))*trans(gradt(u)),trans(grad(v))) );
#endif
  for(auto it : bc_u){
    if(boption("myVerbose") && Environment::isMasterRank() )
      std::cout << "[BC] - Applying " << it.second << " on " << it.first << std::endl;
    a+=on(_range=markedfaces(mesh,it.first), _rhs=l, _element=u, _expr=it.second );
  }
  a.solve(_rhs=l,_solution=u);
  auto e = exporter( _mesh=mesh );
  for(int i = 0; i < 3; i ++){
    for(auto const &it : model.postProcess()["Fields"] )
    {
      if(it == "diffused")
        e->step(i)->add("diffused",u);
      else if(it == "k11")
        e->step(i)->add("k11",k11);
      else if(it == "k12")
        e->step(i)->add("k12",k12);
      else if(it == "k11")
        e->step(i)->add("k22",k22);
#if MODEL_DIM == 3
      else if(it == "k13")
        e->step(i)->add("k13",k13);
      else if(it == "k11")
        e->step(i)->add("k23",k23);
      else if(it == "k33")
        e->step(i)->add("k33",k33);
#endif
    }
    e->save();
  }
  return 0;
}

13. Programming in Feel++

Unresolved directive in 06-programming/README.adoc - include::contributing.adoc[]

13.1. Feel++ Coding Styles

Christophe Prud’homme <@prudhomm> v1.0, 2017/03/07

This is an overview of the coding conventions we use when writing Feel++ code.

13.1.1. Clang Format

clang-format is a powerful tool to reformat your code according to rules defined in a .clang-format file at the toplevel directory of your software.

Feel++ has such a file and define the indentation, space and breaks rules defined later on.

For clang-format to function properly, follow the Comments rules.

To apply Feel++ rules on a file a.cpp in a Feel++ sub-directory, type

clang-format a.cpp

to dump the results of the reformating to the standard output or type

clang-format -i a.cpp

which will replace a.cpp by the reformated file.

be careful when reformating, make sure nobody is working on that file, to avoid creating possibly massive conflicts with the persons currently modifying the same code when they get merged.

13.1.2. Header files and multiple inclusions

To avoid multiple inclusions, wrap every header files using the following technique

// say we have myheader.hpp
#if !defined(FEELPP_MYHEADER_HPP)
#define FEELPP_MYHEADER_HPP 1

// your header here...

#endif // FEELPP_MYHEADER_HPP

more details here

13.1.3. Naming Convention

In Feel++, we basically follow the same naming conventions as in Qt and KDE.

Class names starts with a capital. The rest is in camel case. Function names starts with a lower case, but the first letter of each successive word is capitalized.

Functions and classes should be in the Feel namespace.

The prefix set is used for setters, but the prefix get is not used for accessors. Accessors are simply named with the name of the property they access. The exception is for accessors of a boolean which may start with the prefix is.

Acronyms are lowercased too. Example: Url instead of URL and isFemEnabled() instead of isFEMEnabled()

Accessors should usually be const.

This example shows some possible functions names

class A
{
public:
    void setMatrix(const Matrix& c);
    Matrix matrix() const;
    void setDiagonal(bool b);
    bool isDiagonal() const;
};

13.1.4. Reserved identifiers

13.1.5. Indentation

  • 4 spaces are used for indentation but not in namespace

  • Spaces, not tabs!

  • Suggestion: use emacs and [http://emacswiki.org/emacs/dirvars.el dirvars.el], here is the content of .emacs-dirvars in top Feel++ directory

    indent-tabs-mode: nil
    tab-width: 4
    c-basic-offset: 4
    evaluate: (c-set-offset 'innamespace '0)
    show-trailing-whitespace: t
    indicate-empty-lines: t
    evaluate: (add-hook 'write-file-hooks 'delete-trailing-whitespace)
namespace Feel
{
// no space indentation in namespace
Class A
{
    // 4 spaces indentation
    A() {};
    void f();
};
}

13.1.6. Comments

Use C++ style comment // rather than /* */. It uses less characters but also it is easier to reflow using clang format.

/* Wrong
  the doc
  */

// Correct
// the doc
Doxygen
Doxygen is the tool to document Feel++ and create a reference manual

Use //! to comment function, variables, classes rather than /** */, it allows to reflow comments using clang format.

//!
//! @brief the class
//! @author me <me@email>
//!
class TheClass
{
public:
  //! constructor
  TheClass() {}

private:
  //! member
  int member;
};

//! the function
void thefunction() {}
Feel++ used to promote /** */ but this is no longer the case. The comment style will be updated progressively to match the new style using //!

13.1.7. Declaring variables

  • Declare each variable on a separate line

  • Avoid short (e.g. a, rbarr, nughdeget) names whenever possible

  • Single character variable names are only okay for counters and temporaries, where the purpose of the variable is obvious

  • Wait when declaring a variable until it is needed

// Wrong
int a, b;
char __c, __d;

// Correct
int height;
int width;
char __nameOfThis;
char __nameOfThat;
  • Variables and functions start with a lower-case letter. Each consecutive word in a variable’s or function’s name starts with an upper-case letter

  • Avoid abbreviations

// Wrong
short Cntr; char ITEM_DELIM = '';

// Correct
short counter; char itemDelimiter = '';
  • Classes always start with an upper-case letter.

// Wrong
class meshAdaptation {};

// Correct
class MeshAdaptation {};
  • Non-static data members name of structures and classes always start with M_ . M stands for Member. The rational behind this is for example :

    • to be able to immediately see that the data is a member of a class or a struct

    • to easily search and query-replace

// Wrong
class meshAdaptation { std::vector directions_; };

// Correct
class MeshAdaptation { std::vector M_directions; };
  • Static data members name of structures and classes always start with S_ . S stands for Static. The rational behind this is for example :

    • to be able to immediately see that the data is a static member of a class or a struct

    • to easily search and query-replace

// Wrong
class meshAdaptation { static std::vector directions_; };

// Correct
class MeshAdaptation { static std::vector S_directions; };

13.1.8. Whitespace

  • Use blank lines to group statements together where suited

  • Always use only one blank line

  • Always use a single space after a keyword and before a curly brace.

// Correct
if (foo) { }

// Wrong
if(foo) { }
  • For pointers or references, always use a single space between the type and or &, but no space between the or & and the variable name.

char *x;
const std::string &myString;
const char * const y = "hello";
  • Surround binary operators with spaces.

  • No space after a cast.

  • Avoid C-style casts when possible.

// Wrong
char* blockOfMemory = (char* ) malloc(data.size());

// Correct
char *blockOfMemory = reinterpret_cast(malloc(data.size()));

13.1.9. Braces

  • As a base rule, the left curly brace goes on the same line as the start of the statement:

// Wrong
if (codec) { }

// Correct
if (codec) { }
  • Function implementations and class declarations always have the left brace on the start of a line:

static void foo(int g) { std::cout << g << "" }

class Moo { };
  • Use curly braces when the body of a conditional statement contains more than one line, and also if a single line statement is somewhat complex.

// Wrong
if (address.isEmpty()) { return false; }

for (int i = 0; i < 10; ++i) { std::cout << "i=" << i << ""; }

// Correct
if (address.isEmpty()) return false;

for (int i = 0; i < 10; ++i) std::cout << "=" << i << "";
  • Exception 1: Use braces also if the parent statement covers several lines / wraps

// Correct
if (address.isEmpty() || !isValid() || !codec)
{
    return false;
}
  • Exception 2: Use braces also in if-then-else blocks where either the if-code or the else-code covers several lines

// Wrong
if (address.isEmpty())
  return false;
else
{
  std::cout << address << ""; ++it;
 }

// Correct
if (address.isEmpty())
{
    return false;
}
else
{
    std::cout << address << ""; ++it;
}

// Wrong
if (a) if (b) ... else ...

// Correct
if (a) { if (b) ... else ... }
  • Use curly braces when the body of a conditional statement is empty

// Wrong
while (a);

// Correct
while (a) {}

13.1.10. Parentheses

  • Use parentheses to group expressions:

// Wrong
if (a && b || c)

// Correct
if ((a && b) || c)

// Wrong
a + b & c

// Correct
(a + b) & c

13.1.11. Switch statements

  • The case labels are in the same column as the switch

  • Every case must have a break (or return) statement at the end or a comment to indicate that there’s intentionally no break, unless another case follows immediately.

switch (myEnum)
{
  case Value1:
    doSomething();
    break;
  case Value2:
  case Value3:
    doSomethingElse(); // fall through
  default:
    defaultHandling();
    break;
}

13.1.12. Line breaks

  • Keep lines shorter than 100 characters; insert breaks if necessary.

  • Commas go at the end of a broken line; operators start at the beginning of the new line. An operator at the end of the line is easy to not see if your editor is too narrow.

// Correct
if (longExpression + otherLongExpression + otherOtherLongExpression) { }

// Wrong
if (longExpression + otherLongExpression + otherOtherLongExpression) { }

13.1.13. Inheritance and the virtual keyword

When reimplementing a virtual method, do not put the virtual keyword in the header file.

Quick Reference

14. Quick Reference

In this chapter, we develop a quick reference for the various stages of a simulation using Feel++.

14.1. CMake and Feel++ applications

Feel++ offers a development environment for solving partial differential equations. It uses many tools to cover the different steps pre-processing, processing and post-processing and large range of numerical methods and needs. To this end it is crucial to have a powerful build environment. Feel++ uses CMake from Kitware and provides various macros to help setting up your own application or research project.

14.1.1. CMake macros

Adding a new testcase

For a give application or multiple applications you may define testcases. testcases are difrectory containing a set of files that may include geometry, mesh, cfg or json files.

To define a new testcase case, create a sub-directory where your application, say myapp like in the previous section, stands and copy the required files there.

cd <source directory of my application>
mkdir mytestcase
# copy files (.geo, .msh, .cfg...) to mytestcase
...

then edit the CMakeLists.txt in your application directory and add the following line:

feelpp_add_testcase(mytestcase)

Then type make feelpp_add_testcase_mytestcase in the build directory of your application myapp. It will copy in the build directory of your application the directory mytestcase.

INFO: if you updated the testcase data files, executing make feelpp_testcase_mytestcase will use rsync to update the files that were changed in the source.

The macro feelpp_add_testcase supports options:

PREFIX

(default is feelpp) set the prefix of the target to avoid eg name clash

feelpp_add_testcase(mytestcase PREFIX foo)

then the target is foo_add_testcase_mytestcase.

DEPS

set the dependencies of the testcase

feelpp_add_testcase(mytestcase DEPS myothertestcase)

it allows to update a testcase depending on changes in an other one.

14.2. Setting runtime environment

In this section, we present some tools to initialize and manipulate Feel++ environment.

14.2.1. Initialize Feel++

Environment class is necessary to initialize your application, as seen in FirstApp. Interface is as follows:

Environment env( _argc, _argv, _desc, _about );

None of those parameters are required but it is highly recommended to use the minimal declaration:

Environment env( _argc=argc, _argv=argv,
                 _desc=feel_option(),
                 _about=about(_name="name_of_your_app",
                              _author="your_name",
                              _email="your_email_adress"));
  • _argc and _argv are the arguments of your main function.

  • _desc is a description of your options.

  • _about is a brief description of your application.

14.2.2. Options Description

Adding Options

feel_options() returns a list of default options used in Feel++.

You can create your own list of options as follows:

using namespace Feel;
inline
po::options_description
makeOptions()
{
    po::options_description myappOptions( "My app options" );
    myappOptions.add_options()
      ( "option1", po::value<type1>()->default_value( value1 ), "description1" )
      ( "option2", po::value<type2>()->default_value( value2 ), "description2" )
      ( "option3", po::value<type3>()->default_value( value3 ), "description3" )
      ;
     // Add the default feel options to your list
    return myappOptions.add( feel_options() );
}

makeOptions is the usual name of this routine but you can change it amd myappOptions is the name of you options list.

Parameter

Description

option

the name of parameter

type

the type parameter

value

the default value of parameter

description

the description of parameter

You can then use makeOptions() to initialize the Feel++ Environment as follows

Environment env( _argc=argc, _argv=argv,
                 _desc=makeOptions(),
                 _about=about(_name="myapp",
                              _author="myname",
                              _email="my@email.com") );

Then, at runtime, you can change the parameter as follows

$ ./myapp --option1=alpha --option2=beta --option3=gamma
Excerpt from quickstart/qs_laplacian.cpp
    using namespace Feel;
    using Feel::cout;
        po::options_description laplacianoptions( "Laplacian options" );
        laplacianoptions.add_options()
        ( "no-solve", po::value<bool>()->default_value( false ), "No solve" )
                ;

        Environment env( _argc=argc, _argv=argv,
                   _desc=laplacianoptions,
                   _about=about(_name="qs_laplacian",
                                _author="Feel++ Consortium",
                                _email="feelpp-devel@feelpp.org"));
Options Accessors
Options Description :
Environment::optionsDescription();

Returns options description data structure of type po::options_description po where is a namespace alias defined in Feel++ for boost::program_options.

Variable map

You can access to the parameters of your application environment using the following function:

Environment::vm(_name);

_name is the name of the parameter as seen in the previous paragraph. This function returns a data structure of type po::variable_value. You can then extract the proper parameter value as follows

const double E = Environment::vm(_name="E").template as<double>();
const double nu = Environment::vm(_name="nu").template as<double>();

auto Tfinal =  Environment::vm( _name="test" ).template as<int>()*dt;

14.2.3. Repository

Change changeRepository

You can change the default repository where the results are stored

void changeRepository( _directory, _subdir, _filename );

Parameter

Description

Status

Default value

_directory

directory name

Required

_subdir

true

_filename

logfile

You can use boost format to customize the path as follows:

Environment::changeRepository( boost::format( "doc/manual/laplacian/%1%/%2%-%3%/P%4%/h_%5%/" )
                                   % this->about().appName()
                                   % shape
                                   % Dim
                                   % Order
                                   % meshSize );

Then results will be store in: /doc/manual/laplacian/<appName>/<shape>-<Dim>/P<Order>/h_<meshSize>/

findFile
Interface
std::string findFile( std::string const& filename );

Returns the string containing the filename path.

The lookup is as follows:

  • look into current path

  • look into paths that went through changeRepository(), it means that we look for example into the path from which the executable was run

If the file has an extension .geo or .msh, try also to

  • look into localGeoRepository() which is usually $HOME/feel/geo

  • look into systemGeoRepository() which is usually $FEELPP_DIR/share/feel/geo

If filename is not found, then the empty string is returned.

14.2.4. Utility functions

Communications

A lot of data structures, in fact most of them, in Feel++ are parallel and are associated with a WorldComm data structure which allows us to access and manipulate the MPI communicators. We provide some utility free functions that allow a transparent access to the WorldComm data structure.

We denote by c a Feel++ data structure associated to a WorldComm.

Feel++ Keyword

Description

rank(c)

returns the local MPI rank of the data structure c

globalRank(c)

returns the global MPI rank of the data

For example to print the rank of a mesh data structure

// initialise environment...
auto mesh = makeMesh<Simplex<2,1>>();
std::cout << "local rank : " << rank(mesh) << "\n";

14.3. Using computational meshes

14.3.1. Introduction

Feel++ provides some tools to manipulate meshes. Here is a basic example that shows how to generate a mesh for a square geometry

Excerpt from codes/mymesh.cpp
Unresolved directive in 07-quickref/Mesh/README.adoc - include::../codes/03-mymesh.cpp[tag=mesh]

As always, we initialize the Feel[]+ environment (see section link:[FirstApp] ). The `unitSquare()` will generate a mesh for a square geometry. Feel+ provides several functions to automate the GMSH mesh generation for different topologies. These functions will create a geometry file .geo and a mesh file .msh. We can visualize them in GMSH.

$ gmsh <entity_name>.msh

Finally we use the exporter() (see \ref Exporter) function to export the mesh for post processing. It will create by default a Paraview format file .sos and an Ensight format file .case.

$ paraview <app_name>.case

In this section, we present some of the mesh definition and manipulation tools provided by Feel++. For more information you can also see \ref Gmsh.

14.3.2. Basic Meshes

There is a list of basic geometries you can automatically generate with Feel++ library.

Table 8. Table of function generation basic primitive mesh

Feel++ function

Dim

Description

unitSegment()

1

Build a mesh of the unit segment \(\)[0,1]\(\)

unitSquare()

2

Build a mesh of the unit square \(\)[0,1]^2\(\) using triangles

unitCircle()

2

Build a mesh of the unit circle using triangles

unitHypercube()

3

Build a mesh of the unit hypercube \(\)[0,1]^3\(\) using tetrahedrons

unitSphere()

3

Build a mesh of the unit sphere using tetrahedrons

Examples:

From doc/manual/tutorial/myfunctionspace.cpp

auto mesh = unitSquare();

14.3.3. Load Meshes

loadMesh

You can use this function to:

  • load a .msh file and use the mesh data structure

  • load a .geo file and automatically generate a mesh data structure on this geometrical structure

Interface:

mesh_ptrtype loadMesh(_mesh, _filename, _refine, _update, physical_are_elementary_regions);

Required Parameters:

  • _mesh a mesh data structure.

Optional Parameters:

  • _hsize (double): characteristic size of the mesh. This option will edit the .geo file and change the variable h if defined

    • Default: 0.1

    • Option: gmsh.hsize

  • _geo_variables_list (string): Set a list of variable that may be defined in a .geo file

    • Default: ""

    • Option: gmsh.geo-variables-list

  • _filename (string): filename with extension.

    • Default: feel.geo

    • Option: gmsh.filename

  • _depends (string): list of files (separated by , or ;) on which gmsh.filename depends

    • Default: ""

    • Option: gmsh.depends

  • _refine (boolean): optionally refine with \p refine levels the mesh.

    • Default: 0.

    • Option: gmsh.refine

  • _update (integer): update the mesh data structure (build internal faces and edges).

    • Default: true

  • _physical_are_elementary_regions (boolean): to load specific meshes formats.

    • Default: false.

    • Option: gmsh.physical_are_elementary_regions

  • _straighten (boolean): in case of curvilinear elements, straighten the elements which are not touching with a face the boundary of the domain

    • Default: true

    • Option: gmsh.straighten

  • _partitioner (integer): define the mesh partitioner to use:

    • Default: 1 (if Metis is available) 0 if not (CHACO)

    • Option: gmsh.partitioner

The file you want to load has to be in an appropriate repository. Feel++ looks for .geo and .msh files in the following directories (in this order):

  • current path

  • paths that went through changeRepository(), it means that we look for example into the path from which the executable was run

  • localGeoRepository() which is usually "$HOME/feel/geo" (cf Environment )

  • systemGeoRepository() which is usually "$FEELPP_DIR/share/feel/geo" (cf Environment)

Examples:

Load a mesh data structure from the file "$HOME/feel/mymesh.msh".

auto mesh = loadMesh(_mesh=new mesh_type,
                     _filename="mymesh.msh");

Load a geometric structure from the file ./mygeo.geo and automatically create a mesh data structure.

auto mesh = loadMesh(_mesh=new mesh_type,
                     _filename="mygeo.geo");

Create a mesh data structure from the file ./feel.geo.

auto mesh = loadMesh(_mesh=new Mesh<Simplex< 2 > > );

14.3.4. loadGMSHMesh

In order to load only .msh file, you can also use the loadGMSHMesh.

Interface:

mesh_ptrtype loadGMSHMesh(_mesh, _filename, _refine, _update, _physical_are_elementary_regions);

Required Parameters:

  • _mesh a mesh data structure.

  • _filename filename with extension.

Optional Parameters:

  • _refine optionally refine with \p refine levels the mesh. - Default =0

  • _update update the mesh data structure (build internal faces and edges).

    • Default =true

  • _physical_are_elementary_regions to load specific meshes formats.

    • Default = false

The file you want to load has to be in an appropriate repository. See LoadMesh.

Examples:

From doc/manual/heatns.cpp

 mesh_ptrtype mesh = loadGMSHMesh( _mesh=new mesh_type,
                                   _filename="piece.msh",
                                   _update=MESH_CHECK|MESH_UPDATE_FACES|MESH_UPDATE_EDGES|MESH_RENUMBER );

From applications/check/check.cpp

mesh = loadGMSHMesh( _mesh=new mesh_type,
                     _filename=soption("filename"),
                     _rebuild_partitions=(Environment::worldComm().size() > 1),
                     _update=MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );

14.3.5. Create Meshes

createGMSHMesh

Interface:

mesh_ptrtype createGMSHMesh(_mesh, _desc, _h, _order, _parametricnodes, _refine, _update, _force_rebuild, _physical_are_elementary_regions);

Required Parameters:

  • _mesh mesh data structure.

  • _desc descprition. See further.

Optional Parameters:

  • _h characteristic size.

    • Default = 0.1

  • _order order.

    • Default = 1

  • _parametricnodes

    • Default = 0

  • _refine optionally refine with \p refine levels the mesh.

    • Default =0

  • _update update the mesh data structure (build internal faces and edges).

    • Default =true

  • _force_rebuild rebuild mesh if already exists.

    • Default = false

  • _physical_are_elementary_regions to load specific meshes formats.

    • Default = false

To generate your mesh you need a description parameter. This one can be create by one the two following function.

geo

Use this function to create a description from a .geo file.

Interface:

gmsh_ptrtype geo(_filename, _h, _dim, _order, _files_path);

Required Parameters:

  • filename: file to load.

Optional Parameters:

  • _h characteristic size of the mesh.

    • Default = 0.1.

  • _dim dimension.

    • Default = 3.

  • _order order.

    • Default = 1.

  • _files_path path to the file.

    • Default = localGeoRepository().

The file you want to load has to be in an appropriate repository. See LoadMesh.

Example:

From doc/manual/heat/ground.cpp

mesh = createGMSHMesh( _mesh=new mesh_type,
                       _desc=geo( _filename="ground.geo",
                                  _dim=2,
                                  _order=1,
                                  _h=meshSize ) );

From doc/manual/fd/penalisation.cpp

mesh = createGMSHMesh( _mesh=new mesh_type,
                       _desc=geo( _filename=File_Mesh,
                                  _dim=Dim,
                                  _h=Environment::vm(_name="hsize").template as<double>() ),
                                  _update=MESH_CHECK|MESH_UPDATE_FACES|MESH_UPDATE_EDGES|MESH_RENUMBER );
domain

Use this function to generate a simple geometrical domain from parameters.

Interface:

gmsh_ptrtype domain(_name, _shape, _h, _dim, _order, _convex, \
                    _addmidpoint, _xmin, _xmax, _ymin, _ymax, _zmin, _zmax);

Required Parameters:

  • _name name of the file that will ge generated without extension.

  • _shape shape of the domain to be generated (simplex or hypercube).

Optional Parameters:

  • _h characteristic size of the mesh.

    • Default = 0.1

  • _dim dimension of the domain.

    • Default = 2

  • _order order of the geometry.

    • Default = 1

  • _convex type of convex used to mesh the domain.

    • Default = simplex

  • _addmidpoint add middle point.

    • Default = true

  • _xmin minimum x coordinate.

    • Default = 0

  • _xmax maximum x coordinate.

    • Default = 1

  • _ymin minimum y coordinate.

    • Default = 0

  • _ymax maximum y coordinate.

    • Default = 1.

  • _zmin minimum z coordinate.

    • Default = 0

  • _zmax maximum z coordinate.

    • Default = 1

Example:

From doc/manual/laplacian/laplacian.ccp

mesh_ptrtype mesh = createGMSHMesh( _mesh=new mesh_type,
                                    _desc=domain( _name=( boost::format( "%1%-%2%" ) % shape % Dim ).str() ,
                                                  _usenames=true,
                                                  _shape=shape,
                                                  _h=meshSize,
                                                  _xmin=-1,
                                                  _ymin=-1 ) );

From doc/manual/stokes/stokes.cpp

mesh = createGMSHMesh( _mesh=new mesh_type,
                       _desc=domain( _name=(boost::format("%1%-%2%-%3%")%"hypercube"%convex_type().dimension()%1).str() ,
                                     _shape="hypercube",
                                     _dim=convex_type().dimension(),
                                     _h=meshSize ) );

From doc/manual/solid/beam.cpp

mesh_ptrtype mesh = createGMSHMesh( _mesh=new mesh_type,
                                    _update=MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK,
                                    _desc=domain( _name=( boost::format( "beam-%1%" ) % nDim ).str(),
                                                  _shape="hypercube",
                                                  _xmin=0., _xmax=0.351,
                                                  _ymin=0., _ymax=0.02,
                                                  _zmin=0., _zmax=0.02,
                                                  _h=meshSize ) );
Explanation pending on straightenMesh

14.3.6. Mesh iterators

Feel++ mesh data structure allows to iterate over its entities: elements, faces, edges and points.

The following table describes free-functions that allow to define mesh region over which to operate. MeshType denote the type of mesh passed to the free functions in the table.

MeshType can be a pointer, a shared_pointer or a reference to a mesh type.

For example :

auto mesh = loadMesh( _mesh=Mesh<Simplex<2>>);
auto r1 = elements(mesh); // OK
auto r2 = elements(*mesh); // OK
Table 9. Table of mesh iterators
Type Function Description

elements_t<MeshType>

elements(mesh)

All the elements of a mesh

markedelements_t<MeshType>

markedelements(mesh, id)

All the elements marked by marked id

boundaryelements_t<MeshType>

boundaryelements(mesh)

All the elements of the mesh which share a face with the boundary of the mesh.

internalelements_t<MeshType>

internalelements(mesh)

All the elements of the mesh which share a face with the boundary of the mesh.

pid_faces_t<MeshType>

faces(mesh)

All the faces of the mesh.

markedfaces_t<MeshType>

markedfaces(mesh)

All the faces of the mesh which are marked.

boundaryfaces_t<MeshType>

boundaryfaces(mesh)

All elements that own a topological dimension one below the mesh. For example, if you mesh is a 2D one, boundaryfaces(mesh) will return all the lines (because of dimension \(\)2-1=1\(\)). These elements which have one dimension less, are corresponding to the boundary faces.

internalfaces_t<MeshType>

internalelements(mesh)

All the elements of the mesh which are stricly within the domain that is to say they do not share a face with the boundary.

edges_t<MeshType>

edges(mesh)

All the edges of the mesh.

boundaryedges_t<MeshType>

boundaryedges(mesh)

All boundary edges of the mesh.

points_t<MeshType>

points(mesh)

All the points of the mesh.

markedpoints_t<MeshType>

markedpoints(mesh,id)

All the points marked id of mesh.

boundarypoints_t<MeshType>

boundarypoints(mesh)

All boundary points of the mesh.

internalpoints_t<MeshType>

internalpoints(mesh)

All internal points of the mesh(not on the boundary)

Here are some examples on how to use these functionSpace

auto mesh = ...;

auto r1 = elements(mesh);
// iterate over the set of elements local to the process(no ghost cell selected, see next section)
for ( auto const&  e : r2 )
{
  auto const& elt = unwrap_ref( e );
  // work with element elt
  ...
}

auto r2 = markedelements(mesh,"iron");
// iterate over the set of elements marked iron in the mesh
for ( auto const&  e : r2 )
{
  auto const& elt = unwrap_ref( e );
  // work with element elt
  ...
}

auto r3 = boundaryfaces(mesh);
// iterate over the set of faces on the boundary of the mesh
for ( auto const&  e : r3 )
{
  auto const& elt = unwrap_ref( e );
  // work with element elt
  ...
}

auto r4 = markededges(mesh,"line");
// iterate over the set of edges marked "line" in the mesh
for ( auto const&  e : r4 )
{
  auto const& elt = unwrap_ref( e );
  // work with element elt
  ...
}
Extended set of entities

Feel++ allows also to select an extended sets of entities from the mesh, you can extract entities which belongs to the local process but also ghost entities which satisfy the same property as the local ones.

Actually you can select both or one one of them thanks to the enum data structure entity_process_t which provides the following options

entity_process_t

Description

LOCAL_ONLY

only local entities

GHOST_ONLY

only ghost entities

ALL

both local and ghost entities

Type

Function

Description

ext_elements_t<MeshType>

elements(mesh,entity_process_t)

all the elements of mesh associated to entity_process_t.

ext_elements_t<MeshType>

markedelements(mesh, id, entity_process_t)

all the elements marked id of mesh associated to entity_process_t.

ext_faces_t<MeshType>

faces(mesh,entity_process_t)

all the faces of mesh associated to entity_process_t.

ext_faces_t<MeshType>

markedfaces(mesh, id, entity_process_t)

all the faces marked id of mesh associated to entity_process_t.

ext_edges_t<MeshType>

edges(mesh,entity_process_t)

all the edges of mesh associated to entity_process_t.

ext_edges_t<MeshType>

markededges(mesh, id, entity_process_t)

all the edges marked id of mesh associated to entity_process_t.

The type of the object returned for an entity is always the same, for elements it is ext_elements_t<MeshType> whether the elements are marked or not. The reason is that in fact we have to create a temporary data structure embedded in the range object that stores a reference to the elements which are selected.

Here is how to select both local and ghost elements from a Mesh

auto mesh =...;
auto r = elements(mesh,entity_process_t::ALL);
for (auto const& e : r )
{
  // do something on the local and ghost element
  ...
  // do something special on ghost cells
  if ( unwrap_ref(e).isGhostCell() )
  {...}
}
Concatenate sets of entities

Denote \(\mathcal{E}_{1}, \ldots ,\mathcal{E}_{n}\) \(n\) disjoints sets of the same type of entities (eg elements, faces,edges or points), \(\cup_{i=1}^{n} \mathcal{E}_i\) with \(\cap_{i=0}^{n} \mathcal{E}_i = \emptyset\).

We wish to concatenate these \(n\) sets. To this end, we use concatenate which takes an arbitrary number of disjoints sets.

#include <feel/feelmesh/concatenate.hpp>
...
auto E_1 = internalfaces(mesh);
auto E_2 = markedfaces(mesh,"Gamma_1");
auto E_3 = markedfaces(mesh,"Gamma_2");
auto newset = concatenate( E_1, E_2, E_3 );
cout << "measure of newset = " << integrate(_range=newset, _expr=cst(1.)).evaluate() << std::endl;
Compute the complement of a set of entities

Denote \(\mathcal{E}\) a set of entities, eg. the set of all faces (both internal and boundary faces). Denote \(\mathcal{E}_\Gamma\) a set of entities marked by \(\Gamma\). We wish to build \({\Gamma}^c=\mathcal{E}\backslash\Gamma\). To compute the complement, Feel++ provides a complement template function that requires \(\mathcal{E}\) and a predicate that return true if an entity of \(\mathcal{E}\) belongs to \(\Gamma\), false otherwise. The function returns mesh iterators over \(\Gamma^c\).

#include <feel/feelmesh/complement.hpp>
...
auto E = faces(mesh);
// build set of boundary faces, equivalent to boundaryfaces(mesh)
auto bdyfaces = complement(E,[](auto const& e){return e.isOnBoundary()});
cout << "measure of bdyfaces = " << integrate(_range=bdyfaces, _expr=cst(1.)).evaluate() << std::endl;
// should be the same as above
cout << "measure of boundaryfaces = " << integrate(_range=boundaryfaces(mesh), _expr=cst(1.)).evaluate() << std::endl;
Helper function on entities set

Feel++ provides some helper functions to apply on set of entities. We denote by range_t the type of the entities set.

Table 10. Utility functions

Type

Function

Description

size_type

nelements(range_t,bool)

returns the local number of elements in entities set range_t of bool is false, other the global number which requires communication (default: global number)

WorldComm

worldComm(range_t)

returns the WorldComm associated to the entities set

Create a new range

A range can be also build directly by the user. This customized range is stored in a std container which contains the c++ references of entity object. We use boost::reference_wrapper for take c++ references and avoid copy of mesh data. All entities enumerated in the range must have same type (elements,faces,edges,points). Below we have an example which select all active elements in mesh for the current partition (i.e. identical to elements(mesh)).

auto mesh = ...;
// define reference entity type
typedef boost::reference_wrapper<typename mesh_type::element_type const> element_ref_type;
// store entities in a vector
typedef std::vector<element_ref_type> cont_range_type;
boost::shared_ptr<cont_range_type> myelts( new cont_range_type );
for (auto const& elt : elements(mesh) )
{
    myelts->push_back(boost::cref(elt));
}
// generate a range object usable in feel++
auto myrange = boost::make_tuple( mpl::size_t<MESH_ELEMENTS>(),
                                  myelts->begin(),myelts->end(),myelts );

Next, this range can be used in feel++ language.

double eval = integrate(_range=myrange,_expr=cst(1.)).evaluate()(0,0);

14.3.7. Mesh Markers

Elements and their associated sub-entities can be marked.

A marker is an integer specifying for example a material id, a boundary condition id or some other property associated with the entity.

A dictionary can map string to marker ids.

The dictionary is stored in the Mesh data structures and provides the set of correspondances between strings and ids.

To access a marker, it is necessary to verify that it exists as follows

for( auto const& ewrap : elements(mesh))
{
  auto const& e = unwrap_ref( ewrap );
  if ( e.hasMarker() ) (1)
  {
    std::cout << "Element " << e.id() << " has marker " << e.marker() << std::endl;
  }
  if ( e.hasMarker(5) ) (2)
  {
    std::cout << "Element " << e.id() << " has marker 5 " << e.marker(5) << std::endl;
  }
}
1 check if marker 1 (the default marker) exists, if yes then print it
2 check if marker 5 exists, if yes then print it

14.4. Integration

You should be able to create a mesh now. If it is not the case, get back to the section Mesh.

14.4.1. Integrals

Feel++ provide the integrate() function to define integral expressions which can be used to compute integrals, define linear and bi-linear forms.

Interface
  integrate( _range, _expr, _quad, _geomap );

Please notice that the order of the parameter is not important, these are boost parameters, so you can enter them in the order you want. To make it clear, there are two required parameters and 2 optional and they of course can be entered in any order provided you give the parameter name. If you don’t provide the parameter name (that is to say _range = or the others) they must be entered in the order they are described below.

Required parameters:

  • _range = domain of integration

  • _expr = integrand expression

Optional parameters:

  • _quad = quadrature to use instead of the default one, wich means _Q<integer>() where the integer is the polynomial order to integrate exactely

  • _geomap = type of geometric mapping to use, that is to say:

Feel Parameter

Description

GEOMAP_HO

High order approximation (same of the mesh)

GEOMAP_OPT

Optimal approximation: high order on boundary elements order 1 in the interior

GEOMAP_01

Order 1 approximation (same of the mesh)

Example

From doc/manual/tutorial/dar.cpp

form1( ... ) = integrate( _range = elements( mesh ),
                          _expr = f*id( v ) );

From doc/manual/tutorial/myintegrals.cpp

  // compute integral f on boundary
  double intf_3 = integrate( _range = boundaryfaces( mesh ),
                             _expr = f );

From doc/manual/advection/advection.cpp

  form2( _test = Xh, _trial = Xh, _matrix = D ) +=
    integrate( _range = internalfaces( mesh ),
               _quad = _Q<2*Order>(),
               _expr = ( averaget( trans( beta )*idt( u ) ) * jump( id( v ) ) )
               + penalisation*beta_abs*( trans( jumpt( trans( idt( u ) )) )
               *jump( trans( id( v ) ) ) ),
               _geomap = geomap );

From doc/manual/laplacian/laplacian.cpp

 auto l = form1( _test=Xh, _vector=F );
 l = integrate( _range = elements( mesh ),
                _expr=f*id( v ) ) +
     integrate( _range = markedfaces( mesh, "Neumann" ),
                _expr = nu*gradg*vf::N()*id( v ) );

14.4.2. Computing my first Integrals

This part explains how to integrate on a mesh with Feel++ (source doc/manual/tutorial/myintegrals.cpp ).

Let’s consider the domain \(\Omega=[0,1\)^d] and associated meshes. Here, we want to integrate the following function

\[\begin{aligned} f(x,y,z) = x^2 + y^2 + z^2 \end{aligned}\]

on the whole domain \(\Omega\) and on part of the boundary \(\Omega\).

There is the appropriate code:

int
main( int argc, char** argv )
{
    // Initialize Feel++ Environment
    Environment env( _argc=argc, _argv=argv,
                     _desc=feel_options(),
                     _about=about( _name="myintegrals" ,
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org" ) );

    // create the mesh (specify the dimension of geometric entity)
    auto mesh = unitHypercube<3>();

    // our function to integrate
    auto f = Px()*Px() + Py()*Py() + Pz()*Pz();

    // compute integral of f (global contribution)
    double intf_1 = integrate( _range = elements( mesh ),
                               _expr = f ).evaluate()( 0,0 );

    // compute integral of f (local contribution)
    double intf_2 = integrate( _range = elements( mesh ),
                               _expr = f ).evaluate(false)( 0,0 );

    // compute integral f on boundary
    double intf_3 = integrate( _range = boundaryfaces( mesh ),
                               _expr = f ).evaluate()( 0,0 );

    std::cout << "int global ; local ; boundary" << std::endl
              << intf_1 << ";" << intf_2 << ";" << intf_3 << std::endl;
}

14.4.3. Mean value of a function

Let \(f\) a bounded function on domain \(\Omega\). You can evaluate the mean value of a function thanks to the mean() function :

\[\bar{f}=\frac{1}{|\Omega|}\int_\Omega f=\frac{1}{\int_\Omega 1}\int_\Omega f\]
Interface
  mean( _range, _expr, _quad, _geomap );

Required parameters:

  • _range = domain of integration

  • _expr = mesurable function

Optional parameters:

  • _quad = quadrature to use.

    • Default = _Q<integer>()

  • _geomap = type of geometric mapping.

    • Default = GEOMAP_OPT

Example
Stokes example using mean
int main(int argc, char**argv )
{
    Environment env( _argc=argc, _argv=argv,
                     _about=about(_name="mystokes",
                                  _author="Feel++ Consortium",
                                  _email="feelpp-devel@feelpp.org"));

    // create the mesh
    auto mesh = loadMesh(_mesh=new Mesh<Simplex< 2 > > );


    // function space
    auto Vh = THch<2>( mesh );

    // element U=(u,p) in Vh
    auto U = Vh->element();
    auto u = U.element<0>();
    auto p = U.element<1>();

    // left hand side
    auto a = form2( _trial=Vh, _test=Vh );
    a = integrate(_range=elements(mesh),
                  _expr=trace(gradt(u)*trans(grad(u))) );

    a+= integrate(_range=elements(mesh),
                  _expr=-div(u)*idt(p)-divt(u)*id(p));

    auto syms = symbols<2>();
    auto u1 = parse( option(_name="functions.alpha").as<std::string>(), syms );
    auto u2 = parse( option(_name="functions.beta").as<std::string>(), syms );
    matrix u_exact = matrix(2,1);
    u_exact = u1,u2;
    auto p_exact = parse( option(_name="functions.gamma").as<std::string>(), syms );
        auto f = -laplacian( u_exact, syms ) + grad( p_exact, syms ).transpose();
    LOG(INFO) << "rhs : " << f;

    // right hand side
    auto l = form1( _test=Vh );
    l = integrate(_range=elements(mesh),
                  _expr=trans(expr<2,1,5>( f, syms ))*id(u));
    a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
          _expr=expr<2,1,5>(u_exact,syms));

    // solve a(u,v)=l(v)
    a.solve(_rhs=l,_solution=U);

    double mean_p = mean(_range=elements(mesh),_expr=idv(p))(0,0);
    double mean_p_exact = mean(_range=elements(mesh),_expr=expr(p_exact,syms))(0,0);
    double l2error_u = normL2( _range=elements(mesh), _expr=idv(u)-expr<2,1,5>( u_exact, syms ) );
    double l2error_p = normL2( _range=elements(mesh), _expr=idv(p)-mean_p-(expr( p_exact, syms )-mean_p_exact) );
    LOG(INFO) << "L2 error norm u: " << l2error_u;
    LOG(INFO) << "L2 error norm p: " << l2error_p;

    // save results
    auto e = exporter( _mesh=mesh );
    e->add( "u", u );
    e->add( "p", p );
    e->save();
}

14.4.4. Norms

Let \(f\) a bounded function on domain \(\Omega\).

L2 norms

Let \(\)f \in L^2(\Omega)\(\) you can evaluate the \(\)L^2\(\) norm using the normL2() function:

\(\) \parallel f\parallel_{L^2(\Omega)}=\sqrt{\int_\Omega |f|^2} \(\)

Interface
  normL2( _range, _expr, _quad, _geomap );

or squared norm:

  normL2Squared( _range, _expr, _quad, _geomap );

Required parameters:

  • _range = domain of integration

  • _expr = mesurable function

Optional parameters:

  • _quad = quadrature to use.

    • Default = _Q<integer>()

  • _geomap = type of geometric mapping.

    • Default = GEOMAP_OPT

Example

From doc/manual/laplacian/laplacian.cpp

  double L2error =normL2( _range=elements( mesh ),
                          _expr=( idv( u )-g ) );

From doc/manual/stokes/stokes.cpp

Stokes example using mean
int main(int argc, char**argv )
{
    Environment env( _argc=argc, _argv=argv,
                     _about=about(_name="mystokes",
                                  _author="Feel++ Consortium",
                                  _email="feelpp-devel@feelpp.org"));

    // create the mesh
    auto mesh = loadMesh(_mesh=new Mesh<Simplex< 2 > > );


    // function space
    auto Vh = THch<2>( mesh );

    // element U=(u,p) in Vh
    auto U = Vh->element();
    auto u = U.element<0>();
    auto p = U.element<1>();

    // left hand side
    auto a = form2( _trial=Vh, _test=Vh );
    a = integrate(_range=elements(mesh),
                  _expr=trace(gradt(u)*trans(grad(u))) );

    a+= integrate(_range=elements(mesh),
                  _expr=-div(u)*idt(p)-divt(u)*id(p));

    auto syms = symbols<2>();
    auto u1 = parse( option(_name="functions.alpha").as<std::string>(), syms );
    auto u2 = parse( option(_name="functions.beta").as<std::string>(), syms );
    matrix u_exact = matrix(2,1);
    u_exact = u1,u2;
    auto p_exact = parse( option(_name="functions.gamma").as<std::string>(), syms );
        auto f = -laplacian( u_exact, syms ) + grad( p_exact, syms ).transpose();
    LOG(INFO) << "rhs : " << f;

    // right hand side
    auto l = form1( _test=Vh );
    l = integrate(_range=elements(mesh),
                  _expr=trans(expr<2,1,5>( f, syms ))*id(u));
    a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
          _expr=expr<2,1,5>(u_exact,syms));

    // solve a(u,v)=l(v)
    a.solve(_rhs=l,_solution=U);

    double mean_p = mean(_range=elements(mesh),_expr=idv(p))(0,0);
    double mean_p_exact = mean(_range=elements(mesh),_expr=expr(p_exact,syms))(0,0);
    double l2error_u = normL2( _range=elements(mesh), _expr=idv(u)-expr<2,1,5>( u_exact, syms ) );
    double l2error_p = normL2( _range=elements(mesh), _expr=idv(p)-mean_p-(expr( p_exact, syms )-mean_p_exact) );
    LOG(INFO) << "L2 error norm u: " << l2error_u;
    LOG(INFO) << "L2 error norm p: " << l2error_p;

    // save results
    auto e = exporter( _mesh=mesh );
    e->add( "u", u );
    e->add( "p", p );
    e->save();
}
H1 norm

In the same idea, you can evaluate the H1 norm or semi norm, for any function \(f \in H^1(\Omega)\):

\[\begin{aligned} \parallel f \parallel_{H^1(\Omega)}&=\sqrt{\int_\Omega |f|^2+|\nabla f|^2}\\ &=\sqrt{\int_\Omega |f|^2+\nabla f * \nabla f^T}\\ |f|_{H^1(\Omega)}&=\sqrt{\int_\Omega |\nabla f|^2} \end{aligned}\]

where \(*\) is the scalar product \(\cdot\) when \(f\) is a scalar field and the frobenius scalar product \(:\) when \(f\) is a vector field.

Interface
  normH1( _range, _expr, _grad_expr, _quad, _geomap );

or semi norm:

  normSemiH1( _range, _grad_expr, _quad, _geomap );

Required parameters:

  • _range = domain of integration

  • _expr = mesurable function

  • _grad_expr = gradient of function (Row vector!)

Optional parameters:

  • _quad = quadrature to use.

    • Default = _Q<integer>()

  • _geomap = type of geometric mapping.

    • Default = GEOMAP_OPT

normH1() returns a float containing the \(\)H^1\(\) norm.

Example

With expression:

  auto g = sin(2*pi*Px())*cos(2*pi*Py());
  auto gradg = 2*pi*cos(2* pi*Px())*cos(2*pi*Py())*oneX()
            -2*pi*sin(2*pi*Px())*sin(2*pi*Py())*oneY();
// There gradg is a column vector!
// Use trans() to get a row vector
  double normH1_g = normH1( _range=elements(mesh),
                     _expr=g,
                   _grad_expr=trans(gradg) );

With test or trial function u

  double errorH1 = normH1( _range=elements(mesh),
                    _expr=(u-g),
                  _grad_expr=(gradv(u)-trans(gradg)) );
\(L^\infty\) norm

You can evaluate the infinity norm using the normLinf() function:

\[\parallel f \parallel_\infty=\sup_\Omega(|f|)\]
Interface
  normLinf( _range, _expr, _pset, _geomap );

Required parameters:

  • _range = domain of integration

  • _expr = mesurable function

  • _pset = set of points (e.g. quadrature points)

Optional parameters:

  • _geomap = type of geometric mapping.

    • Default = GEOMAP_OPT

The normLinf() function returns not only the maximum of the function over a sampling of each element thanks to the _pset argument but also the coordinates of the point where the function is maximum. The returned data structure provides the following interface

  • value(): return the maximum value

  • operator()(): synonym to value()

  • arg(): coordinates of the point where the function is maximum

Example
  auto uMax = normLinf( _range=elements(mesh),
                        _expr=idv(u),
                        _pset=_Q<5>() );
  std::cout << "maximum value : " << uMax.value() << std::endl
            <<  "         arg : " << uMax.arg() << std::endl;

14.5. Function Spaces

Prerequisites

The prerequisites are

14.5.1. Notations

We now turn to the next crucial mathematical ingredient: the function space, whose definition depends on \(\Omega_h\) - or more precisely its partitioning \(\mathcal{T}_h\) - and the choice of basis function. Function spaces in Feel++ follow the same definition and Feel++ provides support for continuous and discontinuous Galerkin methods and in particular approximations in \(L^2\), \(H^1\)-conforming and \(H^1\)-nonconforming, \(H^2\), \(H(\mathrm{div})\) and \(H(\mathrm{curl})\)[^1].

We introduce the following spaces

\[\begin{aligned} \mathbb{W}_h &= \{v_h \in L^2(\Omega_h): \ \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{P}_K\},\\ \mathbb{V}_h &= \mathbb{W}_h \cap C^0(\Omega_h)= \{ v_h \in \mathbb{W}_h: \ \forall F \in \mathcal{F}^i_h\ [ v_h ]_F = 0\}\\ \mathbb{H}_h &= \mathbb{W}_h \cap C^1(\Omega_h)= \{ v_h \in \mathbb{W}_h: \ \forall F \in \mathcal{F}^i_h\ [ v_h ]_F = [ \nabla v_h ]_F = 0\}\\ \mathbb{C}\mathbb{R}_h &= \{ v_h \in L^2(\Omega_h):\ \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{P}_1; \forall F \in \mathcal{F}^i_h\ \int_F [ v_h ] = 0 \}\\ \mathbb{R}{a}\mathbb{T}{u}_h &= \{ v_h \in L^2(\Omega_h):\ \forall K \in \mathcal{T}_h, v_h|_K \in \mathrm{Span}\{1,x,y,x^2-y^2\}; \forall F \in \mathcal{F}^i_h\ \int_F [ v_h ] = 0 \}\\ \mathbb{R}\mathbb{T}_h &= \{\mathbf{v}_h \in [L^2(\Omega_h)]^d:\ \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{R}\mathbb{T}_k; \forall F \in \mathcal{F}^i_h\ [{\mathbf{v}_h \cdot \mathrm{n}}]_F = 0 \}\\ \mathbb{N}_h &= \{\mathbf{v}_h \in [L^2(\Omega_h)]^d:\ \forall K \in \mathcal{T}_h, v_h|_K \in \mathbb{N}_k; \forall F \in \mathcal{F}^i_h\ [{\mathbf{v}_h \times \mathrm{n}}]_F = 0 \} \end{aligned}\]

where \(\mathbb{R}\mathbb{T}_k\) and \(\mathbb{N}_k\) are respectively the Raviart-Thomas and Nédélec finite elements of degree \(k\).

The Legrendre and Dubiner basis yield implicitely discontinuous approximations, the Legendre and Dubiner boundary adapted basis, see~\cite MR1696933, are designed to handle continuous approximations whereas the Lagrange basis can yield either discontinuous or continuous (default behavior) approximations. \(\mathbb{R}\mathbb{T}_h\) and \(\mathbb{N}_h\) are implicitely spaces of vectorial functions \(\mathbf{f}\) such that \(\mathbf{f}: \Omega_h \subset \mathbb{R}^d \mapsto \mathbb{R}^d\). As to the other basis functions, i.e. Lagrange, Legrendre, Dubiner, etc., they are parametrized by their values namely Scalar, Vectorial or Matricial.

Products of function spaces must be supported. This is very powerful to describe complex multiphysics problems when coupled with operators, functionals and forms described in the next section. Extracting subspaces or component spaces are part of the interface.
Function Spaces

Function spaces support is provided by the FunctionSpace class

The FunctionSpace class

  • constructs the table of degrees of freedom which maps local (elementwise) degrees of freedom to the global ones with respect to the geometrical entities,

  • embeds the definition of the elements of the function space allowing for a tight coupling between the elements and their function spaces,

  • stores an interpolation data structure (e.g. region tree) for rapid localisation of point sets (determining in which element they reside).

C++ Function

C++ Type

Function Space [1]

Pch<N>(mesh)

Pch_type<MeshType,N>

\(P^N_{c,h}\)

Pchv<N>(mesh)

Pchv_type<MeshType,N>

\([P^N_{c,h}\)^d]

Pdh<N>(mesh)

Pdh_type<MeshType,N>

\(P^N_{d,h}\)

Pdhv<N>(mesh)

Pdhv_type<MeshType,N>

\([P^N_{d,h}\)^d]

THch<N>(mesh)

THch_type<MeshType,N>

\([P^{N+1}_{c,h}\)^d \times P^N_{c,h}]

Dh<N>(mesh)

Dh_type<MeshType,N>

\(\mathbb{R}\mathbb{T}_h\)

Ned1h<N>(mesh)

Ned1h_type<MeshType,N>

\(\mathbb{N}_h\)

[1]: see Notations for the function spaces definitions.

Here are some examples how to define function spaces with Lagrange basis functions.

#include <feel/feeldiscr/pch.hpp>

// Mesh with triangles
using MeshType = Mesh<Simplex<2>>;
// Space spanned by P3 Lagrange finite element
FunctionSpace<MeshType,bases<Lagrange<3>>> Xh;
// is equivalent to (they are the same type)
Pch_type<MeshType,3> Xh;

// using the auto keyword
MeshType mesh = loadMesh( _mesh=new MeshType );
auto Xh = Pch<3>( mesh );
// is equivalent to
auto Xh = FunctionSpace<MeshType,bases<Lagrange<3>>>::New( mesh );
auto Xh = Pch_type<MeshType,3>::New( mesh );
Functions
One important feature in FunctionSpace is that it embeds the definition of element which allows for the strict definition of an Element of a FunctionSpace and thus ensures the correctness of the code.

An element has its representation as a vector, also in the case of product of multiple spaces.

#include <feel/feeldiscr/pch.hpp>

// Mesh with triangles
using MeshType = Mesh<Simplex<2>>;
auto mesh = loadMesh( _mesh=new MeshType );

// define P3 Lagrange finite element space
auto P3ch = Pch<3>(mesh);

// definie an element from P3ch, initialized to 0
auto u = P3ch.element();
// definie an element from P3ch, initialized to x^2+y^2
auto v = P3ch.element(Px()*Px()+Py()*Py());
Components
FunctionSpace<Mesh<Simplex<2> >,
 bases<Lagrange<2,Vectorial>, Lagrange<1,Scalar>,
       Lagrange<1,Scalar> > > P2P1P1;
auto U = P2P1P1.element();
// Views: changing a view changes U and vice versa
// view on element associated to P2
auto u = U.element<0>();
// extract view of first component
auto ux = u.comp(X);
// view on element associated to 1st P1
auto p = U.element<1>();
// view on element associated to 2nd P1
auto q = U.element<2>();

14.5.2. Interpolation

Feel++ has a very powerful interpolation framework which allows to:

  • transfer functions from one mesh to another

  • transfer functions from one space type to another.

this is done seamlessly in parallel. The framework provides a set of C++ classes and C++ free-functions enabled short, concise and expressive handling of interpolation.

Using interpolation operator
Building interpolation operator \(\)I_h : P^1_{c,h} \rightarrow P^0_{td.h}\(\)
using MeshType = Mesh<Simplex<2>>;
auto mesh loadMesh( _mesh=new MeshType );
auto P1h = Pch<1>( mesh );
auto P0h = Pdh<0>( mesh );
auto Ih = I( _domain=P1h, _image=P0h );
De Rahm Diagram

The De Rahm diagram reads as follows: the range of each of the operators coincides with the null space of the next operator in the sequence below, and the last map is a surjection.

\[\begin{array}{ccccccc} H^1(\Omega)& \overset{\nabla}{\longrightarrow}& H^{\mathrm{curl}}(\Omega)& \overset{\nabla \times}{\longrightarrow}& H^{\mathrm{div}}(\Omega)& \overset{\nabla \cdot}{\longrightarrow}& L^2(\Omega) \end{array}\]

An important result is that the diagram transfers to the discrete level

\[\begin{array}{ccccccc} H^1(\Omega)& \overset{\nabla}{\longrightarrow}& H^{\mathrm{curl}}(\Omega)& \overset{\nabla \times}{\longrightarrow}& H^{\mathrm{div}}(\Omega)& \overset{\nabla \cdot}{\longrightarrow}& L^2(\Omega) \\ \left\downarrow\right.\pi_{c,h}& ~ & \left\downarrow\right.\pi_{\mathrm{curl},h}& ~ & \left\downarrow\right.\pi_{\mathrm{div},_h}& ~ & \left\downarrow\right.\pi_{d,h}& ~ \\ U_h& \overset{\nabla}{\longrightarrow}& V_h& \overset{\nabla \times}{\longrightarrow}& W_h& \overset{\nabla \cdot}{\longrightarrow}& Z_h\\ \end{array}\]

The diagram above is commutative which means that we have the following properties:

\[\begin{aligned} \nabla(\pi_{c,h} u) &= \pi_{\mathrm{curl},h}( \nabla u ),\\ \nabla\times(\pi_{\mathrm{curl},h} u) &= \pi_{\mathrm{div},h}( \nabla\times u ),\\ \nabla\cdot(\pi_{\mathrm{div},h} u) &= \pi_{d,h}( \nabla\cdot u ) \end{aligned}\]
The diagram can be restricted to functions satisfying the homogeneous Dirichlet boundary conditions
\[\begin{array}{ccccccc} H^1_0(\Omega)& \overset{\nabla}{\longrightarrow}& H_0^{\mathrm{curl}}(\Omega)& \overset{\nabla \times}{\longrightarrow}& H_0^{\mathrm{div}}(\Omega)& \overset{\nabla \cdot}{\longrightarrow}& L^2_0(\Omega) \end{array}\]

Interpolation operators are provided as is or as shared pointers. The table below presents the alternatives.

Table 11. Table of Interpolation operators

C++ object

C++ Type

C++ shared object

C++ Type

Mathematical operator

I(_domain=Xh,_image=Yh)

I_t<functionspace_type<decltype(Xh)>, functionspace_type<decltype(Xh)>>

IPtr(…​)

I_ptr_t<functionspace_type<decltype(Xh)>, functionspace_type<decltype(Xh)>>

\(\)I: X_h \rightarrow Y_h \(\)

Grad(_domain=Xh,_image=Wh)

Grad_t<functionspace_type<decltype(Xh)>, functionspace_type<decltype(Wh)>>

GradPtr(…​)

Grad_ptr_t<functionspace_type<decltype(Xh)>, functionspace_type<decltype(Wh)>>

\(\)\nabla: X_h \rightarrow W_h \(\)

Curl(_domain=Wh,_image=Vh)

Curl_t<functionspace_type<decltype(Wh)>, functionspace_type<decltype(Vh)>>

CurlPtr(…​)

Curl_ptr_t<functionspace_type<decltype(Wh)>, functionspace_type<decltype(Vh)>>

\(\)\nabla \times : W_h \rightarrow V_h \(\)

Div(_domain=Vh,_image=Zh)

Div_t<functionspace_type<decltype(Vh)>, functionspace_type<decltype(Zh)>>

DivPtr(…​)

Div_ptr_t<functionspace_type<decltype(Vh)>, functionspace_type<decltype(Zh)>>

\(\)\nabla \cdot: V_h \rightarrow Z_h \(\)

Building the discrete operators associated to the De Rahm diagram in Feel++
auto mesh = loadMesh( _mesh=new Mesh<Simplex<Dim>>());
auto Xh = Pch<1>(mesh);
auto Gh = Ned1h<0>(mesh);
auto Ch = Dh<0>(mesh);
auto P0h = Pdh<0>(mesh);
auto Igrad = Grad( _domainSpace = Xh, _imageSpace=Gh );
auto Icurl = Curl( _domainSpace = Gh, _imageSpace=Ch );
auto Idiv = Div( _domainSpace = Ch, _imageSpace=P0h );

auto u = Xh->element(<expr>);
auto w = Igrad(u); // w in Gh
auto x = Icurl(w); // z in Ch
auto y = Idiv(x); // y in P0h

14.5.3. Saving and loading functions on disk

Saving functions on disk

To save a function on disk to use it later, for example in another application, you can use the save function.

The saved file will be named after the name registered for the variable in the constructor (default : u).

auto Vh = Pch<1>( mesh );
auto u = Vh->element("v");
// do something with u
...
// save /path/to/save/v.h5
u.save( _path="/path/to/save", _type="hdf5" );

The path parameter creates a directory at this path to store all the degrees of liberty of this function.

The type parameter can be binary, text or hdf5 . The first two will create one file per processor, whereas "hdf5" will creates only one file.

Loading functions from disk
To load a function, the mesh need to be exactly the same as the one used when saving it.
auto Vh = Pch<1>( mesh );
auto u = Vh->element("v");
// load /path/to/load/v.h5
u.load( _path="/path/to/load/", _type="hdf5" );

The path and type parameters need to be the same as the one used to save the function.

Extended parallel doftable

In some cases, when we use parallel data, informations from other interfaces of partitions are need. To manage this, we can add ghost degree of freedom on ghost elements at these locations. However, we have to know if data have extended parallel doftable to load and use it.

In order to pass above this restriction, the two function load and save has been updated to use hdf5 format. With this format, extended parallel doftable or not, the function will work without any issues. More than that, we can load elements with extended parallel doftable and resave it without it, and vice versa. This last feature isn’t available with other formats than hdf5.

14.6. Bilinear and Linears Forms

We consider in this section bilinear and linear forms \(a: X_h \times X_h \rightarrow \mathbb{R}\) and \(\ell: X_h \rightarrow \mathbb{R}.\)

We suppose in this section that you know how to define your Mesh and your function spaces. You may need integration tools too, see Integrals.

There are Feel++ tools you need to create linear and bilinear forms in order to solve variational formulation.

from now on, u denotes an element from your trial function space (unknown function) and v an element from your test function space

14.6.1. Building Forms

Using form1

To construct a linear form \(\)\ell: X_h \rightarrow \mathbb{R}\(\), proceed as follows

auto mesh = ...;
// build a P1/Q1 approximation space
auto Xh = Pch<1>( mesh );
auto l = form1(_test=Xh);
Table 12. Table of parameters

Name

Parameter

Description

Status

test

function space e.g. Xh

define test function space

Required

Here are some examples taken from the Feel++ tutorial.

// right hand side
auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh), _expr=id(v));

From myadvection.cpp

// right hand side
auto l = form1( _test=Xh );
l+= integrate( _range=elements( mesh ), _expr=f*id( v ) );
The operators += and = are supported by linear and bilinear forms.
auto a1 = form2(_test=Xh,_trial=Xh);
auto a2 = form2(_test=Xh,_trial=Xh);
// operations on a2 ...
// check that they have the same type and
// copy matrix associated to a2 in a1
a1 = a2;
Using form2

To define a bilinear form \(\)a: X_h \times X_h \rightarrow \mathbb{R}\(\), for example \(\)a(u,v)=\int_\Omega uv\(\)

Building form2

The free-function form2 allows you to simply define such a bilinear form using the Feel++ language:

// define function space
auto Xh = ...;
// define a : Xh x Xh -> R
auto a = form2(_trial=Xh, _test=Xh );
// a(u,v) = \int_\Omega u v
a = integrate(_range=elements(mesh), _expr=idt(u)*id(v));
Table 13. Table of parameters

Name

Parameter

Description

Status

test

function space e.g. Xh

define test function space

Required

trial

function space e.g. Xh

define trial function space

Optional

Here are some examples taken from the Feel++ tutorial

From mylaplacian.cpp

// left hand side
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=elements(mesh),
              _expr=gradt(u)*trans(grad(v)) );

From mystokes.cpp:

// left hand side
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=elements(mesh),
              _expr=trace(gradt(u)*trans(grad(u))) );
a+= integrate(_range=elements(mesh),
              _expr=-div(u)*idt(p)-divt(u)*id(p));
see note above on operators += and =
Solving variational formulations

Once you created your linear and bilinear forms you can use the solve() member function of your bilinear form.

The following generic example solves: find \(\)u \in X_h \text{ such that } a(u,v)=l(v) \forall v \in X_h\(\)

Example
auto Xh = ...; // function space
auto u = Xh->element();
auto a = form2(_test=Xh, _trial=Xh);
auto l = form1(_test=Xh);

a.solve(_solution=u, _rhs=l, _rebuild=false, _name="");
Table 14. Table of arguments for solve()

Name

Parameter

Description

Status

_solution

element of domain function space

the solution

Required

_rhs

linear form

right hand side

Required

_rebuild

boolean(Default = false)

rebuild the solver components

Optional

_name

string(Default = "")

name of the associated Backend

Optional

Here are some examples from the Feel++ tutorial.

From laplacian.cpp
// solve the equation  a(u,v) = l(v)
a.solve(_rhs=l,_solution=u);
Using on for Dirichlet conditions

The function on() allows you to add Dirichlet conditions to your bilinear form before using the solve function.

The interface is as follows

Interface
on(_range=..., _rhs=..., _element=..., _expr=...);

Required Parameters:

  • _range domain concerned by this condition (see Integrals ).

  • _rhs right hand side. The linear form.

  • _element element concerned.

  • _expr the condition.

This function is used with += operator.

Here are some examples from the Feel++ tutorial.

From mylaplacian.cpp
// apply the boundary condition
a+=on(_range=boundaryfaces(mesh),
       _rhs=l,
       _element=u,
      _expr=expr(soption("functions.alpha")) );

There we add the condition: \(\)u = 0 \text{ on }\;\partial\Omega \;\(\).

From mystokes.cpp
a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
      _expr=expr<2,1,5>(u_exact,syms));

You can also apply boundary conditions per component:

Component-wise Dirichlet conditions
a+=on(_range=markedfaces(mesh,"top"),
      _element=u[Component::Y],
      _rhs=l,
      _expr=cst(0.))

The notation u[Component:Y] allows to access the Y component of u. Component::X and Component::Z are respectively the X and Z components.

14.7. Algebraic solutions

14.7.1. Definitions

Matrices

Matrix Definition A matrix is a linear transformation between finite dimensional vector spaces.

Assembling a matrix Assembling a matrix means defining its action as entries stored in a sparse or dense format. For example, in the finite element context, the storage format is sparse to take advantage of the many zero entries.

Symmetric matrix \(\)A = A^T\(\)

Definite (resp. semi-definite) positive matrix All eigenvalue are 1. \(\)>0\(\) (resp \(\)\geq 0\(\)) or 2. \(\)x^T A x > 0, \forall\ x \(\) (resp. \(\)x^T\ A\ x \geq 0\, \forall\ x\(\))

Definite (resp. semi-negative) matrix All eigenvalue are 1. \(\)<0\(\) (resp. \(\)\leq 0\(\)) or 2. \(\)x^T\ A\ x < 0\ \forall\ x\(\) (resp. \(\)x^T\ A\ x \leq 0\, \forall\ x\(\))

Indefinite matrix There exists 1. positive and negative eigenvalue (Stokes, Helmholtz) or 2. there exists \(\)x,y\(\) such that \(\)x^TAx > 0 > y^T A y\(\)

Preconditioners
Definition

Let \(\)A\(\) be a \(\)\mathbb{R}^{n\times n}\(\) matrix, \(\)x\(\) and \(\)b\(\) be \(\)\mathbb{R}^n\(\) vectors, we wish to solve \(\)A x = b.\(\)

Definition: A preconditioner \(\)\mathcal{P}\(\) is a method for constructing a matrix (just a linear function, not assembled!) \(\)P^{-1} = \mathcal{P}(A,A_p)\(\) using a matrix \(\)A\(\) and extra information \(\)A_p\(\), such that the spectrum of \(\)P^{-1}A\(\) (left preconditioning) or \(\)A P^{-1}\(\) (right preconditioning) is well-behaved. The action of preconditioning improves the conditioning of the previous linear system.

Left preconditioning: We solve for \(\) (P^{-1} A) x = P^{-1} b \(\) and we build the Krylov space \(\)\{ P^{-1} b, (P^{-1}A) P^{-1} b, (P{-1}A)2 P^{-1} b, \dots\}\(\)

Right preconditioning: We solve for \(\) (A P^{-1}) P x = b \(\) and we build the Krylov space \(\)\{ b, (P^{-1}A)b, (P{-1}A)2b, \dotsc \}\(\)

Note that the product \(\)P^{-1}A\(\) or \(\)A P^{-1}\(\) is never assembled.

Properties

Let us now describe some properties of preconditioners

  • \(\)P^{-1}\(\) is dense, \(\)P\(\) is often not available and is not needed

  • \(\)A\(\) is rarely used by \(\)\mathcal{P}\(\), but \(\)A_p = A\(\) is common

  • \(\)A_p\(\) is often a sparse matrix, the \e preconditioning \e matrix

Here are some numerical methods to solve the system \(\)A x = b\(\)

  • Matrix-based: Jacobi, Gauss-Seidel, SOR, ILU(k), LU

  • Parallel: Block-Jacobi, Schwarz, Multigrid, FETI-DP, BDDC

  • Indefinite: Schur-complement, Domain Decomposition, Multigrid

Preconditioner strategies
Relaxation

Split into lower, diagonal, upper parts: \(\) A = L + D + U \(\).

Jacobi

Cheapest preconditioner: \(\)P{-1}=D{-1}\(\).

# sequential
pc-type=jacobi
# parallel
pc-type=block_jacobi
Successive over-relaxation (SOR)

\(\) \left(L + \frac 1 \omega D\right) x_{n+1} = \left[\left(\frac 1\omega-1\right)D - U\right] x_n + \omega b \\ P^{-1} = \text{\(k\) iterations starting with \(x_0=0\)}\\ \(\)

  • Implemented as a sweep.

  • \(\)\omega = 1\(\) corresponds to Gauss-Seidel.

  • Very effective at removing high-frequency components of residual.

# sequential
pc-type=sor
Factorization

Two phases

  • symbolic factorization: find where fill occurs, only uses sparsity pattern.

  • numeric factorization: compute factors.

LU decomposition
  • preconditioner.

  • Expensive, for \(\)m\times m\(\) sparse matrix with bandwidth \(\)b\(\), traditionally requires \(\)\mathcal{O}(mb^2)\(\) time and \(\)\mathcal{O}(mb)\(\) space.

    • Bandwidth scales as \(\)m^{\frac{d-1}{d}}\(\) in \(\)d\(\)-dimensions.

    • Optimal in 2D: \(\)\mathcal{O}(m \cdot \log m)\(\) space, \(\)\mathcal{O}(m^{3/2})\(\) time.

    • Optimal in 3D: \(\)\mathcal{O}(m^{4/3})\(\) space, \(\)\mathcal{O}(m^2)\(\) time.

  • Symbolic factorization is problematic in parallel.

Incomplete LU
  • Allow a limited number of levels of fill: ILU(\(\)k\(\)).

  • Only allow fill for entries that exceed threshold: ILUT.

  • Usually poor scaling in parallel.

  • No guarantees.

1-level Domain decomposition
Domain size pass:[\(\)]Lpass:[\(\)], subdomain size pass:[\(\)]Hpass:[\(\)], element size pass:[\(\)]hpass:[\(\)]
  • Overlapping/Schwarz

    • Solve Dirichlet problems on overlapping subdomains.

    • No overlap: \(\)\textit{its} \in \mathcal{O}\big( \frac{L}{\sqrt{Hh}} \big)\(\).

    • Overlap \(\)\delta\(\): \(\)\textit{its} \in \big( \frac L {\sqrt{H\delta}} \big)\(\).

pc-type=gasm # has a coarse grid preconditioner
pc-type=asm
  • Neumann-Neumann

    • Solve Neumann problems on non-overlapping subdomains.

    • \(\)\textit{its} \in \mathcal{O}\big( \frac{L}{H}(1+\log\frac H h) \big)\(\).

    • Tricky null space issues (floating subdomains).

    • Need subdomain matrices, not globally assembled matrix.

Notes: Multilevel variants knock off the leading \(\)\frac L H\(\).
Both overlapping and nonoverlapping with this bound.

  • BDDC and FETI-DP

    • Neumann problems on subdomains with coarse grid correction.

    • \(\)\textit{its} \in \mathcal{O}\big(1 + \log\frac H h \big)\(\).

Multigrid

Hierarchy: Interpolation and restriction operators \(\) \Pi^\uparrow : X_{\text{coarse}} \to X_{\text{fine}} \qquad \Pi^\downarrow : X_{\text{fine}} \to X_{\text{coarse}} \(\)

  • Geometric: define problem on multiple levels, use grid to compute hierarchy.

  • Algebraic: define problem only on finest level, use matrix structure to build hierarchy.

Galerkin approximation

Assemble this matrix: \(\)A_{\text{coarse}} = \Pi^\downarrow A_{\text{fine}} \Pi^\uparrow\(\)

Application of multigrid preconditioner (\(\) V \(\)-cycle)

  • Apply pre-smoother on fine level (any preconditioner).

  • Restrict residual to coarse level with \(\)\Pi^\downarrow\(\).

  • Solve on coarse level \(\)A_{\text{coarse}} x = r\(\).

  • Interpolate result back to fine level with \(\)\Pi^\uparrow\(\).

  • Apply post-smoother on fine level (any preconditioner).

Multigrid convergence properties
  • Textbook: \(\)P^{-1}A\(\) is spectrally equivalent to identity

    • Constant number of iterations to converge up to discretization error.

  • Most theory applies to SPD systems

    • variable coefficients (e.g. discontinuous): low energy interpolants.

    • mesh- and/or physics-induced anisotropy: semi-coarsening/line smoothers.

    • complex geometry: difficult to have meaningful coarse levels.

  • Deeper algorithmic difficulties

    • nonsymmetric (e.g. advection, shallow water, Euler).

    • indefinite (e.g. incompressible flow, Helmholtz).

  • Performance considerations

    • Aggressive coarsening is critical in parallel.

    • Most theory uses SOR smoothers, ILU often more robust.

    • Coarsest level usually solved semi-redundantly with direct solver.

  • Multilevel Schwarz is essentially the same with different language

    • assume strong smoothers, emphasize aggressive coarsening.

List of PETSc Preconditioners

See this PETSc page for a complete list.

Table 15. Table of Preconditioners as of PETSc 3.7

PETSc

Description

Parallel

none

No preconditioner

yes

jacobi

diagonal preconditioner

yes

bjacobi

block diagonal preconditioner

yes

sor

SOR preconditioner

yes

lu

Direct solver as preconditioner

depends on the factorization package (e.g.mumps,pastix: OK)

shell

User defined preconditioner

depends on the user preconditioner

mg

multigrid prec

yes

ilu

incomplete lu

icc

incomplete cholesky

cholesky

Cholesky factorisation

yes

asm

Additive Schwarz Method

yes

gasm

Scalable Additive Schwarz Method

yes

ksp

Krylov subspace preconditioner

yes

fieldsplit

block preconditioner framework

yes

lsc

Least Square Commutator

yes

gamg

Scalable Algebraic Multigrid

yes

hypre

Hypre framework (multigrid…​)

bddc

balancing domain decomposition by constraints preconditioner

yes

14.7.2. Principles

Feel++ abstracts the PETSc library and provides a subset (sufficient in most cases) to the PETSc features. It interfaces with the following PETSc libraries: Mat , Vec , KSP , PC , SNES.

  • Vec Vector handling library

  • Mat Matrix handling library

  • KSP Krylov SubSpace library implements various iterative solvers

  • PC Preconditioner library implements various preconditioning strategies

  • SNES Nonlinear solver library implements various nonlinear solve strategies

All linear algebra are encapsulated within backends using the command line option --backend=<backend> or config file option backend=<backend> which provide interface to several libraries

Library

Format

Backend

PETSc

sparse

petsc

Eigen

sparse

eigen

Eigen

dense

eigen_dense

The default backend is petsc.

14.7.3. Somes generic examples

The configuration files .cfg allow for a wide range of options to solve a linear or non-linear system.

We consider now the following example [import:"marker1"](../../codes/mylaplacian.cpp)

To execute this example

# sequential
./feelpp_tut_laplacian
# parallel on 4 cores
mpirun -np 4 ./feelpp_tut_laplacian

As described in section

Direct solver

cholesky and lu factorisation are available. However the parallel implementation depends on the availability of MUMPS. The configuration is very simple.

# no iterative solver
ksp-type=preonly
#
pc-type=cholesky

Using the PETSc backend allows to choose different packages to compute the factorization.

Table 16. Table of factorization package

Package

Description

Parallel

petsc

PETSc own implementation

yes

mumps

MUltifrontal Massively Parallel sparse direct Solver

yes

umfpack

Unsymmetric MultiFrontal package

no

pastix

Parallel Sparse matriX package

yes

To choose between these factorization package

# choose mumps
pc-factor-mat-solver-package=mumps
# choose umfpack (sequential)
pc-factor-mat-solver-package=umfpack

In order to perform a cholesky type of factorisation, it is required to set the underlying matrix to be SPD.

// matrix
auto A = backend->newMatrix(_test=...,_trial=...,_properties=SPD);
// bilinear form
auto a = form2( _test=..., _trial=..., _properties=SPD );
Using iterative solvers
Using CG and ICC(3)

with a relative tolerance of 1e-12:

ksp-rtol=1.e-12
ksp-type=cg
pc-type=icc
pc-factor-levels=3
Using GMRES and ILU(3)

with a relative tolerance of 1e-12 and a restart of 300:

ksp-rtol=1.e-12
ksp-type=gmres
ksp-gmres-restart=300
pc-type=ilu
pc-factor-levels=3
Using GMRES and Jacobi

With a relative tolerance of 1e-12 and a restart of 100:

ksp-rtol=1.e-12
ksp-type=gmres
ksp-gmres-restart 100
pc-type=jacobi
Monitoring linear non-linear and eigen problem solver residuals
# linear
ksp_monitor=1
# non-linear
snes-monitor=1
# eigen value problem
eps-monitor=1

14.7.4. Solving the Laplace problem

We start with the quickstart Laplacian example, recall that we wish to, given a domain \(\)\Omega\(\), find \(\)u\(\) such that

\(\) -\nabla \cdot (k \nabla u) = f \mbox{ in } \Omega \subset \mathbb{R}^{2},\\ u = g \mbox{ on } \partial \Omega \(\)

Monitoring KSP solvers
feelpp_qs_laplacian --ksp-monitor=true
Viewing KSP solvers
shell> mpirun -np 2 feelpp_qs_laplacian --ksp-monitor=1  --ksp-view=1
  0 KSP Residual norm 8.953261456448e-01
  1 KSP Residual norm 7.204431786960e-16
KSP Object: 2 MPI processes
  type: gmres
    GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
     Orthogonalization with no iterative refinement
    GMRES: happy breakdown tolerance 1e-30
  maximum iterations=1000
  tolerances:  relative=1e-13, absolute=1e-50, divergence=100000
  left preconditioning
  using nonzero initial guess
  using PRECONDITIONED norm type for convergence test
PC Object: 2 MPI processes
  type: shell
    Shell:
  linear system matrix = precond matrix:
  Matrix Object:   2 MPI processes
    type: mpiaij
    rows=525, cols=525
    total: nonzeros=5727, allocated nonzeros=5727
    total number of mallocs used during MatSetValues calls =0
      not using I-node (on process 0) routines
Solvers and preconditioners

You can now change the Krylov subspace solver using the --ksp-type option and the preconditioner using --pc-ptype option.

For example,

  • to solve use the conjugate gradient,cg, solver and the default preconditioner use the following

./feelpp_qs_laplacian --ksp-type=cg --ksp-view=1 --ksp-monitor=1
  • to solve using the algebraic multigrid preconditioner, gamg, with cg as a solver use the following

./feelpp_qs_laplacian --ksp-type=cg --ksp-view=1 --ksp-monitor=1 --pc-type=gamg

14.7.5. Block factorisation

Stokes

We now turn to the quickstart Stokes example, recall that we wish to, given a domain \(\)\Omega\(\), find \(\)(\mathbf{u},p) \(\) such that

\[ -\Delta \mathbf{u} + \nabla p = \mathbf{ f} \mbox{ in } \Omega,\\ \nabla \cdot \mathbf{u} = 0 \mbox{ in } \Omega,\\ \mathbf{u} = \mathbf{g} \mbox{ on } \partial \Omega\]

This problem is indefinite. Possible solution strategies are

  • Uzawa,

  • penalty(techniques from optimisation),

  • augmented lagrangian approach (Glowinski,Le Tallec)

Note that The Inf-sup condition must be satisfied. In particular for a multigrid strategy, the smoother needs to preserve it.

14.7.6. General approach for saddle point problems

The Krylov subspace solvers for indefinite problems are MINRES, GMRES. As to preconditioning, we look first at the saddle point matrix \(\)M\(\) and its block factorization \(\)M = LDL^T\(\), indeed we have :

\[M = \begin{pmatrix} A & B \\ B^T & 0 \end{pmatrix} = \begin{pmatrix} I & 0\\ B^T C & I \end{pmatrix} \begin{pmatrix} A & 0\\ 0 & - B^T A^{-1} B \end{pmatrix} \begin{pmatrix} I & A^{-1} B\\ 0 & I \end{pmatrix}\]
  • Elman, Silvester and Wathen propose 3 preconditioners:

\[P_1 = \begin{pmatrix} \tilde{A}^{-1} & B\\ B^T & 0 \end{pmatrix}, \quad P_2 = \begin{pmatrix} \tilde{A}^{-1} & 0\\ 0 & \tilde{S} \end{pmatrix},\quad P_3 = \begin{pmatrix} \tilde{A}^{-1} & B\\ 0 & \tilde{S} \end{pmatrix}\]

where \(\tilde{S} \approx S^{-1} = B^T A^{-1} B\) and \(\tilde{A}^{-1} \approx A^{-1}\)

14.7.7. Preconditioner strategies

Relaxation

Split into lower, diagonal, upper parts: \(A = L + D + U\).

Jacobi

Cheapest preconditioner: \(P^{-1}=D^{-1}\).

# sequential
pc-type=jacobi
# parallel
pc-type=block_jacobi
Successive over-relaxation (SOR)
\[\left(L + \frac 1 \omega D\right) x_{n+1} = \left[\left(\frac 1\omega-1\right)D - U\right] x_n + \omega b \\ P^{-1} = \text{pass:[\(k\)] iterations starting with pass:[\(x_0=0\)]}\]
  • Implemented as a sweep.

  • \(\omega = 1\) corresponds to Gauss-Seidel.

  • Very effective at removing high-frequency components of residual.

# sequential
pc-type=sor
Factorization

Two phases

  • symbolic factorization: find where fill occurs, only uses sparsity pattern.

  • numeric factorization: compute factors.

LU decomposition
  • preconditioner.

  • Expensive, for \(m\times m\) sparse matrix with bandwidth \(b\), traditionally requires \(\mathcal{O}(mb^2)\) time and \(\mathcal{O}(mb)\) space.

    • Bandwidth scales as \(m^{\frac{d-1}{d}}\) in \(\)d\(\)-dimensions.

    • Optimal in 2D: \(\mathcal{O}(m \cdot \log m)\) space, \(\mathcal{O}(m^{3/2})\) time.

    • Optimal in 3D: \(\mathcal{O}(m^{4/3})\) space, \(\mathcal{O}(m^2)\) time.

  • Symbolic factorization is problematic in parallel.

Incomplete LU
  • Allow a limited number of levels of fill: ILU(\(k\)).

  • Only allow fill for entries that exceed threshold: ILUT.

  • Usually poor scaling in parallel.

  • No guarantees.

1-level Domain decomposition
Domain size pass:[\(\)]Lpass:[\(\)], subdomain size pass:[\(\)]Hpass:[\(\)], element size pass:[\(\)]hpass:[\(\)]
  • Overlapping/Schwarz

    • Solve Dirichlet problems on overlapping subdomains.

    • No overlap: \(\textit{its} \in \mathcal{O}\big( \frac{L}{\sqrt{Hh}} \big)\).

    • Overlap \(\delta: \textit{its} \in \big( \frac L {\sqrt{H\delta}} \big)\).

  • Neumann-Neumann

    • Solve Neumann problems on non-overlapping subdomains.

    • \(\textit{its} \in \mathcal{O}\big( \frac{L}{H}(1+\log\frac H h) \big)\).

    • Tricky null space issues (floating subdomains).

    • Need subdomain matrices, not globally assembled matrix.

Multilevel variants knock off the leading \(\frac L H\).
Both overlapping and nonoverlapping with this bound.
  • BDDC and FETI-DP

    • Neumann problems on subdomains with coarse grid correction.

    • \(\textit{its} \in \mathcal{O}\big(1 + \log\frac H h \big)\).

Multigrid
Introduction

Hierarchy: Interpolation and restriction operators \(\Pi^\uparrow : X_{\text{coarse}} \to X_{\text{fine}} \qquad \Pi^\downarrow : X_{\text{fine}} \to X_{\text{coarse}} \)

  • Geometric: define problem on multiple levels, use grid to compute hierarchy.

  • Algebraic: define problem only on finest level, use matrix structure to build hierarchy.

Galerkin approximation

Assemble this matrix: \(A_{\text{coarse}} = \Pi^\downarrow A_{\text{fine}} \Pi^\uparrow\)

Application of multigrid preconditioner (\(V\)-cycle)

  • Apply pre-smoother on fine level (any preconditioner).

  • Restrict residual to coarse level with \(\Pi^\downarrow\).

  • Solve on coarse level \(A_{\text{coarse}} x = r\).

  • Interpolate result back to fine level with \(\)\Pi^\uparrow\(\).

  • Apply post-smoother on fine level (any preconditioner).

Multigrid convergence properties
  • Textbook: \(P^{-1}A\) is spectrally equivalent to identity

    • Constant number of iterations to converge up to discretization error.

  • Most theory applies to SPD systems

    • variable coefficients (e.g. discontinuous): low energy interpolants.

    • mesh- and/or physics-induced anisotropy: semi-coarsening/line smoothers.

    • complex geometry: difficult to have meaningful coarse levels.

  • Deeper algorithmic difficulties

    • nonsymmetric (e.g. advection, shallow water, Euler).

    • indefinite (e.g. incompressible flow, Helmholtz).

  • Performance considerations

    • Aggressive coarsening is critical in parallel.

    • Most theory uses SOR smoothers, ILU often more robust.

    • Coarsest level usually solved semi-redundantly with direct solver.

  • Multilevel Schwarz is essentially the same with different language

    • assume strong smoothers, emphasize aggressive coarsening.

List of PETSc Preconditioners

See this PETSc page for a complete list.

Table 17. Table of Preconditioners as of PETSc 3.7

PETSc

Description

Parallel

none

No preconditioner

yes

jacobi

diagonal preconditioner

yes

bjacobi

block diagonal preconditioner

yes

sor

SOR preconditioner

yes

lu

Direct solver as preconditioner

depends on the factorization package (e.g.mumps,pastix: OK)

shell

User defined preconditioner

depends on the user preconditioner

mg

multigrid prec

yes

ilu

incomplete lu

icc

incomplete cholesky

cholesky

Cholesky factorisation

yes

asm

Additive Schwarz Method

yes

gasm

Scalable Additive Schwarz Method

yes

ksp

Krylov subspace preconditioner

yes

fieldsplit

block preconditioner framework

yes

lsc

Least Square Commutator

yes

gamg

Scalable Algebraic Multigrid

yes

hypre

Hypre framework (multigrid…​)

bddc

balancing domain decomposition by constraints preconditioner

yes

14.7.8. Algebra Backends

Non-Linear algebra backends are crucial components of Feel++. They provide a uniform interface between Feel++ data structures and underlying the linear algebra libraries used by Feel++.

Libraries

Feel++ interfaces the following libraries:

  • PETSc : Portable, Extensible Toolkit for Scientific Computation

  • SLEPc : Eigen value solver framework based on PETSc

  • Eigen3

Backend

Backend is a template class parametrized by the numerical type providing access to

  • vector

  • matrix : dense and sparse

  • algorithms : solvers, preconditioners, …​

PETSc provides sequential and parallel data structures whereas Eigen3 is only sequential.

To create a Backend, use the free function backend(…​) which has the following interface:

backend(_name="name_of_backend",
        _rebuild=... /// true|false,
        _kind=..., // type of backend,
        _worldcomm=... // communicator
        )

All these parameters are optional which means that the simplest call reads for example:

auto b = backend();

They take default values as described in the following table:

parameter

type

default value

_name

string

"" (empty string )

_rebuild

boolean

false

_kind

string

"petsc"

_worldcomm

WorldComm

Environment::worldComm()

_name

Backends are store in a Backend factory and handled by a manager that allows to keep track of allocated backends. They a registered with respect to their name and kind. The default name value is en empty string ("") which names the default Backend. The _name parameter is important because it provides also the name for the command line/config file option section associated to the associated Backend.

Only the command line/config file options for the default backend are registered. Developers have to register the option for each new Backend they define otherwise failures at run-time are to be expected whenever a Backend command line option is accessed.

Consider that you create a Backend name projection in your code like this

auto b = backend(_name="projection");

to register the command line options of this Backend

Environment env( _argc=argc, _argv=argv,
                 _desc=backend_options("projection") );
_kind

Feel++ supports three kind of Backends:

  • petsc : PETSC Backend

  • eigen_dense

  • eigen_sparse

SLEPc uses the PETSc Backend since it is based on PETSc.

The kind of Backend can be changed from the command line or configuration file thanks to the "backend" option.

auto b = backend(_name="name",
                 _kind=soption(_prefix="name",_name="backend"))

and in the config file

[name]
backend=petsc
backend=eigen_sparse
_rebuild

If you want to reuse a Backend and not allocate a new one plus add the corresponding option to the command line/configuration file, you can rebuild the Backend in order to delete the data structures already associated to this Backend and in particular the preconditioner. It is mandatory to do that when you solve say a linear system first with dimensions \(m\times m\) and then solve another one with different dimension \(n \times n\) because in that case the Backend will throw an error saying that the dimensions are incompatible. To avoid that you have to rebuild the Backend.

auto b = backend(_name="mybackend");
// solve A x = f
b->solve(_solution=x,_matrix=A,_rhs=f);
// rebuild: clean up the internal Backend data structure
b = backend(_name="mybackend",_rebuild=true);
// solve A1 x1 = f1
b->solve(_solution=x1,_matrix=A1,_rhs=f1);
Although this feature might appear useful, you have to make sure that the solving strategy applies to all problems because you won’t be able to customize the solver/preconditioner for each problem. If you have different problems to solve and need to have custom solver/preconditioner it would be best to have different Backends.
_worldComm

One of the strength of Feel++ is to be able to change the communicator and in the case of Feel++ the WorldComm. This allows for example to

  • solve sequential problems

  • solve a problem on a subset of MPI processes

For example passing a sequential WorldComm to backend() would then in the subsequent use of the Backend generate sequential data structures (e.g. IndexSet, Vector and Matrix) and algorithms (e.g. Krylov Solvers).

 // create a sequential Backend
 auto b = backend(_name="seq",
                  _worldComm=Environment::worldCommSeq());
auto A = b->newMatrix(); // sequential Matrix
auto f = b->newVector(); // sequential Vector

Info The default WorldComm is provided by Environment::worldComm() and it conresponds to the default MPI communicator MPI_COMM_WORLD.

14.7.9. Eigen Problem

To solve standard and generalized eigenvalue problems, Feel++ interfaces SLEPc. SLEPc is a library which extends PETSc to provide the functionality necessary for the solution of eigenvalue problems. It comes with many strategies for both standard and generalized problems, Hermitian or not.

We want to find \((\lambda_i,x_i)\) such that \(Ax = \lambda x\). To do that, most eigensolvers project the problem onto a low-dimensional subspace, this is called a Rayleigh-Ritz projection. + Let \(V_j=[v_1,v_2,...,v_j\)] be an orthogonal basis of this subspace, then the projected problem reads:

Find \((\theta_i,s_i)\) for \(i=1,\ldots,j\) such that \(B_j s_i=\theta_i s_i\) where \(B_j=V_j^T A V_j\).

Then the approximate eigenpairs \((\lambda_i,x_i)\) of the original problem are obtained as: \(\lambda_i=\theta_i\) and \(x_i=V_j s_i\).

The eigensolvers differ from each other in the way the subspace is built.

Code

In Feel++, there is two functions that can be used to solve this type of problems, eigs and veigs.

Here is an example of how to use veigs.

auto Vh = Pch<Order>( mesh );
auto a = form2( _test=Vh, _trial=Vh );
// fill a
auto b = form2( _test=Vh, _trial=Vh );
// fill b
auto eigenmodes = veigs( _formA=a, _formB=b );

where eigenmodes is a std::vector<std::pair<value_type, element_type> > with value_type the type of the eigenvalue, and element_type the type of the eigenvector, a function of the space Vh.

The eigs function does not take the bilinear forms but two matrices. Also, the solver used, the type of the problem, the position of the spectrum and the spectral transformation are not read from the options.

auto Vh = Pch<Order>( mesh );
auto a = form2( _test=Vh, _trial=Vh );
// fill a
auto matA = a.matrixPtr();
auto b = form2( _test=Vh, _trial=Vh );
// fill b
auto matB = b.matrixPtr();
auto eigenmodes = eigs( _matrixA=aHat,
                        _matrixB=bHat,
                        _solver=(EigenSolverType)EigenMap[soption("solvereigen.solver")],
                        _problem=(EigenProblemType)EigenMap[soption("solvereigen.problem")],
                        _transform=(SpectralTransformType)EigenMap[soption("solvereigen.transform")],
                        _spectrum=(PositionOfSpectrum)EigenMap[soption("solvereigen.spectrum")]
                         );
auto femodes = std::vector<decltype(Vh->element())>( eigenmodes.size(), Vh->element() );
int i = 0;
for( auto const& mode : modes )
    femodes[i++] = *mode.second.get<2>();

where eigenmodes is a std::map<real_type, eigenmodes_type> with real_type of the magnitude of the eigenvalue. And eigenmodes_type is a boost::tuple<real_type, real_type, vector_ptrtype> with the first real_type representing the real part of the eigenvalue, the second real_type the imaginary part and the vector_ptrtype is a vector but not an element of a functionspace.

The two functions take a parameter _nev that tel how many eigenpair to compute. This can be set from the command line option --solvereigen.nev. + Another important parameter is _ncv which is the size of the subspace, j above. This parameter should always be greater than nev. SLEPc recommends to set it to max(nev+15, 2*nev). This can be set from the command line option --solvereigen.ncv.

Problem type

The standard formulation reads :

Find \(\lambda\in \mathbb{R}\) such that \(Ax = \lambda x\)

where \(\lambda\) is an eigenvalue and \(x\) an eigenvector.

But in the case of the finite element method, we will deal with the generalized form :

Find \(\lambda\in\mathbb{R}\) such that \(Ax = \lambda Bx\)

A standard problem is Hermitian if the matrix \(\)A\(\) is Hermitian (\(A=A^*\)). + A generalized problem is Hermitian if the matrices \(A\) and \(B\) are Hermitian and if \(B\) is positive definite. + If the problem is Hermitian, then the eigenvalues are real. A special case of the generalized problem is when the matrices are not Hermitian but \(B\) is positive definite.

The type of the problem can be specified using the EigenProblemType, or at run time with the command line option --solvereigen.problem and the following value :

Table 18. Table of problem type
Problem type EigenProblemType command line key

Standard Hermitian

HEP

"hep"

Standard non-Hermitian

NHEP

"nhep"

Generalized Hermitian

GHEP

"ghep"

Generalized non-Hermitian

GNHEP

"gnhep"

Positive definite Generalized non-Hermitian

PGNHEP

"pgnhep"

Position of spectrum

You can choose which eigenpairs will be computed. The user can set it programmatically with PositionOfSpectrum or at run time with the command line option --solvereigen.spectrum and the following value :

Table 19. Table of position of spectrum
Position of spectrum PositionOfSpectrum command line key

Largest magnitude

LARGEST_MAGNITUDE

"largest_magnitude"

Smallest magnitude

SMALLEST_MAGNITUDE

"smallest_magnitude"

Largest real

LARGEST_REAL

"largest_real"

Smallest real

SMALLEST_REAL

"smallest_real"

Largest imaginary

LARGEST_IMAGINARY

"largest_imaginary"

Smallest imaginary

SMALLEST_IMAGINARY

"smallest_imaginary"

Spectral transformation

It is observed that the algorithms used to solve the eigenvalue problems find solutions at the extremities of the spectrum. To improve the convergence, one need to compute the eigenpairs of a transformed operator. Those spectral transformations allow to compute solutions that are not on the boundary of the spectrum.

There are 3 types of spectral transformation:

Shift

\(A-\sigma I\) or \(B^{-1}A-\sigma I\)

Shift and invert

\((A-\sigma I)^{-1}\) or \((A-\sigma B)^{-1}B\)

Cayley

\((A-\sigma I)^{-1}(A+\nu I)\) or \((A-\sigma B)^{-1}(A+\nu B)\)

By default, shift and invert is used. You can change it with --solvereigen.transform.

Table 20. Table of spectral transformation
Spectral transformation SpectralTransformationType command line key

Shift

SHIFT

shift

Shift and invert

SINVERT

shift_invert

Cayley

CAYLEY

cayley

Eigensolvers

The details of the implementation of the different solvers can be found in the SLEPc Technical Reports.

The default solver is Krylov-Schur, but can be modified using EigenSolverType or the option --solvereigen.solver.

Table 21. Table of eigensolver
Solver EigenSolverType command line key

Power

POWER

power

Lapack

LAPACK

lapack

Subspace

SUBSPACE

subspace

Arnoldi

Arnoldi

arnoldi

Lanczos

LANCZOS

lanczos

Krylov-Schur

KRYLOVSCHUR

krylovschur

Arpack

ARPACK

arpack

Be careful that all solvers can not compute all the problem types and positions of the spectrum. The possibilities are summarize in the following table.

Table 22. Supported problem type for the eigensolvers
Solver Position of spectrum Problem type

Power

Largest magnitude

any

Lapack

any

any

Subspace

Largest magnitude

any

Arnoldi

any

any

Lanczos

any

standard and generalized Hermitian

Krylov-Schur

any

any

Arpack

any

any

Special cases of spectrum
Computing a large portion of the spectrum

In the case where you want compute a large number of eigenpairs, the rule for ncv implies a huge amount of memory to be used. To improve the performance, you can set the mpd parameter, which will limit the dimension of the projected problem.

You can set it via the command line with --solvereigen.mpd <mpd>.

Computing all the eigenpairs in an interval

If you want to compute all the eigenpairs in a given interval, you need to use the option --solvereigen.interval-a to set the beginning of the interval and --solvereigen.interval-b to set the end.

In this case, be aware that the problem need to be generalized and hermitian. The solver will be set to Krylov-Schur and the transformation to shift and invert. Beside, you’ll need to use a linear solver that will compute the inertia of the matrix, this is set to Cholesky, with mumps if you can use it. + For now, this method is only implemented in the eigs function.

Appendix A: Feel++ File Formats

Feel++ Formats

For performance reasons and allow fast checkpoint restart of simulations, we have develop our own mesh and data file format in parallel.

Format

Description

Mode

Type

json+hdf5

Feel++ parallel file format

Read/Write

Metadata & Binary

The format is decomposed into two files : (i) a json file (.json file extension) which contain some metadata information on the mesh and (ii) a hdf5 file (.h5 file extension) which contains the mesh data structure.

Pre-Processing formats

Feel++ supports various file formats that can be used as input mesh file formats.

Format

Description

Mode

Type

acumesh

Acusim(ALTAIR) mesh file format

Read

Ascii

gmsh

Gmsh mesh file format

Read/Write

Ascii/Binary

json+hdf5

Feel++ parallel file format

Read/Write

Metadata & Binary

med

MED(Salome) mesh file format

Read/Write

Ascii/Binary

mesh

MEDIT(INRIA) mesh file format

Read/Write

Ascii

Post Processing formats

Feel++ supports various file formats that can be used as output mesh and data file formats for post-processing.

Format

Description

Mode

Type

gmsh

Gmsh mesh file format

Read/Write

Ascii/Binary

ensightgold

Ensight Gold case format

Write

Binary

h3d

H3D file format

Read

Database

xdmf

XML/HDF5 file format

Write

VTK

The H3D file format requires that you have the Altair Hypermesh software installed.

15. Ressources

15.1. Licenses

Copyright © 2010-2017 by Feel++ Consortium
Copyright © 2005-2015 by Université Joseph Fourier (Grenoble, France)
Copyright © 2005-2015 by University of Coimbra (Portugal)
Copyright © 2011-2015 by Université de Strasbourg (France)
Copyright © 2011-2015 by CNRS (France)
Copyright © 2005-2006 by Ecole Polytechnique Fédérale de Lausanne (EPFL, Switzerland)

Free use of this software is granted under the terms of the L License.

See the LICENSE file for details

This book is part of Feel++ and is licensed under cc LGPL a.

15.2. Authors

There are many other contributors.

Feel++ is currently managed by Christophe Prud’homme, Professor in applied mathematic and scientific computing at the University of Strasbourg, France.

15.3. Funding

Feel++ has been funded by various sources and especially

logo anr logo amies logo irmia logo prace

15.3.1. Current funding

EU E-INFRA H2020
ANR projects
PlasticOmnium
  • Contract (2016-2017)

Holo3
  • Contract (2016-2017)

AMIES
  • PEPS Holo3

  • PEPS Solodem

  • PEPS NS2++

IRMIA
  • Hifimagnet (2012-2018)

  • 4fastsim (2016-2017)

15.3.2. Past funding

ANR
  • HAMM - (Cosinus call - 2010-2014)

  • OPUS - (TLOG call - 2008-2011)

  • Funding for Cemosis

FRAE
  • RB4FASTSIM - 2010-2014

PRACE projects
  • HP-FEEL++ 2015-2016

  • HP-FEEL++ 2013-2014

  • HP-PDE{1,2} 2012-2014

Rhônes-Alpes region
  • cluster ISLE [fn:2] and the project CHPID (2009-2011)

15.4. Contributors

Feel++ benefits from the many discussions and close research collaborations with the following persons: Mourad Ismail, Zakaria Belhachmi, Silvia Bertoluzza, Micol Pennacchio, Marcela Szopos, Giovanna Guidoboni, Riccardo Sacco, Gonçalo Pena.

Finally Feel++ also benefits from discussions within collaborative projects with many people (in no particular order):

Yannick Hoarau, Philippe Gilotte, Benjamin Surowiec, Yoann Eulalie, Stephie Edwige, Marion Spreng, Benjamin Vanthong, Thomas Lantz, Mamadou Camara, Camille Boulard, Pierre Gerhard, Frédéric Hecht, Michel Fouquembergh, Denis Barbier, Jean-Marc Gratien, Daniele Di Pietro.

15.5. Consortium

Feel++ was initially developed at École Polytechnique Fédérale de Lausanne(Suisse) and is now a joint effort between Université de Strasbourg, Université Grenoble-Alpes, CNRS, LNCMI and Cemosis.

logo cemosis logo uga logo cnrs logo imati logo uds

Glossary

16. Glossary

boundaryelements

Free-function to apply to a mesh to retrieve the iterators over elements touching the boundary of the mesh stored on the current processor with an face, edge or point.

boundaryfaces

Free-function to apply to a mesh to retrieve the iterators over boundary faces of the mesh stored on the current processor.

Cmake

The tool that configures Feel++ build environment and generate Makefiles by default.

edges

Free-function to apply to a mesh to retrieve the iterators over the edges of the mesh stored on the current processor

Eigen3

Eigen is a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms.

elements

Free-function to apply to a mesh to retrieve the iterators over the elements of the mesh stored on the current processor

faces

Free-function to apply to a mesh to retrieve the iterators over the faces of the mesh stored on the current processor

globalRank

MPI global rank of a data structure

integrate

Free-function to define integral expressions entering the definition of integrals, linear and bi-linear forms.

internalelements

Free-function to apply to a mesh to retrieve the iterators over elements which are not touching with a point, edge or face the boundary of the mesh stored on the current processor

Make

A tool that builds Feel++ code from Makefiles generated by Cmake.

marked2elements

Free-function to apply to a mesh to retrieve the iterators over marked elements (by a string or an integer id) with marker2 of the mesh stored on the current processor

marked3elements

Free-function to apply to a mesh to retrieve the iterators over marked elements (by a string or an integer id) with marker3 of the mesh stored on the current processor

markededges

Free-function to apply to a mesh to retrieve the iterators over marked edges (by a string or an integer id) of the mesh stored on the current processor

markedelements

Free-function to apply to a mesh to retrieve the iterators over marked elements (by a string or an integer id) of the mesh stored on the current processor

markedfaces

Free-function to apply to a mesh to retrieve the iterators over marked faces (by a string or an integer id) of the mesh stored on the current processor

marker

Marker for mesh element, faces, edges or point. Element marker are often associated to material properties

marker2

Marker for mesh element, faces, edges or point. It is used for example to iterate over element thanks to a particular piecewise constant field

marker3

Marker for mesh element, faces, edges or point. It is used for example to iterate over element thanks to a particular piecewise constant field

mean

Free-function to compute the average value of a function.

MUMPS

A parallel sparse direct solvers

normH1

Free-function to compute the \(H^1\) norm of an expression

normL2

Free-function to compute the \(L^2\) norm of an expression

normLinf

Free-function to compute the \(L^{\infty}\) norm of an expression

Pastix

PaStiX (Parallel Sparse matriX package) is a scientific library that provides a high performance parallel solver for very large sparse linear systems based on direct methods. Numerical algorithms are implemented in single or double precision (real or complex) using LLt, LDLt and LU with static pivoting (for non symmetric matrices having a symmetric pattern). This solver provides also an adaptive blockwise iLU(k) factorization that can be used as a parallel preconditioner using approximated supernodes to build a coarser block structure of the incomplete factors. See http://pastix.gforge.inria.fr/.

PETSc

A library for High Performance Computing providing parallel data structures and numerical methods linear and non-linear algebraic problems arising for example PDE discretisation. PETSc is the main solver strategy provider for FEEL++.

project

Free-function to project an expression \(e\) over a nodal function space \(X_h\). It would typically return the interpolant \(\Pi_h e \in X_h\) of the expression in the function space.

rank

MPI local rank of a data structure

SLEPc

A library based on PETSc providing a framework to solve eigenvalue problems.

SPD

Symmetric Positive Definite

UMFPACK

UMFPACK /ˈʌmfpæk/ is a set of routines for solving sparse linear systems of the form Ax=b, using the Unsymmetric MultiFrontal method (Matrix A is not required to be symmetric) [source: https://en.wikipedia.org/wiki/UMFPACK]