Issue with intersection of two meshes with curved boundaries

[Najwa Alshehri]

Dear developers and users,

 

I have two meshes one is immersed in the other. I wanted to find the intersection between the two meshes, so I used the following function.

 

    NonMatching::compute_intersection(omega_grid_tools_cache,

                                      omega2_grid_tools_cache,

                                      4,      // degree

                                      1e-20); // to

include header

source code

This function uses CGAL to find a quadrature formula to integrate exactly on the intersection of two meshes which neglecte intersections of areas with tolerance smaller than “tol“ that one chooses and gives a quadrature formula on the triangulation of the intersection area that integrates exactly polynomials of a specific degree ( which allows maximum degree of 4).

Say that I want to integrate, on the intersection area, a polynomial of order 4.

 

I noticed that If I am considering a circular immersed domain (unlike square or L-shaped domains), after a few cycles, the quadrature formula is not accurate enough. To be precise, I find that the sum of the weights of the quadrature formula defined on the triangulation of the exact intersection of the two meshes does not sum up (up to a tolerance) to the measure of the domain. When this occurs, the solution that is solved by evaluating the integral considering the quadrature formula on the exact intersection is no longer correct and the error starts to diverge in the later cycles after this point.

 

Moreover, the difference gets large suddenly, in one cycle, the difference was relatively smaller (1e-13), and in the next, it is much larger (1e-8) as can be seen in the attached plot (Plot shows the difference between the sum of the weights on the whole domain defined on the triangulation of the exact intersection of the two meshes and the measure of the domain under uniform refinement of the mesh). 

 

Any suggestion on what could be the issue and what should I do to fix it?

 

Thanks

Najwa

[Najwa Alshehri]

Dear team,

I am adding the two domain meshes (cycle 0) in .vtk format for a better understanding of what type of grid I am working with

[Marco Feder]

After debugging this particular instance, it turns out the issue is related to the CGAL kernels used by deal.II. In particular, with CGAL it's possible to rely on the so-called "exact computation paradigm" (a brief explanation is available here: https://www.cgal.org/exact.html). In the dealii::CGALWrappers namespace this kernel is **not** used for the quad-quad intersections, which is the one relevant for this example. Switching to exact kernels allows to keep the error around 1e-13 for most of the refinement cycles with this particular configuration.

Maybe it could be nice to provide a policy that allows the used to simply tweak the CGAL kernels that are employed.

Best,

Marco