Giter Site home page Giter Site logo

fdrmrc / polydeal Goto Github PK

View Code? Open in Web Editor NEW
0.0 1.0 0.0 81.32 MB

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

Home Page: https://fdrmrc.github.io/Polydeal/

License: Other

CMake 0.31% C++ 98.07% Dockerfile 0.14% Shell 0.19% GLSL 1.09% Python 0.21%
finite-element-methods scientific-computing discontinuous-galerkin agglomeration cpp-library polygonal-elements meshes polygonal-meshes

polydeal's People

Contributors

fdrmrc avatar luca-heltai avatar

Watchers

 avatar

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.

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.

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.

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:
immagine

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.

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)$:

Screenshot 2023-11-18 at 13 13 16

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

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    ๐Ÿ–– Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. ๐Ÿ“Š๐Ÿ“ˆ๐ŸŽ‰

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google โค๏ธ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.