From: abn Date: Thu, 25 Jan 2024 19:37:39 +0000 (+0100) Subject: [TetraIntersec] Several bug fixes in Grandy's implementation X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=47ece939902ff5b427fa8bf54a490a80946149c1;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 back to 20 (not 500) as in Grandy + 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 + fixed surface-edge intersection test: triple product equality (and zeroing) must be checked with care + removed TriangleTransformedInline.hxx to merge with header (helping IDEs!) --- diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index 596cd90e4..b3478ae9c 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -127,7 +127,8 @@ namespace INTERP_KERNEL _coords[5*Q + 4] = 1 - q[0] - q[1]; _coords[5*R + 4] = 1 - r[0] - r[1]; - resetNearZeroCoordinates(); + // Handle degenerated case where PQR is (almost) inside XYZ plane + handlePQRInXZYPlane(); // initialise rest of data preCalculateDoubleProducts(); @@ -785,7 +786,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 841ec476f..1c7aa1fd2 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -69,6 +69,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 @@ -94,18 +97,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 @@ -139,126 +138,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 handlePQRInXZYPlane(); 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 @@ -281,6 +244,8 @@ 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 double _tripleProducts[4]; @@ -339,11 +304,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; @@ -353,10 +318,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 +389,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 @@ -436,7 +401,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; } // ---------------------------------------------------------------------------------- @@ -632,8 +600,6 @@ namespace INTERP_KERNEL return ( cPQ*cQR > 0.0 ) && ( cPQ*cRP > 0.0 ); } - } - #endif diff --git a/src/INTERP_KERNEL/TransformedTriangleMath.cxx b/src/INTERP_KERNEL/TransformedTriangleMath.cxx index 72f42d1ba..a2bbf2ff1 100644 --- a/src/INTERP_KERNEL/TransformedTriangleMath.cxx +++ b/src/INTERP_KERNEL/TransformedTriangleMath.cxx @@ -62,24 +62,39 @@ 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 PQR is (almost) in XYZ plane. + // We follow Grandy's suggestion and perturb slightly to have exactly h=0 when h is very small (Grandy p.447) + // Note that if PQR == 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) + void TransformedTriangle::handlePQRInXZYPlane() { - for (int i=0; i<15; i++) - if (fabs(_coords[i]) THRESHOLD_F*TransformedTriangle::MULT_PREC_F) // Here, it is important to use same threshold as for DP computation! + cond = false; + + if (cond) // put the PQR triangle exactly inside the XYZ plane (h=0) + for (int i=0; i<3; i++) + { + const double correc = hPQR[i]/3.; // this should be really small! + _coords[i*5+0] += correc; + _coords[i*5+1] += correc; + _coords[i*5+2] += correc; + _coords[i*5+XYZ] = 0.0; // <=> -= h + } } // ---------------------------------------------------------------------------------- @@ -92,7 +107,7 @@ namespace INTERP_KERNEL * and it is thus the "stable" double products that are stored. * */ - void TransformedTriangle::preCalculateDoubleProducts(void) + void TransformedTriangle::preCalculateDoubleProducts() { if(_is_double_products_calculated) return; @@ -101,7 +116,10 @@ namespace INTERP_KERNEL for(TriSegment seg = PQ ; seg <= RP ; seg = TriSegment(seg + 1)) { for(DoubleProduct dp = C_YZ ; dp <= C_10 ; dp = DoubleProduct(dp + 1)) - _doubleProducts[8*seg + dp] = calcUnstableC(seg, dp); + { + const int idx = 8*seg + dp; + _doubleProducts[idx] = calcUnstableC(seg, dp, _deltas[idx]); + } } std::map distances; @@ -120,33 +138,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 +163,11 @@ namespace INTERP_KERNEL } #endif - - _doubleProducts[8*seg + dp] = 0.0; - + _doubleProducts[idx] = 0.0; } } } - + _is_double_products_calculated = true; } @@ -176,30 +180,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 +269,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 +289,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 +305,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 @@ -304,12 +318,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; } @@ -372,12 +383,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 @@ -440,23 +454,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); @@ -468,6 +482,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); @@ -477,15 +507,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/SingleElementTetraTests.hxx b/src/INTERP_KERNELTest/SingleElementTetraTests.hxx index 2d63a672a..dd974d994 100644 --- a/src/INTERP_KERNELTest/SingleElementTetraTests.hxx +++ b/src/INTERP_KERNELTest/SingleElementTetraTests.hxx @@ -33,6 +33,8 @@ namespace INTERP_TEST { CPPUNIT_TEST_SUITE( SingleElementTetraTests ); +// CPPUNIT_TEST( tetraAdri ); + CPPUNIT_TEST( tetraReflexiveUnit ); CPPUNIT_TEST( tetraReflexiveGeneral ); CPPUNIT_TEST( tetraNudgedSimpler ); @@ -53,6 +55,15 @@ namespace INTERP_TEST public: + void tetraAdri() + { +// _testTools->intersectMeshes("bh1_14", "bh2_14", 48.6591695501729); + + _testTools->intersectMeshes("bh2_14", "bh1_14", 48.6591695501729); + + } + + /// Unit tetrahedron mesh intersecting itself /// \brief Status : pass void tetraReflexiveUnit() diff --git a/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.cxx b/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.cxx index ca6cbb6bd..4c8328303 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" @@ -2128,9 +2129,160 @@ namespace INTERP_TEST } -} // NAMESPACE + // Extract from BoxHexa1.med (or was it BoxHexa2.med ?) + void TransformedTriangleIntersectTest::testTriangle_vol1() + { + LOG(1, "+++++++ Testing testTriangle_vol1" ); + + 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); + } + + +} // NAMESPACE diff --git a/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.hxx b/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.hxx index 1fa23e6f0..1aa358e3a 100644 --- a/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.hxx +++ b/src/INTERP_KERNELTest/TransformedTriangleIntersectTest.hxx @@ -47,6 +47,14 @@ 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_SUITE_END(); typedef INTERP_KERNEL::TransformedTriangle::TriSegment TriSegment; @@ -55,40 +63,28 @@ 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(); + }; } - - - - - #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();