optimad / bitpit

Open source library for scientific HPC
http://optimad.github.io/bitpit/
GNU Lesser General Public License v3.0
117 stars 34 forks source link

Problem with bitpit::CGElem::distancePointPolygon when compiling bitpit in Debug Mode #355

Closed kgkara closed 1 year ago

kgkara commented 1 year ago

The problem regards the usage of bitpit::CGElem::distancePointPolygon when compiling bitpit in Debug mode. An assertion is triggered by line /home/kostas/Optimad/bitpit/src/CG/CG_elem.cpp:66 because of line 572 of the same file. Here we need to project a point on the polygon in pink contour of this picture image when computing the minimum distance between the point and the surface. On the right hand side you might see the vertex coordinates comprising this polygon. The problem is that the formed triangle of the 3 pink colored points of this figure image has a zero surface area. You might see their coordinates on the right hand side. However, when it comes to computation of the minimum distance of a point from the surface, this normally should not make the code to stop, because the neighbour boundary interfaces don't lead to triangles of zero surface area. So, they could provide the minimum distance instead and the code could just neglect this problematic triangle. Alternatively, the triangles could be created with the barycenter of the polygon.

andrea-iob commented 1 year ago

All algorithms in CG assume that the elements are non-degenerate. There are some assertion that check if the elements are valid. The behaviour you are describing is therefore not a bug, but it's by design.

The coordinates of that triangle doesn't seems so close that the assertion should be triggered. Could you print the values of v0, v1, and v2 from "bool validTriangle(const array3D &P0, const array3D &P1, const array3D &P2 )"?

kgkara commented 1 year ago

All algorithms in CG assume that the elements are non-degenerate. There are some assertion that check if the elements are valid. The behaviour you are describing is therefore not a bug, but it's by design.

The coordinates of that triangle doesn't seems so close that the assertion should be triggered. Could you print the values of v0, v1, and v2 from "bool validTriangle(const array3D &P0, const array3D &P1, const array3D &P2 )"?

image

andrea-iob commented 1 year ago

MicrosoftTeams-image

The area is 1e-19 and that's why the triangle is marked as invalid. I don't understand why with such big differences in point positions the area is so small. Do you see an error in this check?

array3D n = crossProduct( v0, -1.*v2);
if( utils::DoubleFloatingEqual()( norm2(n), 0. ) ){
        return false;
}

If I'm not mistaken, its intent should be check if the area is close to zero.

kgkara commented 1 year ago

MicrosoftTeams-image

The area is 1e-19 and that's why the triangle is marked as invalid. I don't understand why with such big differences in point positions the area is so small. Do you see an error in this check?

array3D n = crossProduct( v0, -1.*v2);
if( utils::DoubleFloatingEqual()( norm2(n), 0. ) ){
        return false;
}

If I'm not mistaken, its intent should be check if the area is close to zero.

If I understand your question, the problem is that the 3 points are practically collinear, so they cannot define a triangle.

andrea-iob commented 1 year ago

Now I get it. I thought the edge was parallel to the y axis. However, it is slightly inclined.

The three points are collinear. The module CG will not see this as a triangle.

It is possible to add support to degenerate elements, but it need to be added to the whole CG module. Otherwise is a bit uselsee, because it may fix this problem, but then the same assertion will be triggered in other parts of the code.

We can add a new function that evaluates distances with respect to generic polygons, and that function can do a preliminary check to figure out what's the polygon type being passed and call the proper functions.

andrea-iob commented 1 year ago

I just realized that you are already calling distancePointPolygon (for unknown reasons I though you were using projectPointTriangle).

The algorithm that triangulates the polygon produces and invalid triangle. I will figure out how to fix it.

Unfortunately we have two version of the tesselation algorithm, one in CG an the other in PatchKernel. The one in CG is simpler/faster and can produce invalid triangles.