|
HArD::Core3D
Hybrid Arbitrary Degree::Core 3D - Library to implement 3D schemes with vertex, edge, face and cell polynomials as unknowns
|
HArD::Core (sources: https://github.com/jdroniou/HArDCore) provides a suite of C++ tools to implement numerical schemes whose unknowns are polynomials in the cells, on the edges, and on the faces. The focus is on dealing on generic polytopal meshes. This documentation addresses the 3D version of HArD::Core, but similar principles are valid for the 2D version. Transferring a scheme's implementation from 3D to 2D or vice-versa is very straightforward, provided that the scheme's mathematical definition does not depend on the dimension and that the generic types provided in basis.hpp and MeshObject.hpp are used; see readme file of the HArD::Core github depository https://github.com/jdroniou/HArDCore.
You will find on this page the following sections:
A tutorial describes the principles for implementing a scheme using the HArDCore3D tools.
A template for the main CMakeLists is provided in the repository. Copy CMakeLists-template.txt into CMakeLists.txt and edit the later according to your compilation needs.
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 libboost-dev, libboost-filesystem-dev, libboost-program-options-dev, libboost-chrono-dev and libboost-timer-dev from your package manager.
The provided CMakeLists.txt file should take care of installing these minimal libraries in case they are not found on your system; see INSTALL_NOTES.txt for details.
The default solver for linear systems is the LU implementation of Eigen, but alternative can easily be selected via the class LinearSolver. Some of these alternative (UMFPACK, SUPERLU, Pastix...) require additional libraries, which can be installed by passing the proper arguments to cmake (see INSTALL_NOTES.txt and the provided CMakeLists file).
A couple of schemes have procedures to compute condition numbers of matrices; these require the Spectra library (https://spectralib.org/).
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. hho-diffusion) to run the schemes. These executables need to access the meshes, which they should naturally find if you have left the 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 docs and the LaTeX version is generated inside docs/latex and can be compiled using the generated Makefile.
After it is loaded, the mesh is represented by typedefs of MeshObject describing a Vertex, an Edge, a Face, and a 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, faces and cells (or a particular vertex, edge, face or cell).
For example, if mesh_ptr is a pointer to a Mesh instance, the lines
store the coordinates of the fifth vertex into the Eigen vector vert_coord. As a generic rule, all geometrical vectors are Eigen::Vector3d. We also use Eigen::Vector{3,X}d and Eigen::Matrix{3,X}d for objects on which linear algebraic operations are performed. Lists (e.g. of cells, of functions...) are usually instances of std::vector<...>. Finally, Eigen::multi_array is used for lists of values of functions at quadrature nodes.
Here is an example that loops over all cells, grabs all the faces of the cell, and loops over these faces to output their diameter. Here, mesh_ptr is a pointer to the mesh.
There is no direct access from a high-level geometrical entity to elements purely associated with lower-level entities. For example, if mesh_ptr is a pointer to the mesh, there is no direct method to access the coordinates of the i-th 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 such a direct access is easy...
HArDCore3D can read meshes in RF format. Previous versions could read TG and MSH files and were based on G. Manzini's mesh library https://github.com/gmanzini-LANL/PDE-Mesh-Manager. From Version 4.1, HArDCore3D uses an independent mesh reader written by L. Yemm.
A mesh file must be read using an instance of the meshbuilder class, and then built using build_the_mesh. 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 execution, the mesh destructor is called (which destroys all cells, faces, edges, vertices...). Some 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.
Note: the mesh formats allow for meshes with very generic polygonal cells, including non-convex cells. However, the builder assumes that each cell is star-shaped 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 star-shaped with respect to its center of mass.
Several other codes (in src/Mesh/Transforms) are provided to manipulate meshes:
The main classes in the Common describe polynomial basis functions on a cell, face, or edge. These could be bases of full polynomial spaces \(\mathbb{P}^k\), or other related spaces (vector-valued polynomials, subspaces, image of gradient, image of curl, complements, etc.). The underlying basis functions are monomial (MonomialScalarBasisCell and MonomialScalarBasisFace), but derived bases (or set of non-necessarily linearly independent polynomial functions) can be handled through various classes, such as Family, TensorizedVectorFamily, GradientBasis, etc.
Free functions are also available to compute basis functions at quadrature nodes (using the Quadrature rules module), orthonormalise basis functions, and compute Gram-like matrices between various families of functions. These matrices are essential in the design of high-order methods on polytopal meshes. Again, see examples in the HybridCore, DDRCore and the various schemes built on them.
This module also contains:
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/face/edge. For the cells, for example, this is done by partitioning each cell into tetrahedra and by using classical quadrature rules on tetrahedras. The choice was also made not to pre-compute all quadrature rules for all cells, faces 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 Quadratures module provides routines to do that. The key method, generate_quadrature_rule(CFE,doe), calculates quadrature nodes and weights, exact up to the polynomial degree doe, for an integration over cell/face/edge CFE (passed as a reference of Cell, Face or Edge class). At present, the quadrature rules available in the code support a total degree of exactness in the cells up to 14 in the cells and 20 on the faces and the edges (the quadrature rules on the faces come from John Burkardt's implementation of the Dunavant rules). The generated quadrature rule is stored in a structure QuadratureRule, which is a vector of quadrature nodes (weights and position).
The evaluate_quad template function evaluate basis functions (their value, gradients, etc.) at provided quadrature nodes. The boost::multi_array provided by this function can then be passed to compute_gram_matrix to create a matrix of inner products of two families of basis functions. Here is an example.
Note the usage of the type VectorRd defined in basis.hpp, which enables for a dimension-independent piece of code (easier to adapt to the 2D case). This procedure can also be applied, e.g., to cell basis functions on face quadrature nodes, etc. Additionally, the values at quadrature nodes obtained via evaluate_quad can be transformed using transform_values_quad (see also scalar_product and vector_product); this gives for example an easy way of constructing the values at quadrature nodes on a face of normal or tangential traces of cell polynomials.
Because of the generality of some polytopal elements and the resulting need to generate quadrature rules by splitting the elements into tetrahedra, quadrature rules can very quickly reach a very large number of nodes (more than 1000 on some Voronoi cells for higher degrees of exactness). The calculation of integrals, and in particular of Gram matrices, then represent a very large amount of the computational time. To remediate this, we have implemented in HArDCore3D the HNI approach for monomials of https://doi.org/10.1007/s00466-015-1213-7. Since all the polynomial bases in the Common module are built from monomial polynomials, this technique actually enables us to very efficiently compute Gram matrices of many of these bases.
The function GramMatrix precisely takes care of computing gram matrices of bases and derived sets of cell- and face-polynomials. At the moment, it properly manages all the following cases.
| Scalar bases: | MonomialScalarBasisCell || MonomialScalarBasisCell | |
| DivergenceBasis<T> || MonomialScalarBasisCell | where T=TensorizedVectorFamily or RolyComplBasisCell | |
| DivergenceBasis<T1> || DivergenceBasis<T2> | where T1,T2=TensorizedVectorFamily or RolyComplBasisCell | |
| Vector bases: | TensorizedVectorFamily || T | where T is any basis with rank=Vector (TensorizedVectorFamily, GradientBasis, GolyComplBasisCell, etc.) |
| GradientBasis<T1> || GradientBasis<T2> | where T1, T2 are any scalar bases | |
| CurlBasis<T1> || CurlBasis<T2> | where T1,T2=TensorizedVectorFamily or GolyComplBasisCell | |
| GolyComplBasisCell || GolyComplBasisCell | ||
| RolyComplBasisCell || RolyComplBasisCell | ||
| Matrix bases: | MatrixFamily<T1, N> || MatrixFamily<T2, N> | where T1, T2 are any scalar bases, and N is an integer |
| DivergenceBasis<MatrixFamily<T1, dimspace>> || TensorizedVectorFamily<T2, dimspace> | where T1, T2 are any scalar bases | |
| GradientBasis<TensorizedVectorFamily<T1, dimspace>> || MatrixFamily<T2, dimspace> | where T1, T2 are any scalar bases | |
| GradientBasis<TensorizedVectorFamily<T1, N>> || GradientBasis<TensorizedVectorFamily<T2, N>> | where T1, T2 are any scalar bases, N is an integer |
| Scalar bases: | MonomialScalarBasisFace || MonomialScalarBasisFace | |
| DivergenceBasis<T> || MonomialScalarBasisFace | where T=TangentFamily or RolyComplBasisFace (and derived) | |
| Vector bases: | TensorizedVectorFamily || TensorizedVectorFamily | |
| TangentFamily || TangentFamily | ||
| CurlBasis || CurlBasis | ||
| CurlBasis || TangentFamily | ||
| TangentFamily || RolyComplBasisFace | ||
| RolyComplBasisFace || RolyComplBasisFace | ||
| GolyComplBasisCell || GolyComplBasisCell |
This module encapsulates routines to create bases of polynomial spaces in each cell, on each face, and on each edge, and to manipulate discrete functions through the class UVector.
The instantiation of an HybridCore class creates basis functions for the full polynomial spaces \(\mathbb{P}^k\) in the cells, on the faces and on the edges, specifying the maximum degree required for each geometric entity. The basis functions are elements of the Family class and are accessed via the CellBasis, FaceBasis and EdgeBasis method. For example, the following piece of code initialises an HybridCore instance with degrees \(K+1\) in the cells, \(K\) on the faces, and no edge basis functions, and access the value at some Eigen::Vector3d x of the i-th basis function on face iF, and the gradient of the j-th basis function in cell iT.
The basis functions are hierarchical, which means that they are constructed by increasing degree. Hence, for example, if cell basis functions up to degree \(K+1\) have been generated, a basis of the space of polynomials of degree \(K\) in the cell is thus obtained by selecting the first PolynomialSpaceDimension<Cell>::Poly(K) cell basis functions.
The UVector class describes coefficients on cell and face basis functions. The first coefficients correspond to cell basis functions, ordered by the cells themselves, and the last coefficients correspond to face basis functions. The methods in this class provide the surrounding structure to make sense of these coefficients (degrees of considered polynomial functions in cells/on faces, restrictions to a certain cell and its faces, etc.).
As explained above, the evaluate_quad template enables the evaluation of basis functions (or their gradient, curl, or divergence) at pre-computed quadrature nodes. In the HybridCore module, however, the quadrature rules and values of basis functions (and gradients) at the quadrature nodes can be conveniently computed and stored using the ElementQuad 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 HHO3D module provides typedefs, a class and functions to implement Hybrid High Order (HHO) schemes. Rules (functions) to create local bilinear forms (matrices) and loading terms (vectors) are passed to the HHO3D class, that takes care of the global assembly and solving of the system.
The LocalDOFSpace class provides methods to manage degrees of freedom (DOFs) attached to a cell, its faces, its edges and its vertices. The global version is GlobalDOFSpace, and there is also a version, VariableDOFSpace, to handle cases where the number of DOFs varies, e.g., from one vertex to the other, from one edge to the other, etc. These classes provide methods to access the local or global degrees of freedom (determining the correct offset), etc.
The DOFs in these classes are organised by increasing order of the mesh entities dimensions (DOFs linked to all vertices, then DOFs linked to all edges, then DOFs linked to all faces, and finally all DOFs linked to cells). These DOFs only make sense when bases of polynomial spaces have been chosen on the geometric entities (such as some of the bases created by the HHOSpace or DDRCore classes below), and correspond then to the coefficients on these bases.
The following modules are based on these classes.
The HHOSpace module provides classes and functions to construct the local basis functions and operators for HHO space, including vector-valued spaces.
The DDRCore module provides classes and functions to implement the discrete de Rham (DDR) complex. This complex is based on spaces with unknowns on all geometric entities (vertices, edges, faces and cells), and discrete differential operators acting between these spaces. It is based on the principles detailed in https://doi.org/10.1007/s10208-021-09542-8 (see also the founding work https://doi.org/10.1142/S0218202520500372, and the arxiv versions https://arxiv.org/abs/2101.04940, https://arxiv.org/abs/1911.03616).
The main elements of the DDRCore module are:
Note that DDRSpace could potentially be used for generic schemes (not just based on the discrete de Rham sequence), perhaps with a different ordering of the DOFs.
Important note: a directory "DDRCore-orth" can be found in the repository. It corresponds to the DDR spaces using orthogonal complements to the images of the gradient and curl operators on polynomial spaces, as described in https://doi.org/10.1142/S0218202520500372. This directory is not commented here, and its code is not maintained any longer. The directory "DDRCore" is the one referred to in this documentation, is based on the Kozsul complements of the images of gradient and curl, as in https://arxiv.org/abs/2101.04940, and is the one that is still maintained. DDRCore-orth is only provided for comparison purposes; the complements in this directory are much more expensive to create and manipulated, as explained in https://arxiv.org/abs/2101.04946. To compile a scheme using the orthogonal complements, simply modify the main CMakeLists.txt and change all "DDRCore" into "DDRCore-orth".
The following schemes are currently available in HArD::Core3D. The Hybrid High-Order schemes follow the implementation principles described in Appendix B of the book available at https://hal.archives-ouvertes.fr/hal-02151813.
The directory runs contains BASH scripts 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, shows the convergence rate in the standard console output, and creates a pdf file rate.pdf describing the rates of convergence in various energy norms (you will need pdflatex to create this pdf file; commenting out the corresponding line is fine, the pdf will simply not be create).
1.9.8