fdrmrc/polydeal

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

polydeal's Issues

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.

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.

Run actions when opening PR

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

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.

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.

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.

Give flags in proper way, without hardcoding

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

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.

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.

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

Initialization of AgglomerationHandler

Constructor should take a Cached tria.

Missing 3D tests

Check things work also in 3D.

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).

Polytopic output

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

Penalty parameter

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

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.

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.

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

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.

Use Signals to invalidate data