##### https://hal.archives-ouvertes.fr/hal-03517921

README.md

```
## README
#### Description
ParaSkel++ is a C++ platform for the high-performance, arbitrary-order, 2/3D numerical approximation of PDEs on general polytopal meshes using skeletal Galerkin methods.
Skeletal Galerkin methods are a vast family of numerical methods for the approximation of PDE-based models that satisfy the following two building principles:
* the degrees of freedom (DOF) of the method split into (i) skeleton DOF, attached to the geometric entities (vertices, edges, faces) composing the mesh skeleton and common to all cells sharing the geometric entity in question, which prescribe the conformity properties of the underlying discrete functional space, and (ii) bulk DOF (if need be), attached to the interior of the cells, which play no role in the prescription of the conformity properties of the underlying discrete functional space;
* the global discrete bilinear form of the problem (potentially after linearization, if the problem is nonlinear) writes as the sum over the mesh cells of cell-wise (referred to as local) bilinear contributions.
The very structure underpinning skeletal methods grants them the property of being amenable to static condensation, i.e. locally to each cell, bulk DOF can be eliminated in terms of the local skeleton DOF. The final global system to solve thus writes in terms of the skeleton DOF only.
Examples of skeletal methods include standard FE methods and virtual-like Galerkin methods (VEM, HHO, HDG...). Notice that (plain vanilla) DG methods do not enter the skeletal framework.
ParaSkel++ offers a high-performance factorized C++ architecture for the implementation of arbitrary-order skeletal methods on general 2/3D polytopal meshes.
The ParaSkel++ code is built around the core class FiniteElement, which embeds all what concerns the local approximation space (on a given cell).
To describe the local space, one starts from the postulate that the DOF can be attached to precisely (up to) 3 different types of geometric entities in 2D (vertex, face, cell), and 4 different types
in 3D (vertex, edge, face, cell).
The fact that "edges" are termed "faces" in 2D enables to perform a unified 2/3D implementation of hybrid methods (like HHO or HDG).
The global approximation space (class DiscreteSpace) is then essentially built as a list of nbr_cells objects FiniteElement.
Locally, a complete latitude is left on the shape of the cell, the nature of the local approximation space, the polynomial degree(s), etc.
At the local level, static condensation is performed by means of a purely algebraic Schur complement routine, that is common to all methods.
The global assembly is performed in the class MatrixPattern independently of the method, taking advantage of the generic description (in 3 or 4 types of DOF) of the local spaces.
The local_to_global map (that maps, on every cell, the local index of a skeleton DOF to its index in the final global system), which is used for assembly, is built by looping over the mesh cells, and by indexing one by one the DOF groups attached to the geometric entities composing the cell boundary, merging duplicates.
The global assembly is then performed by sum over the mesh cells of the local (cell-wise) matrix and vector contributions.
Adding a new method to the code is thus essentially just a matter of implementing the local traits.
Up to now (v1), a sequential version of the following skeletal methods is available. The order of the method is denoted by k >= 1.
For the Poisson problem:
* Pk Lagrange (on simplices)
* P1 Crouzeix-Raviart (on simplices)
* HHO(b,k-1)
* VEM(b,k), and enhanced VEM(b,k) in 3D
For the Stokes problem:
* Taylor-Hood Pk - P(k-1) for k > 1 (on simplices)
* P1 Crouzeix-Raviart - P0 (on simplices)
* HHO(b,k-1) - discontinuous P(k-1)
Eventually, the ParaSkel++ platform is expected to possess five main assets: (i) a unified 2/3D implementation, (ii) the native support of any type of DOF (vertex-, edge-, face-, and cell-based), (iii) an ultra-factorized architecture (with common-to-all-methods local elimination and global assembly steps), (iv) the use of efficient quadrature formulas on general polytopal cells (without the need for subtessellation), and (v) the embedding of parallelism facilities.
#### Installation
The following softwares/libraries are required for compilation:
* C++(>=17) compiler
* CMake (version >= 3.5)
* Eigen (version >= 3.3.4)
To compile the code, you may need to configure the file CMakeLists.txt, in particular with the path to the Eigen library.
When done, you should create a build directory. The compilation will create a bin directory.
mkdir build <br>
cd build <br>
cmake .. <br>
make
There are three compiling modes (Debug, RelWithDebInfo, and Release), which can be selected with "ccmake ..".
#### Getting started
The dimension (2 or 3) and the skeletal method must be specified (in the main file) before the compilation.
To run the code (from the build directory), you need at least to specify a mesh file with the flag -m:
* ./bin/exec -m input_mesh_file_full_name
The following mesh formats are supported (specifying the extension when running is mandatory).
In 2D: the Gmsh format "*.msh", or the FVCA5 formats "*.typ1" and "*.typ2".
In 3D: the Gmsh format "*.msh", or the FVCA6 format "*.typ6".
A few 2D and 3D meshes are provided in the Meshes directory.
Beware of the fact that some methods (see above) only support simplicial meshes.
To change the order k of the method, use the option -s:
* ./bin/exec -m input_mesh_file_full_name -s val
Note that for the HHO method, the order of the method is k = s+1.
For all other methods, the order is k = s.
You will find some help about the different options of the code by running it with the flag -h.
Available options:
* -m input_mesh_file_full_name: mandatory
* -s skeleton_polynomial_degree: for all methods (except P1 CR)
* -b bulk_polynomial_degree: for methods HHO(b,s) and (enhanced) VEM(b,s)
* -v (optional: output_vtk_file_full_name, default value exists): creates a VTK output file for visualization
* -e (optional: output_info_file_full_name, default value exists): exports errors and solver information
* -D: prints the ambient dimension
* -h: prints the help message
Main files are available in the Schemes directory. In particular, the two demo files are perfect to get started.
* ./bin/demo_poisson -m ../Meshes/gmsh_meshes/3Dcube_1_737tetra.msh
* ./bin/demo_stokes -m ../Meshes/gmsh_meshes/3Dcube_1_737tetra.msh
The examples below solve the Poisson problem with different methods (P3 Lagrange, P1 Crouzeix-Raviart, HHO(1,1), or VEM(1,2)). These methods are respectively of order 3, 1, 2, and 2.
* ./bin/pklagrange_poisson -m ../Meshes/fvca5_meshes/Lshape_tri2.typ2 -s 3 -v
* ./bin/p1cr_poisson -m ../Meshes/fvca5_meshes/Lshape_tri2.typ2 -e
* ./bin/hho_poisson -m ../Meshes/fvca6_meshes/dkershaw08.typ6 -b 1 -s 1 -v
* ./bin/vem_poisson -m ../Meshes/fvca6_meshes/dkershaw08.typ6 -b 1 -s 2
More information on the structure of the code is available using the Doxyfile in the doc directory (note that among the main files, only the demo ones are documented).
```