From: abn Date: Thu, 25 Jan 2024 19:37:39 +0000 (+0100) Subject: [TetraIntersec] Several bug fixes in Grandy's implementation X-Git-Tag: V9_13_0a1~23 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=c84e8386ca771a4e6450ec178263961feefcb0bb;p=tools%2Fmedcoupling.git [TetraIntersec] Several bug fixes in Grandy's implementation + triple product inconsistency was not properly detected (now using original deltas from douple prod computation) + f/F factor set to 20 (not 500) as in Grandy's article -> (will be set to 100, see next commit) + epsilonEqual used when necessary: when it provides more points in polygon A or B (risk is completely missing a point in polygon!) + better handling of degenerated case where PQR triangle is in h=0 plane, or when P,Q or R close to a tet corner + fixed surface-edge intersection test: triple product equality (and zeroing) must be checked with care + using 'long double' is not necessary -> double are enough --- diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index d10d28ca7..5e1a08606 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -122,13 +122,15 @@ namespace INTERP_KERNEL _coords[5*Q + 3] = 1 - q[0] - q[1] - q[2]; _coords[5*R + 3] = 1 - r[0] - r[1] - r[2]; + // Handle degenerated case where one of the seg of PQR is (almost) inside XYZ plane, + // and hence by extension when the whole PQR triangle is in the XYZ plane + handleDegenerateCases(); + // H coordinate _coords[5*P + 4] = 1 - p[0] - p[1]; _coords[5*Q + 4] = 1 - q[0] - q[1]; _coords[5*R + 4] = 1 - r[0] - r[1]; - resetNearZeroCoordinates(); - // initialise rest of data preCalculateDoubleProducts(); @@ -785,7 +787,7 @@ namespace INTERP_KERNEL { // coordinate to check const int coord = static_cast(facet); - return (_coords[5*P + coord] == _coords[5*Q + coord]) && (_coords[5*P + coord] == _coords[5*R + coord]); + return (epsilonEqual(_coords[5*P + coord], _coords[5*Q + coord])) && (epsilonEqual(_coords[5*P + coord], _coords[5*R + coord])); } /** diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index bd658c9b4..b2ce44315 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -70,6 +70,9 @@ namespace INTERP_KERNEL * coordinates of the corners of the triangle. It copies their coordinates and then proceeds to pre-calculating certain * entities used in the intersection calculation : the double products, triple products and the values of the function E * (Grandy, [53]). + * It is also at this point in constructor that: + * - the special case of PQR included in the XYZ plane is treated + * - the inconsistencies between double products/triple products computation is handled * * calculateIntersectionVolume() : * This is the only method in the public interface. It calculates the volume under the intersection polygons @@ -95,18 +98,14 @@ namespace INTERP_KERNEL * When an intersection point has been detected it is calculated with a corresponding calc* - method in the cases where it * is not known directly. It is then added to the polygon A and/or B as necessary. * - * OPTIMIZE : - * If OPTIMIZE is defined, a large number of methods will be prefixed with inline and some optimizations concerning the tests - * with zero double products will be used. */ class INTERPKERNEL_EXPORT TransformedTriangle { public: - - friend class INTERP_TEST::TransformedTriangleTest; friend class INTERP_TEST::TransformedTriangleIntersectTest; + friend class INTERP_TEST::TransformedTriangleTest; /* * Enumerations representing the different geometric elements of the unit tetrahedron * and the triangle. The end element, NO_* gives the number of elements in the enumeration @@ -140,126 +139,90 @@ namespace INTERP_KERNEL double calculateIntersectionVolume(); double calculateIntersectionSurface(TetraAffineTransform* tat); - void dumpCoords() const; // Queries of member values used by UnitTetraIntersectionBary - const double* getCorner(TriCorner corner) const { return _coords + 5*corner; } - const std::vector& getPolygonA() const { return _polygonA; } - double getVolume() const { return _volume; } protected: - TransformedTriangle() { } // ---------------------------------------------------------------------------------- // High-level methods called directly by calculateIntersectionVolume() // ---------------------------------------------------------------------------------- void calculateIntersectionAndProjectionPolygons(); - void calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter); - void sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter); - double calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter); // ---------------------------------------------------------------------------------- // High-level methods called directly by calculateIntersectionSurface() // ---------------------------------------------------------------------------------- void calculateIntersectionPolygon(); - double calculateSurfacePolygon(); // ---------------------------------------------------------------------------------- // Detection of degenerate triangles // ---------------------------------------------------------------------------------- - bool isTriangleInPlaneOfFacet(const TetraFacet facet) const; - bool isTriangleParallelToFacet(const TetraFacet facet) const; - int isTriangleInclinedToFacet(const TetraFacet facet) const; - bool isTriangleBelowTetraeder() const; // ---------------------------------------------------------------------------------- // Intersection test methods and intersection point calculations // ---------------------------------------------------------------------------------- - inline bool testSurfaceEdgeIntersection(const TetraEdge edge) const; - void calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const; - inline bool testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const; - void calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const; - bool testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const; - void calcIntersectionPtSegmentEdge(const TriSegment seg, const TetraEdge edge, double* pt) const ; - bool testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const ; - inline bool testSurfaceRayIntersection(const TetraCorner corner) const; - bool testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edg); - void calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const; - bool testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const; - inline bool testCornerInTetrahedron(const TriCorner corner) const; - inline bool testCornerOnXYZFacet(const TriCorner corner) const; - inline bool testCornerAboveXYZFacet(const TriCorner corner) const; // ---------------------------------------------------------------------------------- // Utility methods used in intersection tests // ---------------------------------------------------------------------------------- - bool testTriangleSurroundsEdge(const TetraEdge edge) const; - inline bool testEdgeIntersectsTriangle(const TetraEdge edge) const; - inline bool testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const; - inline bool testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const; - bool testSegmentIntersectsHPlane(const TriSegment seg) const; - bool testSurfaceAboveCorner(const TetraCorner corner) const; - bool testTriangleSurroundsRay(const TetraCorner corner) const; // ---------------------------------------------------------------------------------- // Double and triple product calculations // ---------------------------------------------------------------------------------- - - void resetNearZeroCoordinates(); - + void handleDegenerateCases(); bool areDoubleProductsConsistent(const TriSegment seg) const; - - void preCalculateDoubleProducts(void); - + void preCalculateDoubleProducts(); inline void resetDoubleProducts(const TriSegment seg, const TetraCorner corner); - double calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const; - - void preCalculateTripleProducts(void); - + void preCalculateTripleProducts(); double calculateAngleEdgeTriangle(const TetraEdge edge) const; - inline double calcStableC(const TriSegment seg, const DoubleProduct dp) const; - inline double calcStableT(const TetraCorner corner) const; + inline double calcUnstableC(const TriSegment seg, const DoubleProduct dp, double & delta) const; + double calcTByDevelopingRow(const TetraCorner corner, const int row, const bool project) const; - inline double calcUnstableC(const TriSegment seg, const DoubleProduct dp) const; - - double calcTByDevelopingRow(const TetraCorner corner, const int row = 1, const bool project = false) const; + // ---------------------------------------------------------------------------------- + // Debug + // ---------------------------------------------------------------------------------- + inline const std::string& strTC(TetraCorner tc) const; + inline const std::string& strTE(TetraEdge te) const; + inline const std::string& strTF(TetraFacet tf) const; + inline const std::string& strTriC(TriCorner tc) const; + inline const std::string& strTriS(TriSegment tc) const; // ---------------------------------------------------------------------------------- // Member variables @@ -282,8 +245,12 @@ namespace INTERP_KERNEL /// following order in enumeration DoubleProduct double _doubleProducts[24]; + double _deltas[24]; + /// Array containing the 4 triple products. /// order : t_O, t_X, t_Y, t_Z + /// For example t_O represent the signed volume of the tetrahedron OPQR, and is positive if PQR is oriented clockwise + // when seen from the vertex O. double _tripleProducts[4]; /// Vector holding the points of the intersection polygon A. @@ -340,11 +307,11 @@ namespace INTERP_KERNEL // by a given row static const DoubleProduct DP_FOR_DETERMINANT_EXPANSION[12]; - // values used to decide how imprecise the double products + // values used to decide how/when imprecise the double products // should be to set them to 0.0 - static const long double MACH_EPS; // machine epsilon - static const long double MULT_PREC_F; // precision of multiplications (Grandy : f) - static const long double THRESHOLD_F; // threshold for zeroing (Grandy : F/f) + static const double MACH_EPS; // machine epsilon + static const double MULT_PREC_F; // precision of multiplications (Grandy : f) + static const double THRESHOLD_F; // threshold for zeroing (Grandy : F/f) static const double TRIPLE_PRODUCT_ANGLE_THRESHOLD; @@ -354,10 +321,10 @@ namespace INTERP_KERNEL // signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION static const double SIGN_FOR_SEG_FACET_INTERSECTION[12]; - + // coordinates of corners of tetrahedron static const double COORDS_TET_CORNER[12]; - + // indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION // for the calculation of the coordinates (x,y,z) of the intersection points // for Segment-Facet and Segment-Edge intersections @@ -424,7 +391,7 @@ namespace INTERP_KERNEL return _tripleProducts[corner]; } - inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp) const + inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp, double& delta) const { // find the points of the triangle // 0 -> P, 1 -> Q, 2 -> R @@ -435,7 +402,10 @@ namespace INTERP_KERNEL const int off1 = DP_OFFSET_1[dp]; const int off2 = DP_OFFSET_2[dp]; - return _coords[5*pt1 + off1] * _coords[5*pt2 + off2] - _coords[5*pt1 + off2] * _coords[5*pt2 + off1]; + const double prd1 = _coords[5*pt1 + off1] * _coords[5*pt2 + off2], + prd2 = _coords[5*pt1 + off2] * _coords[5*pt2 + off1]; + delta = std::fabs(prd1) + std::fabs(prd2); + return prd1 - prd2; } // ---------------------------------------------------------------------------------- @@ -517,8 +487,7 @@ namespace INTERP_KERNEL inline bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const { - // correspondence edge - triple products - // for edges OX, ..., ZX (Grandy, table III) + // correspondence edge - triple products for edges OX, ..., ZX (Grandy, table III) static const TetraCorner TRIPLE_PRODUCTS[12] = { X, O, // OX @@ -533,9 +502,17 @@ namespace INTERP_KERNEL const double t1 = calcStableT(TRIPLE_PRODUCTS[2*edge]); const double t2 = calcStableT(TRIPLE_PRODUCTS[2*edge + 1]); - //? should equality with zero use epsilon? - LOG(5, "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 ); - return (t1*t2 <= 0.0) && !epsilonEqual(t1 - t2, 0.0); // tuleap26461 + // [ABN] Okayyy: if either t1 or t2 exactly equal zero, then it can mean two things: + // - either PQR is very close to the corner -> this is OK, further computation of intersection point between + // surface and edge will produce a correct result + // - or, if the other triple prod is also very small, then this is a degenerate case: the edge is almost in PQR -> this is bad + // and leads to weird intersection point computation -> we avoid this. + // PS : here was written "// tuleap26461" -> whoo this helps :-) + if (t1 == 0.0 || t2 == 0.0) + if (std::fabs(t1+t2) < THRESHOLD_F*MULT_PREC_F) + return false; + + return (t1*t2 <= 0.0) && !epsilonEqual(t1,t2, (double)MULT_PREC_F); } inline bool TransformedTriangle::testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const @@ -583,19 +560,20 @@ namespace INTERP_KERNEL inline bool TransformedTriangle::testSurfaceAboveCorner(const TetraCorner corner) const { - // ? There seems to be an error in Grandy -> it should be C_XY instead of C_YZ in [28]. - // ? I haven't really figured out why, but it seems to work. - const double normal = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY); - - LOG(6, "surface above corner " << corner << " : " << "n = " << normal << ", t = [" << calcTByDevelopingRow(corner, 1, false) << ", " << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) ); - LOG(6, "] - stable : " << calcStableT(corner) ); - - //? we don't care here if the triple product is "invalid", that is, the triangle does not surround one of the - // edges going out from the corner (Grandy [53]) - if(!_validTP[corner]) - return ( calcTByDevelopingRow(corner, 1, false) * normal ) >= 0.0; - else - return ( calcStableT(corner) * normal ) >= 0.0; + // There is an error in Grandy -> it should be C_XY instead of C_YZ in [28]. + // + // Idea: the nz value (Grandy [28] corrected!) can be interpreted as a special variant of the triple product t_O + // where the z coordinate has been set to 1. It represents the signed volume of the tet OP'Q'R' where P', Q' and R' are the + // projection of P, Q, R on the z=1 plane. + // Comparing the sign of this triple product with t_X (or t_Y, ... dep on the corner) indicates whether the corner is in the + // half-space above or below the PQR triangle, similarly to what is explained in [16]. + // (this works even for the Z corner, since we could have chosen z=24 (instead of z=1) this would not have changed the final sign test). + const double nz = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY); + + // Triple product might have not been computed, but here we need one: + const double tp = _validTP[corner] ? calcStableT(corner) : calcTByDevelopingRow(corner, 1, false); + + return tp*nz >= 0.0; } inline bool TransformedTriangle::testTriangleSurroundsRay(const TetraCorner corner) const @@ -651,5 +629,4 @@ namespace INTERP_KERNEL } - #endif diff --git a/src/INTERP_KERNEL/TransformedTriangleMath.cxx b/src/INTERP_KERNEL/TransformedTriangleMath.cxx index 6e972b66d..bac06d5cd 100644 --- a/src/INTERP_KERNEL/TransformedTriangleMath.cxx +++ b/src/INTERP_KERNEL/TransformedTriangleMath.cxx @@ -62,24 +62,55 @@ namespace INTERP_KERNEL }; /// The machine epsilon, used in precision corrections - const long double TransformedTriangle::MACH_EPS = std::numeric_limits::epsilon(); + const double TransformedTriangle::MACH_EPS = std::numeric_limits::epsilon(); /// 4.0 * the machine epsilon, represents the precision of multiplication when performing corrections corrections ( f in Grandy ) - const long double TransformedTriangle::MULT_PREC_F = 4.0 * TransformedTriangle::MACH_EPS; + const double TransformedTriangle::MULT_PREC_F = 4.0 * TransformedTriangle::MACH_EPS; /// Threshold for resetting double and triple products to zero; ( F / f in Grandy ) - const long double TransformedTriangle::THRESHOLD_F = 500.0; + const double TransformedTriangle::THRESHOLD_F = 20.0; /// Threshold for what is considered a small enough angle to warrant correction of triple products by Grandy, [57] const double TransformedTriangle::TRIPLE_PRODUCT_ANGLE_THRESHOLD = 0.1; - // after transformation to the U-space, coordinates are inaccurate - // small variations around zero should not be taken into account - void TransformedTriangle::resetNearZeroCoordinates() + // Handle cases where one of the segment (or all) is (almost) in XYZ plane. + // We follow Grandy's suggestion and perturb slightly to have exactly h=0 for the segment (Grandy p.447) + // Note that if PQR is == to the upper facet of the unit tetra (XYZ), the tetra-corner-inclusion test should take it in, + // thanks to Grandy [21] and the fact that S_x test is "<=0" (not <0) + // After that, we also snap P,Q,R to the corners if they're very close. + void TransformedTriangle::handleDegenerateCases() { - for (int i=0; i<15; i++) - if (fabs(_coords[i]) distances; @@ -120,33 +154,21 @@ namespace INTERP_KERNEL distances.insert( std::make_pair( dist, corner ) ); } - // first element -> minimum distance + // first element -> minimum distance (because map is sorted) const TetraCorner minCorner = distances.begin()->second; resetDoubleProducts(seg, minCorner); distances.clear(); } } - + // -- (2) check that each double product satisfies Grandy, [47], else set to 0 for(TriSegment seg = PQ ; seg <= RP ; seg = TriSegment(seg + 1)) { for(DoubleProduct dp = C_YZ ; dp <= C_10 ; dp = DoubleProduct(dp + 1)) { - // find the points of the triangle - // 0 -> P, 1 -> Q, 2 -> R - const int pt1 = seg; - const int pt2 = (seg + 1) % 3; - - // find offsets - const int off1 = DP_OFFSET_1[dp]; - const int off2 = DP_OFFSET_2[dp]; + const int idx = 8*seg+dp; - const double term1 = _coords[5*pt1 + off1] * _coords[5*pt2 + off2]; - const double term2 = _coords[5*pt1 + off2] * _coords[5*pt2 + off1]; - - const long double delta = MULT_PREC_F * ( std::fabs(term1) + std::fabs(term2) ); - - if( epsilonEqual(_doubleProducts[8*seg + dp], 0.0, (double)(THRESHOLD_F * delta))) + if( epsilonEqual(_doubleProducts[idx], 0.0, THRESHOLD_F * MULT_PREC_F * _deltas[idx])) { // debug output #if LOG_LEVEL >= 5 @@ -157,13 +179,11 @@ namespace INTERP_KERNEL } #endif - - _doubleProducts[8*seg + dp] = 0.0; - + _doubleProducts[idx] = 0.0; } } } - + _is_double_products_calculated = true; } @@ -176,30 +196,45 @@ namespace INTERP_KERNEL */ bool TransformedTriangle::areDoubleProductsConsistent(const TriSegment seg) const { - const double term1 = _doubleProducts[8*seg + C_YZ] * _doubleProducts[8*seg + C_XH]; - const double term2 = _doubleProducts[8*seg + C_ZX] * _doubleProducts[8*seg + C_YH]; - const double term3 = _doubleProducts[8*seg + C_XY] * _doubleProducts[8*seg + C_ZH]; + // Careful! Here doubleProducts have not yet been corrected for roundoff errors! + // So we need to epsilon-adjust to correctly identify zeros: + static const DoubleProduct DP_LST[6] = {C_YZ, C_XH, + C_ZX, C_YH, + C_XY, C_ZH}; + double dps[6]; + for (int i = 0; i < 6; i++) + { + const double dp = _doubleProducts[8*seg + DP_LST[i]]; + dps[i] = dp; + } + + const double term1 = dps[0] * dps[1]; + const double term2 = dps[2] * dps[3]; + const double term3 = dps[4] * dps[5]; LOG(2, "for seg " << seg << " consistency " << term1 + term2 + term3 ); LOG(2, "term1 :" << term1 << " term2 :" << term2 << " term3: " << term3 ); + // Test for "== 0.0" here is OK since doubleProduct has been fixed for rounding to zero already. const int num_zero = (term1 == 0.0 ? 1 : 0) + (term2 == 0.0 ? 1 : 0) + (term3 == 0.0 ? 1 : 0); const int num_neg = (term1 < 0.0 ? 1 : 0) + (term2 < 0.0 ? 1 : 0) + (term3 < 0.0 ? 1 : 0); const int num_pos = (term1 > 0.0 ? 1 : 0) + (term2 > 0.0 ? 1 : 0) + (term3 > 0.0 ? 1 : 0); assert( num_zero + num_neg + num_pos == 3 ); - // calculated geometry is inconsistent if we have one of the following cases + // Calculated geometry is inconsistent if we have one of the following cases // * one term zero and the other two of the same sign // * two terms zero // * all terms positive // * all terms negative - if(((num_zero == 1 && num_neg != 1) || num_zero == 2 || (num_neg == 0 && num_zero != 3) || num_neg == 3 )) - { + const bool inconsist = (num_zero == 1 && num_neg != 1) || + num_zero == 2 || + (num_neg == 0 && num_zero != 3) || + num_neg == 3; + if(inconsist) { LOG(4, "inconsistent dp found" ); } - return !((num_zero == 1 && num_neg != 1) || num_zero == 2 || (num_neg == 0 && num_zero != 3) || num_neg == 3 ); - + return !inconsist; } /** @@ -250,12 +285,10 @@ namespace INTERP_KERNEL * the problem of errors due to cancellation. * */ - void TransformedTriangle::preCalculateTripleProducts(void) + void TransformedTriangle::preCalculateTripleProducts() { if(_is_triple_products_calculated) - { - return; - } + return; // find edge / row to use -> that whose edge makes the smallest angle to the triangle // use a map to find the minimum @@ -272,7 +305,7 @@ namespace INTERP_KERNEL // get edge by using correspondence between Double Product and Edge TetraEdge edge = TetraEdge(dp); - + // use edge only if it is surrounded by the surface if( _triangleSurroundsEdgeCache[edge] ) { @@ -288,13 +321,10 @@ namespace INTERP_KERNEL const int minRow = anglesForRows.begin()->second; if(minAngle < TRIPLE_PRODUCT_ANGLE_THRESHOLD) - { - _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, true); - } + _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, true); else - { - _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, false); - } + _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, false); + _validTP[corner] = true; } else @@ -303,12 +333,9 @@ namespace INTERP_KERNEL LOG(6, "Triple product not calculated for corner " << corner ); _tripleProducts[corner] = -3.14159265; _validTP[corner] = false; - } anglesForRows.clear(); - } - _is_triple_products_calculated = true; } @@ -371,12 +398,15 @@ namespace INTERP_KERNEL } /** - * Calculates triple product associated with the given corner of tetrahedron, developing + * Calculates triple product associated with the given corner of tetrahedron, developing * the determinant by the given row. The triple product gives the signed volume of * the tetrahedron between this corner and the triangle PQR. If the flag project is true, * one coordinate is projected out in order to eliminate errors in the intersection point * calculation due to cancellation. * + * Consistency with the double product computation and potential cancellation is also done here. + * + * * @pre double products have already been calculated * @param corner corner for which the triple product is calculated * @param row row (1 <= row <= 3) used to calculate the determinant @@ -439,23 +469,23 @@ namespace INTERP_KERNEL const double cPQ = calcStableC(PQ, dp); double alpha = 0.0; - + // coordinate to use for projection (Grandy, [57]) with edges // OX, OY, OZ, XY, YZ, ZX in order : // (y, z, x, h, h, h) // for the first three we could also use {2, 0, 1} static const int PROJECTION_COORDS[6] = { 1, 2, 0, 3, 3, 3 } ; - + const int coord = PROJECTION_COORDS[ dp ]; - + // coordinate values for P, Q and R const double coordValues[3] = { _coords[5*P + coord], _coords[5*Q + coord], _coords[5*R + coord] }; - + if(project) { // products coordinate values with corresponding double product const double coordDPProd[3] = { coordValues[0] * cQR, coordValues[1] * cRP, coordValues[2] * cPQ }; - + const double sumDPProd = coordDPProd[0] + coordDPProd[1] + coordDPProd[2]; const double sumDPProdSq = dot(coordDPProd, coordDPProd); @@ -467,6 +497,22 @@ namespace INTERP_KERNEL const double cRPbar = cRP * (1.0 - alpha * coordValues[1] * cRP); const double cPQbar = cPQ * (1.0 - alpha * coordValues[2] * cPQ); + // [ABN] Triple product cancellation logic: + // This part is not well described in Grandy (end of p.446) : + // "We use a method analogous to (47) to remove imprecise triple products,..." + // + // Our algo for cancelling a triple product: + // - retrieve the deltas associated with each DP involved (because a DP itself is a sum of two terms - see [42] + // - multiply them by the coordinate coming from the determinant expansion + // - and finally sum the 3 corresponding terms of the developement + // + // Using directly the DP (as was done before here) leads to issues, since this DP might have been cancelled + // already earlier on, and we lost the delta information -> doing this, we risk not cancelling the triple prod + // when we should have. + const double cQRbar_delta = _deltas[8*QR + dp], + cRPbar_delta = _deltas[8*RP + dp], + cPQbar_delta = _deltas[8*PQ + dp]; + // check sign of barred products - should not change // assert(cQRbar * cQR >= 0.0); //assert(cRPbar * cRP >= 0.0); @@ -476,15 +522,13 @@ namespace INTERP_KERNEL const double q_term = _coords[5*Q + offset] * cRPbar; const double r_term = _coords[5*R + offset] * cPQbar; - // check that we are not so close to zero that numerical errors could take over, - // otherwise reset to zero (cf Grandy, p. 446) -#ifdef FIXED_DELTA - const double delta = FIXED_DELTA; -#else - const long double delta = MULT_PREC_F * (std::fabs(p_term) + std::fabs(q_term) + std::fabs(r_term)); -#endif + const double p_delta = std::fabs(_coords[5*P + offset] * cQRbar_delta), + q_delta = std::fabs(_coords[5*Q + offset] * cRPbar_delta), + r_delta = std::fabs(_coords[5*R + offset] * cPQbar_delta); + + const double delta = p_delta + q_delta + r_delta; - if( epsilonEqual( p_term + q_term + r_term, 0.0, (double)(THRESHOLD_F * delta)) ) + if( epsilonEqual( p_term + q_term + r_term, 0.0, THRESHOLD_F * MULT_PREC_F * delta) ) { LOG(4, "Reset imprecise triple product for corner " << corner << " to zero" ); return 0.0; diff --git a/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.cxx b/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.cxx index ca6cbb6bd..41c12e1c7 100644 --- a/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.cxx +++ b/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.cxx @@ -19,6 +19,7 @@ #include "TransformedTriangleIntersectTest.hxx" #include +#include #include "Log.hxx" @@ -1194,7 +1195,7 @@ namespace INTERP_TEST CPPUNIT_ASSERT_EQUAL(false, tri->testSurfaceEdgeIntersection(TT::XY)); // surface-ray (3 possibilities) - CPPUNIT_ASSERT_EQUAL(true , tri->testSurfaceRayIntersection(TT::X)); + CPPUNIT_ASSERT_EQUAL(true, tri->testSurfaceRayIntersection(TT::X)); CPPUNIT_ASSERT_EQUAL(false, tri->testSurfaceRayIntersection(TT::Y)); CPPUNIT_ASSERT_EQUAL(false, tri->testSurfaceRayIntersection(TT::Z)); @@ -2127,10 +2128,229 @@ namespace INTERP_TEST delete tri; } + // Extract from BoxHexa1.med (or was it BoxHexa2.med ?) + void TransformedTriangleIntersectTest::testTriangle_vol1() + { + LOG(1, "+++++++ Testing testTriangle_vol1" ); -} // NAMESPACE + double coords[9] = { + -0.8088235294117645, 0, 0.55882352941176472, + 0.44117647058823506, 0, 0.55882352941176483, + -0.89215686274509864, 1.3333333333333339, 0.55882352941176483}; + + double refVol = 0.054383777732546296; + + TransformedTriangle tri(&coords[0], &coords[3], &coords[6]); + const double vol = tri.calculateIntersectionVolume(); + + LOG(3, " --- Final points in polygon A"); + for(const auto& pt: tri._polygonA) + LOG(3,vToStr(pt)); + LOG(3, " --- Final points in polygon B"); + for(const auto& pt: tri._polygonB) + LOG(3,vToStr(pt)); + + CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10); + } + + // Extract from BoxHexa1.med (or was it BoxHexa2.med ?) + void TransformedTriangleIntersectTest::testTriangle_vol2() + { + LOG(1, "+++++++ Testing testTriangle_vol2" ); + + double coords[9] = { + 0.44117647058823506, 0, 0.55882352941176483, + -0.55882352941176472, 0, 1.5588235294117649, + -0.89215686274509864, 1.3333333333333339, 0.55882352941176483 }; + + double refVol = -0.06869529818848; + + TransformedTriangle tri(&coords[0], &coords[3], &coords[6]); + const double vol = tri.calculateIntersectionVolume(); + + LOG(3, " --- Final points in polygon A"); + for(const auto& pt: tri._polygonA) + LOG(3,vToStr(pt)); + LOG(3, " --- Final points in polygon B"); + for(const auto& pt: tri._polygonB) + LOG(3,vToStr(pt)); + + CPPUNIT_ASSERT(tri.isTriangleInPlaneOfFacet(TransformedTriangle::XYZ)); + CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10); + } + + // Extract from BoxModSmall1.med + void TransformedTriangleIntersectTest::testTriangle_vol3() + { + LOG(1, "+++++++ Testing testTriangle_vol3" ); + + double coords[9] = { + 0, -4.4408920985006262e-16, 1, + -1.2062474433365091, -0.037350951323461778, 2.1879983126221099, + 0.49877186496532655, 0.59827776894780405, 0.79353793765518521 + }; + double refVol = -0.051135429735185; + + TransformedTriangle tri(&coords[0], &coords[3], &coords[6]); + const double vol = tri.calculateIntersectionVolume(); + + LOG(3, " --- Final points in polygon A"); + for(const auto& pt: tri._polygonA) + LOG(3,vToStr(pt)); + LOG(3, " --- Final points in polygon B"); + for(const auto& pt: tri._polygonB) + LOG(3,vToStr(pt)); + + CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10); + } + + // Taken from Box1Moderate.med + void TransformedTriangleIntersectTest::testTriangle_vol4() + { + LOG(1, "+++++++ Testing testTriangle_vol4" ); + + double coords[9] = { + 1, 3.552713678800501e-15, 0, + 2.022774182629973, -1.020222639063029, -0.01375178680446254, + 0.7495960843059706, 0.1125313911637846, 0.7430770879625861 + }; + double refVol = -0.00060846166394417; + + TransformedTriangle tri(&coords[0], &coords[3], &coords[6]); + const double vol = tri.calculateIntersectionVolume(); + + LOG(3, " --- Final points in polygon A"); + for(const auto& pt: tri._polygonA) + LOG(3,vToStr(pt)); + LOG(3, " --- Final points in polygon B"); + for(const auto& pt: tri._polygonB) + LOG(3,vToStr(pt)); + CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10); + } + + // Taken from Box1Moderate.med too + void TransformedTriangleIntersectTest::testTriangle_vol5() + { + LOG(1, "+++++++ Testing testTriangle_vol5" ); + + double coords[9] = { + 1, 3.552713678800501e-15, 0, + -3.552713678800501e-15, 0, 0.9999999999999982, + 0, 1.000000000000004, -8.881784197001252e-16 + }; + double refVol = -1/6.; + + TransformedTriangle tri(&coords[0], &coords[3], &coords[6]); + const double vol = tri.calculateIntersectionVolume(); + LOG(3, " --- Final points in polygon A"); + for(const auto& pt: tri._polygonA) + LOG(3,vToStr(pt)); + LOG(3, " --- Final points in polygon B"); + for(const auto& pt: tri._polygonB) + LOG(3,vToStr(pt)); + + CPPUNIT_ASSERT(tri.isTriangleInPlaneOfFacet(TransformedTriangle::XYZ)); + CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10); + } + + // Taken from Box1Moderate.med again + void TransformedTriangleIntersectTest::testTriangle_vol6() + { + LOG(1, "+++++++ Testing testTriangle_vol6" ); + + double coords[9] = { + 1.000000000000004, 0, 0, + 0, 0, 0.9999999999999929, + 3.552713678800501e-15, 1, 0}; + double refVol = -1/6.; + + TransformedTriangle tri(&coords[0], &coords[3], &coords[6]); + const double vol = tri.calculateIntersectionVolume(); + + LOG(3, " --- Final points in polygon A"); + for(const auto& pt: tri._polygonA) + LOG(3,vToStr(pt)); + LOG(3, " --- Final points in polygon B"); + for(const auto& pt: tri._polygonB) + LOG(3,vToStr(pt)); + + CPPUNIT_ASSERT(tri.isTriangleInPlaneOfFacet(TransformedTriangle::XYZ)); + CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10); + } + + // Taken from MEDCouplingRemapperTest.py::testCellToNodeReverse3D() + void TransformedTriangleIntersectTest::testTriangle_vol7() + { + LOG(1, "+++++++ Testing testTriangle_vol7" ); + + double coords[9] = { + 0.9999999999999999, 3.700743415417188e-17, 3.700743415417188e-17, + -5.551115123125783e-17, 0, 1, + 3.700743415417188e-17, 0.9999999999999999, 3.700743415417188e-17 + }; + + double refVol = -1/6.; + + TransformedTriangle tri(&coords[0], &coords[3], &coords[6]); + const double vol = tri.calculateIntersectionVolume(); + + LOG(3, " --- Final points in polygon A"); + for(const auto& pt: tri._polygonA) + LOG(3,vToStr(pt)); + LOG(3, " --- Final points in polygon B"); + for(const auto& pt: tri._polygonB) + LOG(3,vToStr(pt)); + + CPPUNIT_ASSERT(tri.isTriangleInPlaneOfFacet(TransformedTriangle::XYZ)); + CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10); + + // Polygon A should not contain anything else than tetra corners: + double barycenter[3]; + tri.calculatePolygonBarycenter(TransformedTriangle::A, barycenter); + tri.sortIntersectionPolygon(TransformedTriangle::A, barycenter); + const double eps = TransformedTriangle::THRESHOLD_F * TransformedTriangle::MULT_PREC_F; + for (auto& p: tri._polygonA) + for (int d=0; d<3; d++) + { + CPPUNIT_ASSERT(fabs(p[d]) < eps || fabs(p[d]-1) < eps); + } + } + + // This test is useful if TransformedTriangle::THRESHOLD_F = 20 (like in Grandy's article) + // With this value, it is an almost degenerate case : the triangle is very close to XYZ plane, but this still + // falls in the 'non-degenerate' category. Total volume is 1/6 but coming from + // contributions from A and B. + // With current value of threshold (=100, chosen because of P1P1 intersector), triangle becomes completely degenerated. + void TransformedTriangleIntersectTest::testTriangle_vol8() + { + LOG(1, "+++++++ Testing testTriangle_vol8" ); + + double coords[9] = { + 0.9999999999999787, 7.105427357601002e-15, -3.552713678800501e-15, + -1.4210854715202e-14, 0, 0.9999999999999964, + -7.105427357601002e-15, 1.000000000000014, -3.552713678800501e-15 + }; + + double refVol = -1/6.; + + TransformedTriangle tri(&coords[0], &coords[3], &coords[6]); + const double vol = tri.calculateIntersectionVolume(); + + LOG(3, " --- Final points in polygon A"); + for(const auto& pt: tri._polygonA) + LOG(3,vToStr(pt)); + LOG(3, " --- Final points in polygon B"); + for(const auto& pt: tri._polygonB) + LOG(3,vToStr(pt)); + +// CPPUNIT_ASSERT (! tri.isTriangleInPlaneOfFacet(TransformedTriangle::XYZ) ); // not degenerate if threshold == 20.0 + CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10); + } + + +} // NAMESPACE diff --git a/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.hxx b/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.hxx index 1fa23e6f0..8d7ef5e64 100644 --- a/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.hxx +++ b/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.hxx @@ -47,6 +47,17 @@ namespace INTERP_TEST CPPUNIT_TEST( testTriangle12 ); CPPUNIT_TEST( testTriangle13 ); + + // Tests for degenerated cases where PQR is almost in XYZ plane: + CPPUNIT_TEST( testTriangle_vol1 ); + CPPUNIT_TEST( testTriangle_vol2 ); + CPPUNIT_TEST( testTriangle_vol3 ); + CPPUNIT_TEST( testTriangle_vol4 ); + CPPUNIT_TEST( testTriangle_vol5 ); + CPPUNIT_TEST( testTriangle_vol6 ); + CPPUNIT_TEST( testTriangle_vol7 ); + CPPUNIT_TEST( testTriangle_vol8 ); + CPPUNIT_TEST_SUITE_END(); typedef INTERP_KERNEL::TransformedTriangle::TriSegment TriSegment; @@ -55,40 +66,30 @@ namespace INTERP_TEST public: void testTriangle1(); - void testTriangle2(); - void testTriangle3(); - void testTriangle4(); - void testTriangle5(); - void testTriangle6(); - void testTriangle7(); - void testTriangle8(); - void testTriangle9(); - void testTriangle10(); - void testTriangle11(); - void testTriangle12(); - void testTriangle13(); - private: - + void testTriangle_vol1(); + void testTriangle_vol2(); + void testTriangle_vol3(); + void testTriangle_vol4(); + void testTriangle_vol5(); + void testTriangle_vol6(); + void testTriangle_vol7(); + void testTriangle_vol8(); + }; } - - - - - #endif diff --git a/src/INTERP_KERNELTest/TransformedTriangleTest.cxx b/src/INTERP_KERNELTest/TransformedTriangleTest.cxx index bfe66cb29..856e06086 100644 --- a/src/INTERP_KERNELTest/TransformedTriangleTest.cxx +++ b/src/INTERP_KERNELTest/TransformedTriangleTest.cxx @@ -148,15 +148,15 @@ namespace INTERP_TEST double c_vals[3 * 8]; for(TriSegment seg = TransformedTriangle::PQ ; seg <= TransformedTriangle::RP ; seg = TriSegment(seg + 1)) { - - c_vals[seg*8 + 0] = tri1->calcUnstableC(seg, TransformedTriangle::C_XY); - c_vals[seg*8 + 1] = tri1->calcUnstableC(seg, TransformedTriangle::C_YZ); - c_vals[seg*8 + 2] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZX); - c_vals[seg*8 + 3] = tri1->calcUnstableC(seg, TransformedTriangle::C_XH); - c_vals[seg*8 + 4] = tri1->calcUnstableC(seg, TransformedTriangle::C_YH); - c_vals[seg*8 + 5] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZH); - c_vals[seg*8 + 6] = tri1->calcUnstableC(seg, TransformedTriangle::C_01); - c_vals[seg*8 + 7] = tri1->calcUnstableC(seg, TransformedTriangle::C_10); + double dnu; + c_vals[seg*8 + 0] = tri1->calcUnstableC(seg, TransformedTriangle::C_XY, dnu); + c_vals[seg*8 + 1] = tri1->calcUnstableC(seg, TransformedTriangle::C_YZ, dnu); + c_vals[seg*8 + 2] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZX, dnu); + c_vals[seg*8 + 3] = tri1->calcUnstableC(seg, TransformedTriangle::C_XH, dnu); + c_vals[seg*8 + 4] = tri1->calcUnstableC(seg, TransformedTriangle::C_YH, dnu); + c_vals[seg*8 + 5] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZH, dnu); + c_vals[seg*8 + 6] = tri1->calcUnstableC(seg, TransformedTriangle::C_01, dnu); + c_vals[seg*8 + 7] = tri1->calcUnstableC(seg, TransformedTriangle::C_10, dnu); } @@ -253,12 +253,13 @@ namespace INTERP_TEST // find unstable products to check for consistency (Grandy [46]) for(TriSegment seg = TransformedTriangle::PQ ; seg <= TransformedTriangle::RP ; seg = TriSegment(seg + 1)) { - const double c_xy = tri2->calcUnstableC(seg, TransformedTriangle::C_XY); - const double c_yz = tri2->calcUnstableC(seg, TransformedTriangle::C_YZ); - const double c_zx = tri2->calcUnstableC(seg, TransformedTriangle::C_ZX); - const double c_xh = tri2->calcUnstableC(seg, TransformedTriangle::C_XH); - const double c_yh = tri2->calcUnstableC(seg, TransformedTriangle::C_YH); - const double c_zh = tri2->calcUnstableC(seg, TransformedTriangle::C_ZH); + double dnu; + const double c_xy = tri2->calcUnstableC(seg, TransformedTriangle::C_XY, dnu); + const double c_yz = tri2->calcUnstableC(seg, TransformedTriangle::C_YZ, dnu); + const double c_zx = tri2->calcUnstableC(seg, TransformedTriangle::C_ZX, dnu); + const double c_xh = tri2->calcUnstableC(seg, TransformedTriangle::C_XH, dnu); + const double c_yh = tri2->calcUnstableC(seg, TransformedTriangle::C_YH, dnu); + const double c_zh = tri2->calcUnstableC(seg, TransformedTriangle::C_ZH, dnu); const int num_zero = (c_yz*c_xh == 0.0 ? 1 : 0) + (c_zx*c_yh == 0.0 ? 1 : 0) + (c_xy*c_zh == 0.0 ? 1 : 0); const int num_neg = (c_yz*c_xh < 0.0 ? 1 : 0) + (c_zx*c_yh < 0.0 ? 1 : 0) + (c_xy*c_zh < 0.0 ? 1 : 0); diff --git a/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx b/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx index 4ae487e87..991b93c76 100644 --- a/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx +++ b/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx @@ -300,6 +300,41 @@ namespace INTERP_TEST CPPUNIT_ASSERT_DOUBLES_EQUAL(6944.4444444444443,volume,1e-9); } + /** + * Extract from a single tetra from BoxHexa1.med and BoxHexa2.med. + * Symetry test was failing. + */ + void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_14() + { + double S[] = { + 25.0, 96.0, 0.0, + 25.0, 120.0, 0.0, + 37.5, 120.0, 0.0, + 25.0, 120.0, 11.333333333333339255}; + + double T[] = { + 25.0, 90.0, 6.333333333333339255, + 25.0, 120.0, 6.3333333333333357018, + 41.6666666666666714036, 120.0, 6.3333333333333348136, + 25.0, 120.0, 17.6666666666666714036}; + + mcIdType conn[4] = { 0,1,2,3 }; + const double* tnodes[4]={ T, T+3, T+6, T+9 }; + const double* snodes[4]={ S, S+3, S+6, S+9 }; + const double refVol = 48.6591695501729; + + __MESH_DUMMY dummyMesh; + SplitterTetra<__MESH_DUMMY> src( dummyMesh, snodes, conn ); + double volume = src.intersectTetra( tnodes ); + CPPUNIT_ASSERT_DOUBLES_EQUAL(refVol,volume,1e-9); + + // Now the other way round: + SplitterTetra<__MESH_DUMMY> tgt( dummyMesh, tnodes, conn ); + double volume2 = tgt.intersectTetra( snodes ); + CPPUNIT_ASSERT_DOUBLES_EQUAL(refVol,volume2,1e-9); + } + + void UnitTetraIntersectionBaryTest::test_TetraAffineTransform_reverseApply() { double nodes[12] = { -4.0, 9.0, 3.0, diff --git a/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.hxx b/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.hxx index 0f217bedb..ab7dcecbd 100644 --- a/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.hxx +++ b/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.hxx @@ -37,6 +37,7 @@ namespace INTERP_TEST class INTERPKERNELTEST_EXPORT UnitTetraIntersectionBaryTest : public CppUnit::TestFixture { CPPUNIT_TEST_SUITE( UnitTetraIntersectionBaryTest ); + CPPUNIT_TEST( test_UnitTetraIntersectionBary_14 ); CPPUNIT_TEST( test_UnitTetraIntersectionBary_13 ); CPPUNIT_TEST( test_UnitTetraIntersectionBary_12 ); CPPUNIT_TEST( test_UnitTetraIntersectionBary_1 ); @@ -69,6 +70,7 @@ namespace INTERP_TEST void test_UnitTetraIntersectionBary_11(); void test_UnitTetraIntersectionBary_12(); void test_UnitTetraIntersectionBary_13(); + void test_UnitTetraIntersectionBary_14(); void test_TetraAffineTransform_reverseApply(); void test_barycentric_coords(); void test_cuboid_mapped_coords_3D();