fdrmrc / polydeal

C++ implementation of Polygonal Discontinuous Galerkin method within the deal.II Finite Element library.

polyDEAL (Polytopal Discontinuous Galerkin in deal.II)

PolyDEAL is an open source library which aims to provide building blocks for the developement of Polytopal Discontinuous Galerkin methods, using the Finite Element library deal.II in 2D and 3D. It is written in C++ using the C++17 standard. The parallel implementation builds on top of the Message Passing Interface (MPI) communication model.

Getting started and prerequisites

We require:

• cmake version greater than 2.8.
• One of the following compilers:
  • gcc version >= 9.4.0
  • clang version >= 15
  • icc (Intel compiler) 2021.2
• openMPI version >= 4.0.3
• deal.II version >= 9.5

The library polyDEAL employs deal.II as main third-party library. As deal.II itself depends on other external libraries for many functionalities, we strongly suggest to download and install deal.II following the instructions available at https://www.dealii.org/download.html and https://www.dealii.org/developer/readme.html. The minimal set of other external libraries that we require are: METIS, p4est, Trilinos. All of them should be compiled against MPI during the installation phase of deal.II.

While METIS is generally used to partition a triangulation among several processors, in the context of polytopal methods it is heavily employed as an agglomeration strategy to build polytopic elements out of fine grids composed by standard shapes. Trilinos (in particular its multilevel solvers and distributed matrices) is employed as main parallel linear algebra library. We also support novel agglomeration strategies based on the R-tree spatial data structure. The associated preprint can be found online on arXiv.

To enable to computation of some quality metrics, mostly of theoretical interests and not really relevant in application codes, the external library CGAL is required. As this is a dependency of deal.II as well, it is sufficient to configure deal.II with it.

We currently support the following features:

• Unified interface for 2D and 3D problems.
• Distributed-memory implementation through MPI.
• Discontinuous Galerkin spaces of order $p$.
• Different agglomeration strategies.
• Multigrid support (still experimental)

Building polyDEAL

Assuming deal.II is installed on your machine and meets the requirements above, all is required to do is:

git clone [email protected]:fdrmrc/Polydeal.git
cd Polydeal/
mkdir build
cd build/
cmake -DDEAL_II_DIR=/path/to/deal.II ..
make -j<N>

being N is the number of jobs you want to use to compile.

About grids and post-processing

Polygonal or polyhedral grids are generated through agglomeration. To give an easy example in 2D, here's a polygonal grid associated to the unit ball on which you can define a Discontinous Galerkin space. The 3D case is completely analougous.

The (discontinuous) Finite Element space is defined on the bounding box of each agglomerate. In order to visualize the solution, carry out convergence tests, and other post-processing issues, the solution computed on the polytopal grid is interpolated onto the underlying fine mesh composed of classical quadrilaterals or hexahedra.

Examples

Some example applications are shown in the examples/ directory. To build and run one of the examples, say diffusion_reaction.cc, it is sufficient the following:

// assume you have a build generated as above
cd build/examples
make
mpirun -np<N> ./diffusion_reaction

where N is the number of processors you want to use.

Documentation

A Doxygen generated documentation is built and deployed at each merge to the main branch. You can find the latest documentation here: https://fdrmrc.github.io/Polydeal/.

Authors and Contact

This project is developed and maintained by:

• Marco Feder (@fdrmrc), SISSA, IT

under the supervision of

• Prof. Andrea Cangiani (@andreacangiani), SISSA, IT
• Prof. Luca Heltai (@luca-heltai), University of Pisa, IT

Feel free to start a discussion or open an issue, especially if you want to contribute. For any other inquiries or special requests, you can directly contact [email protected].

polydeal's Issues

Provide geometrical metrics for polygons

This paper https://arxiv.org/abs/1906.01627 provides interesting quality metrics for polygons. We should really compute some of them to assess the goodness of the implementation. In particular, I am planning to compute uniformity factor and circle_ratio for a general 2D polygon.

Missing 3D tests

Check things work also in 3D.

Bug when polygon's faces are not axis-aligned

With bilinear DG elements (i.e. DGQ(1)) the L2 error is not zero for the linear solution $u=x+y$ if the agglomerates have boundary faces that are not axis-aligned. I'm debugging this on a 4x4 grid where agglomerates are a bit distorted. The error is around 1e-5.

That's the grid:

I'm debugging quadrature points on faces and associated normals for this minimal example but they look okay.

Add proper build toolchain

In order to provide several app/examples that can be linked against polydeal and run without changing examples/CMakeLists.txt manually.

Bug in connectivity (neighboring faces of agglomerates)

Sanity checks fail because basis functions involved in the jump terms are evaluated at wrong quadrature points. Indeed, printing explicitely the evaluation points seen from the two cells sharing a same face $f$ between two agglomerates, they are sometimes different:

// print real qpoints from current_cell
qpoints: 
0.0140877 0.125
0.0625 0.125
0.110912 0.125
// print real qpoints from neighboring cell
qpoints :
0.125 0.139088
0.125 0.1875
0.125 0.235912

I'm working on that.

Give flags in proper way, without hardcoding

https://github.com/fdrmrc/AggloDeal/blob/00f56528a4624a5daee618e5aec892477c47d2bb/src/agglomeration_handler.cc#L215

Physical BBox not yet implemented in 3D

    Assert(dim == 2 && spacedim == 2, ExcNotImplemented()); //@todo Not working in 3D

https://github.com/fdrmrc/Agglomeration_Suite/blob/d4d067c43cec176676e515e0423b2a176dbdc944/include/agglomeration_handler.h#L372

Bounding boxes vector should be smaller and contiguous

The minimal change in #64 highlights something serious. In order to map points back from real to unit space, we should not go through the many calls of the Mapping class, but just use BoundingBox::real_to_unit(). Surprisingly, this increased the computational times.

The reason is that the std::vector of bounding boxes is long as the number of elements of the underlying triangulation. This should really be contiguous.

Iterators and Accessors classes

While the current status works, there is a major design aspect which has to be considered.

We should not let the user know explicitely that only one deal.II cell (the "master" one) of an agglomerate carries the DoFs. In our assembly routines, we effectively loop through all the master cells (i.e. all the polygons), but they are by all means simple deal.II cells. This continuous identification between polyongs and master cells is certainly confusing from a user's perspective.

The cleanest solution would be to design an Accessors and Iterators classes in the spirit of the Particles namespace from deal.II (https://www.dealii.org/current/doxygen/deal.II/namespaceParticles.html).

Floating point arithmetic

If we interpret each cell made by one element as a polygon itself, quantites like weights are scaled with the jacobian of the mapping associated to the bounding box and then multiplied back. This in necessary in order to conform with deal.II structures.

Consider the test of measuring the volume of a domain $\Omega$ just by integrating on $\Omega$ the unit function. If cells are perfectly axis-aligned squares describing the unit square, then considering classical deal.II cells as polygons doesn't make any difference. On the contrary, with a ball, the two approaches differ by the unit roundoff $\mu$. This has to be expected, as floating point arithmetic is not associative.

The simplest solution is to just use classical FEValues objects if we know this cell is made by just one element. I can confirm this give results that are equal up to machine precision

Run actions when opening PR

Make sure it compiles against latest version of deal.II, and run all tests.

Use new interface to agglomerate elements.

After #106, a new and more compact way to agglomerate elements is available. We should be consistent and use it everywhere.

The class CellsAgglomerator should be used. Moreover, it allows extracting the polytopal hierarchy in user codes.

Initialization of AgglomerationHandler

Constructor should take a Cached tria.

Order of quadrature rules on polytopes

In several examples we use quadrature rules of order $2p+1$ inside each element of the sub-tessellation, which are definitely an overkill, especially in 3D.
I've verified locally that decreasing the order (say to $p+1$ only) does not affect accuracy in our examples and gives much better assembly times thanks to better usage of caches.

Interpolation on fine, distributed, grid

After #80, the only missing step in order to provide a full working example is to provide a way to interpolate the (distributed) solution vector on the fine (distributed) triangulation, i.e. to implemente a parallel version of AgglomerationHandler<dim, spacedim>::setup_output_interpolation_matrix(), the rectangular matrix that provides the action of the nodal interpolant on the finer finite element space.

I already have a working version, but it needs to be polished in order to be used both in serial and parallel programs.

FEAggloValues?

AgglomerationHandler should not be responsible for the initialization and setup of FEValues-like objects, in the same way as DoFHandler does not evaluate any basis function. A dedicated class, with name something along the lines of FEAggloValues or FEPolyValues, should be the one that calls these methods:

• AgglomerationHandler<dim,spacedim>::reinit(polytope) ,
• AgglomerationHandler<dim,spacedim>::reinit(polytope,face)
• ...

Its constructor should be identical to the one of FEValuesBase with the twist that instead of using a DoFHandler, and AgglomerationHandler is working under the hood.

Penalty parameter

Convergence is guaranteed for a proper choice of the penalty function $\sigma(x)$:

Polytopic output

Given the solution u defined on the BBOxes of each polytope, interpolate in order to visualize it.

Use Signals to invalidate data

Unique face shared by neighboring elements

Whenever we have neighboring elements sharing lots of background faces we enumerate all of them and store this information.

• Consequentlely, if we have two neighboring cells that are sharing many facets we call reinit(polytope, f) for every face index f, which implies many calls. In practice, that means we create and reinit an object of type FEImmersedSurfaceValues (plus the related quadratures and normals needed by ImmersedSurfaceQuadrature) for every (background) face.
  On the other hand, the basis functions evaluated on these faces are always the same so it makes more sense to uniquely identify all of the faces sharing the same neighboring polygon as one single face.

Consider the following mesh made by two elements, and assume that K1 shares with K2 many straight little faces. With the latter approach:

• K1 has two faces: one shared with K2 and its boundary (composed of three lines),
• similarly for K2.

- - - - - - - - - - - - - -
|            |            |
|       K1   |     K2     |
- - - - - - - - - - - - - -

Many tests will fail due to diffs in the golden output.

• Moreover, as stated in the documentation of the class:

Thus, an FEImmersedSurfaceValues object can, typically, not be reused for different cells.

This suggests that it's better to group together these faces, since the cell we have now is conceptually just one polygonal cell.