Fix algorithm in Manifold::normal_vector #5393
Conversation
|
I should add that the new algorithm uses the first combination of vertices that fulfills all three conditions, so it should not take much longer than the old algorithm. |
|
It seems this creates some issues at other places, let me look into that first before going forward with this PR. |
gassmoeller
force-pushed
the
fix_manifold_normal_vector
branch
from
d43fa9f
to
3279cf5
Compare
|
Ok, I figured it out, there were still issues with round-off and floating point comparisons. This should now work as intended. |
|
/run-tests |
|
Test |
source/grid/manifold.cc
Outdated
| else | ||
| // if p is too close to vertex i try again with different i and j | ||
| if ((p - face->vertex(i)).norm_square() < | ||
| std::numeric_limits<double>::epsilon() * (p - face->vertex(j)).norm_square()) |
There was a problem hiding this comment.
I think you should err on the safe side and make this 100 * eps or some such. In fact, since you're comparing squares, you may even want to use 10000 * eps.
source/grid/manifold.cc
Outdated
| if ((p-face->center()).norm_square() < min_distance) | ||
| // Look for a combination of tangent vectors that | ||
| // are of approximately equal length and not linearly dependent | ||
| for (unsigned int i=0,j=1 ; i<4 && j<4; ++j) |
There was a problem hiding this comment.
Can you use GeometryInfo<2>::vertices_per_cell here?
source/grid/manifold.cc
Outdated
|
|
||
| // if p is too close to vertex j try again with different j | ||
| if ((p - face->vertex(j)).norm_square() < | ||
| std::numeric_limits<double>::epsilon() * (p - face->vertex(i)).norm_square()) |
source/grid/manifold.cc
Outdated
| normal = cross_product_3d(t1,t2); | ||
|
|
||
| // if t1 and t2 are (nearly) linearly dependent try again with different j / t2 | ||
| if (normal.norm_square() < std::numeric_limits<double>::epsilon() * |
tests/manifold/flat_manifold_08.cc
Outdated
| @@ -0,0 +1,64 @@ | |||
| //---------------------------- manifold_id_01.cc --------------------------- | |||
| // Copyright (C) 2011 - 2015 by the mathLab team. | |||
source/grid/manifold.cc
Outdated
| const Tensor<1,spacedim> reference_normal = cross_product_3d(reference_t1,reference_t2); | ||
|
|
||
| if (reference_normal * normal < 0.0) | ||
| normal *= -1; |
There was a problem hiding this comment.
Can you explain why you need this? Is it because you've lost track of whether vertices i and j are oriented in clockwise orientation with an angle of less than 180 degrees when seen from point p? That's when you get the sign flip. But I wonder whether you can detect that case without having to call get_tangent_vector twice as that might be an expensive operation.
source/grid/manifold.cc
Outdated
| continue; | ||
|
|
||
| t1 = get_tangent_vector(p, face->vertex(i)); | ||
| t2 = get_tangent_vector(p, face->vertex(j)); |
There was a problem hiding this comment.
get_tangent_vector may be an expensive operation. Do you think you could cache these vectors and only update them when you change i or j (and have decided that you actually need them)?
gassmoeller
force-pushed
the
fix_manifold_normal_vector
branch
2 times, most recently
from
bcb4048
to
f922ef2
Compare
|
After a number of iterations I hope the new algorithm improves on the predecessors. In particular it should be nearly as fast, does not fail for special cases anymore, and preserves existing test results. Unfortunately, it does not preserve existing functionality completely. For points inside of faces of 3D triangulations that do not form a plane (i.e. they are twisted) this algorithm sometimes uses different vertices than the previous implementation to compute tangent vectors, and thus the resulting normal vectors are different. I added test numerics/no_flux_14 to document this behavior, it passes on this branch but fails on master. |
gassmoeller
force-pushed
the
fix_manifold_normal_vector
branch
from
f922ef2
to
faeece0
Compare
|
Yes, excellent. Very nice and clean code, thank you! |
I think the current algorithm to compute normal vectors in the
Manifoldclass does not work, at least for some special cases (and this is the reason for the bug in geodynamics/aspect#1990).Consider the following flat face (i.e. it is not a perspective drawing) of a 3D cell with vertices 0,1,2,3, and assume we want to compute the normal vector for the point at the center of the bottom line:
The current algorithm separates the face into a diamond close to the center, and triangles closest to the vertices. Our point in question here is closer to the center than to any vertex. The current algorithm then looks for the closest vertex (vertex 2 in this case), and based on this index selects two other vertices as the ones to form tangent lines to. But in the current version it will select vertices 0 and 1, which leads to the case that tangent lines t1 and t2 are linearly dependent. It then takes the cross product of t1 and t2, which leads to a zero normal vector, which leads to a division by zero when normalizing the normal vector.
I here propose a new algorithm that does not try to geometrically guess which vertices to use, but instead tries all available vertex combinations, and only ensures that the following conditions are met by the normal vector:
The included test reconstructs exactly the case I have drawn above and crashes with a division by zero on current master, but produces the expected result with this patch.