HArD::Core2D
Hybrid Arbitrary Degree::Core 2D  Library to implement 2D schemes with edge and cell polynomials as unknowns

HArD::Core (sources: https://github.com/jdroniou/HArDCore2D) provides a suite of tools, in C++, to implement numerical schemes whose unknowns are polynomials in the cells and on the edges (in 2D) or faces (in 3D). The focus is on dealing on generic polytopal meshes. This documentation addresses the 2D version of HArD::Core, but similar principles are valid for the 3D version.
To build the libraries and implemented schemes, the minimal requirements are:
Make sure that you have the development version of boost installed. On Linux, install libboostdev
, libboostfilesystemdev
, libboostprogramoptionsdev
, libboostchronodev
and libboosttimerdev
from your package manager.
The linear systems resulting from the assembled scheme are solved using the BiCGStab implementation of Eigen. An alternative (currently commented out in the schemes' implementations) is to use the MA41 solver of the HSL library. To use this alternative you will need:
Once you have installed all of the required dependencies, set up the build directory and generate the build files by running the following from the repository root:
After this, build/Schemes
will contain the executables (e.g. hhodiffusion
) to run the schemes. These executable need to access the typ2 meshes, which they should naturally find if you put the typ2_meshes
directory at the root of the project's files.
The mesh documentation is built with Doxygen (http://www.stack.nl/~dimitri/doxygen/). If you are reading this then somebody has already built it for you. If you modify the code and wish to rebuild the documentation, simply run doxygen
from the root directory. The HTML version of the documentation is generated inside documentation/html
and the LaTeX version is generated inside documentation/latex
and can be compiled using the generated Makefile.
After it is loaded, the mesh is represented by classes describing a vertex, an edge, and a cell: Vertex, Edge, and Cell. Each of these classes contains methods to access useful information for the corresponding element, including other geometrical quantities it is related to. The mesh itself is represented by an element of the Mesh class with methods to access all the vertices, edges and cells (or a particular vertex, edge or cell). In this class, each cell has a unique identifier, numbered from 0.
For example, if mesh_ptr
is a pointer to a Mesh class,
would store the coordinates of the fifth vertex into the Eigen vector vert_coord. As a generic rule, all geometrical vectors are Eigen::Vector2d
. We also use Eigen::Vector{2,X}d
and Eigen::Matrix{2,X}d
for objects on which linear algebraic operations are performed. Lists (e.g. of cells, of functions...) are instances of std::vector<...>
. Finally, Eigen::Array
is used for lists on which componentwise operators are performed (typically, values of some functions at the quadrature nodes, that will be multiplied componentwise by the corresponding quadrature weights and then summed together).
Here is an example that loops over all cells, grabs all the edges of the cell, and loops over these edges to output their length. Here, mesh_ptr
is a pointer to the mesh.
The mesh classes and other auxilliary classes are located inside the namespace HArDCore2D.
There is no direct access from a highlevel geometrical entity to elements purely associated with lowerlevel entities. For example, if mesh_ptr
is a pointer to the mesh, there is no direct method to access the coordinates of the ith vertex of the mesh (no mesh_ptr>coords_vertex()
exists). Instead, this is done through mesh_ptr>vertex(i)>coords()
. This choice is deliberate as it preserves the logical organisation of the data structure, and facilitates the memorisation of the various methods. Of course, writing a wrapper providing a direct access is easy...
HArDCore2D can currently read meshes in the typ2
format designed for the FVCA5 Benchmark. A short documentation describing this format is provided in the typ2_meshes
directory (see README.pdf there). Several meshes can also be found in this directory.
A mesh file must be read using an instance of the MeshReaderTyp2
class, and then built using MeshBuilder
. A working example is given below (assuming the executable will be in build/Schemes
for example).
Note that the builder returns a unique_ptr
for the mesh. This ensures that, at the end of the code, the mesh destructor is called (which destroys all cells, edges, vertices...). Classes and functions use a raw pointer to the mesh, so the .get()
method should be used when passing the mesh as argument to the class constructors or functions (see the HybridCore
example above).
Note: the typ2
format allows for meshes with very generic polygonal cells, including nonconvex cells. However, the builder assumes that each cell is starshaped with respect to the isobarycenter of its vertices – otherwise, the calculation of the center of mass may be incorrect. Similarly, the quadrature rules (see Quadrature rules) assume that each cell is starshaped with respect to its center of mass.
The HybridCore structure encapsulates routines to create bases of polynomial spaces in each cell and on each edge, to integrate functions on these mesh entities, and to evaluate functions defined through their coefficients on the cell and edge basis functions. Start by including the structure as follows:
To initialise a HybridCore structure, we must first have a mesh loaded as a pointer (see loading a mesh). The structure is then created by specifying the desired degree of the polynomial spaces on the edges and in the cells.
Initialising the HybridCore structure constructs basis functions for the polynomial spaces on the edges (up to degree ) and in the cells (up to degree ). The choice of the degree in the cells corresponds to the needs of certain highorder methods, such as the HHO method that requires the reconstruction of a polynomial of degree in each cell. It is also assumed that .
The basis functions are accessed through the methods cell_basis(iT,i) and edge_basis(iE,i) which return the ith basis function on the iTth cell or iEth edge. The cell gradients are available from cell_gradient(iT,i); these gradients are indexed to correspond with their basis functions, which means that the first gradient will always identically be the zero vector, since it corresponds to the constant basis function.
The basis functions are hierarchical, which means that they are constructed by increasing degree. A basis of the space of polynomials of degree in the cell is thus obtained by selecting the first cell basis functions.
When a scheme has polynomial unknowns of degree on the edges and in the cells, these unknowns can be represented as vectors of coefficients on the basis functions, for example by listing all the coefficients on the basis functions in the first cell, then all the coefficients on the basis functions in the second cell, etc., and then listing all the coefficients on the basis functions in the first edge, etc. This is the choice adopted in HHOdiffusion.
A number of convenient quantities relating to the basis functions are available in the HybridCore structure follows.
Symbol name  Meaning 

nlocal_cell_dofs  The number of degrees of freedom of a cell polynomial (ie. the dimension of the space of polynomials of degree in two variables) 
nlocal_face_dofs  The number of degrees of freedom of a face polynomial (ie. the dimension of the space of polynomials of degree in one variable) 
nhighorder_dofs  The number of degrees of freedom of a degree cell polynomial (ie. the dimension of the space of polynomials of degree in two variables) 
ngradient_dofs  The dimension of the gradient space of degree cell polynomials 
ntotal_cell_dofs  The total number of degrees of freedom over all cell polynomials over the entire mesh (i.e. the number of cells times the dimension of the space of polynomials of degree in two variables) 
ntotal_face_dofs  The total number of degrees of freedom over all face polynomials over the entire mesh (i.e. the number of edges times the dimension of the space of polynomials of degree in one variables) 
ninternal_face_dofs  The total number of degrees of freedom over all internal faces over the entire mesh 
ntotal_dofs  The total number of cell and face degrees of freedom over the entire mesh 
The HybridCore structure provides routines to integrate generic functions on the cells and the edges. These routines are however expensive as they recompute the quadrature nodes and weights every time they are called. They should therefore only be used with parsimony; computing quadrature nodes and values of basis functions at these nodes is more efficient, see Quadratures rules.
For example, to integrate over the cell number iT:
Basis functions can also be integrated:
HArD::Core deals with quite arbitrary cell geometries. As a consequence, no reference element can be used, and the quadrature rules have to be adjusted to each particular cell. This is done by partitioning each cell into triangles and by using John Burkardt's implementation of the Dunavant rules. The choice was also made not to precompute all quadrature rules for all cells and edges, but rather to compute them – with a locally chosen degree of exactness – when needed in the code. To reduce the computational cost, quadrature rules – and the values of basis functions at quadrature nodes – should only be computed once when looping over each cell, before being used, e.g., to form mass matrices.
The HybridCore structure provides routines to do that. The method cell_qrule(iT,doe) calculates quadrature nodes and weights, exact up to the polynomial degree doe
for an integration over cell number iT
; see edge_qrule(iE,doe) for the equivalent over an edge. At present, the quadrature rules available in the code support a total degree of exactness in the cells or on the edges up to 20. This quadrature rule, stored for example in quadTE
, is then be provided to basis_quad(type,i,quadTE,deg) which computes the values of the basis functions in cell/edge (depending on type
=T or some other character) number i
up to the specified degree deg
; see also grad_basis_quad(i,quadTE,deg) to compute the gradients of the basis functions at the quadrature nodes.
Typically, these values are then passed on to gram_matrix(f_quad,g_quad,Nf,Ng,quadTE,sym,weight) which computes a ``generalised'' Gram matrix of the families of functions and , provided at the quadrature notes by f_quad
and g_quad
(here, sym
is a boolean indicating if the matrix is expected to be symmetric). This Gram matrix method is useful to compute mass and stiffness matrices.
Here is an example.
The quadrature rules and values of basis functions (and gradients) at the quadrature nodes can be conveniently computed and stored using the CellEdgeQuad class. Instantiating an element of this class on a cell loads these rules and values once, that can then be passed to several functions in charge of various calculations (e.g. one function computes the local cell contribution to the diffusion term, another function is in charge of computing the load term associated to the cell, etc.). This prevents recomputing these rules and values when needed by various functions. It works the following way:
The following schemes are currently available in HArD::Core2D. The Hybrid HighOrder schemes follow the implementation principles described in Appendix B of the book available at https://hal.archivesouvertes.fr/hal02151813.
The directory runs
contains BASH to run series of tests on families of meshes. The files data.sh
describe the parameters of the test cases (polynomial degrees, boundary conditions, mesh families, etc.). The script produces results in the output
directory, including a pdf file rate.pdf
describing the rates of convergence in various energy norms.
To run the scripts as they are, you will need pdflatex
and a FORTRAN compiler, and to adjust the Makefile
to your compiler, to run compute_rates.f90
and compute the rates of convergence in the various norms.