From: abn Date: Mon, 5 Feb 2024 16:07:36 +0000 (+0100) Subject: [TetraIntersect] Formatting and including what's inline really inline! X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=88df74325e9fa59c9ad9d36f6f29bfad372e5beb;p=tools%2Fmedcoupling.git [TetraIntersect] Formatting and including what's inline really inline! (nothing new in this commit) + this helps IDEs to parse code ... + and various other minor formattings ... --- diff --git a/src/INTERP_KERNEL/Log.hxx b/src/INTERP_KERNEL/Log.hxx index 0c056c139..09e39744c 100644 --- a/src/INTERP_KERNEL/Log.hxx +++ b/src/INTERP_KERNEL/Log.hxx @@ -28,10 +28,10 @@ * LOG_LEVEL (which can be defined at compile-time for each file by passing option -DLOG_LEVEL=x to gcc) * than the message is logged. * - * - * */ + +//#define LOG_LEVEL 4 /// define LOG_LEVEL here if it is not already defined #ifndef LOG_LEVEL #define LOG_LEVEL 0 @@ -54,12 +54,4 @@ - - - - - - - - #endif diff --git a/src/INTERP_KERNEL/SplitterTetra.hxx b/src/INTERP_KERNEL/SplitterTetra.hxx index 9ba8f7b0d..80f704873 100644 --- a/src/INTERP_KERNEL/SplitterTetra.hxx +++ b/src/INTERP_KERNEL/SplitterTetra.hxx @@ -441,7 +441,7 @@ namespace INTERP_KERNEL * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace, * the corresponding array element will be true if and only if it was true when the method was called and pt lies in the halfspace. * - * @param pt double[3] containing the coordiantes of a transformed point + * @param pt double[3] containing the coordinates of a transformed point * @param isOutside bool[8] which indicate the results of earlier checks. */ template diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index 6ae42053b..9b631fd7b 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -216,12 +216,8 @@ namespace INTERP_KERNEL // halfspace filtering check // NB : might not be beneficial for caching of triangles for(int i = 0; i < 8; ++i) - { - if(isOutside[i]) - { - isTargetOutside = true; - } - } + if(isOutside[i]) + isTargetOutside = true; double totalVolume = 0.0; @@ -873,6 +869,7 @@ namespace INTERP_KERNEL TransformedTriangle tri(nodes[faceNodes[0]], nodes[faceNodes[1]], nodes[faceNodes[2]]); double vol = tri.calculateIntersectionVolume(); + LOG(1, "ii = " << ii << " Volume=" << vol) totalVolume += vol; } diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index fdaeb6562..596cd90e4 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -147,14 +147,10 @@ namespace INTERP_KERNEL TransformedTriangle::~TransformedTriangle() { // delete elements of polygons - for(std::vector::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it) - { - delete[] *it; - } - for(std::vector::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it) - { - delete[] *it; - } + for(auto& it: _polygonA) + delete[] it; + for(auto& it: _polygonB) + delete[] it; } /** @@ -206,6 +202,11 @@ namespace INTERP_KERNEL if(_polygonA.size() > 2) { LOG(2, "---- Treating polygon A ... "); +#if LOG_LEVEL > 0 + LOG(3, " --- Final points in polygon A"); + for(const auto& pt: _polygonA) + LOG(3,vToStr(pt)); +#endif calculatePolygonBarycenter(A, barycenter); sortIntersectionPolygon(A, barycenter); volA = calculateVolumeUnderPolygon(A, barycenter); @@ -214,17 +215,27 @@ namespace INTERP_KERNEL double volB = 0.0; // if triangle is not in h = 0 plane, calculate volume under B - if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ)) + if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ)) // second condition avoids double counting in case triangle fully included in h=0 facet { LOG(2, "---- Treating polygon B ... "); - +#if LOG_LEVEL > 0 + LOG(3, " --- Final points in polygon B"); + for(const auto& pt: _polygonB) + LOG(3,vToStr(pt)); +#endif calculatePolygonBarycenter(B, barycenter); sortIntersectionPolygon(B, barycenter); volB = calculateVolumeUnderPolygon(B, barycenter); LOG(2, "Volume is " << sign * volB); } - LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********"); +#if LOG_LEVEL >= 2 + LOG(2, "############ Triangle :") + dumpCoords(); + LOG(2, "vol A = " << volA); + LOG(2, "vol B = " << volB); + LOG(2, "TOTAL = " << sign*(volA+volB)); +#endif return _volume = sign * (volA + volB); @@ -265,7 +276,7 @@ namespace INTERP_KERNEL } // ---------------------------------------------------------------------------------- - // TransformedTriangle PRIVATE + // TransformedTriangle PROTECTED // ---------------------------------------------------------------------------------- /** @@ -278,6 +289,11 @@ namespace INTERP_KERNEL */ void TransformedTriangle::calculateIntersectionAndProjectionPolygons() { +#if LOG_LEVEL >= 3 + std::cout << " @@@@@@@@ COORDS @@@@@@ " << std::endl; + dumpCoords(); +#endif + assert(_polygonA.size() == 0); assert(_polygonB.size() == 0); // avoid reallocations in push_back() by pre-allocating enough memory @@ -293,15 +309,15 @@ namespace INTERP_KERNEL double* ptA = new double[3]; calcIntersectionPtSurfaceEdge(edge, ptA); _polygonA.push_back(ptA); - LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A "); + LOG(3,"Surface-edge (edge " << strTE(edge) << "): " << vToStr(ptA) << " added to A "); if(edge >= XY) { double* ptB = new double[3]; copyVector3(ptA, ptB); _polygonB.push_back(ptB); - LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B "); + LOG(3,"Surface-edge (edge " << strTE(edge) << "): " << vToStr(ptB) << " added to B "); } - + } } @@ -313,7 +329,7 @@ namespace INTERP_KERNEL double* ptB = new double[3]; copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); _polygonB.push_back(ptB); - LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B"); + LOG(3,"Surface-ray (corner " << strTC(corner) << "): " << vToStr(ptB) << " added to B"); } } @@ -323,33 +339,31 @@ namespace INTERP_KERNEL bool isZero[NO_DP]; - // check beforehand which double-products are zero + // check beforehand which double-products are zero. for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1)) - { - isZero[dp] = (calcStableC(seg, dp) == 0.0); - } + isZero[dp] = (calcStableC(seg, dp) == 0.0); // segment - facet for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1)) { // is this test worth it? const bool doTest = - !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] && - !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] && - !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]]; + !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] && + !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] && + !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]]; if(doTest && testSegmentFacetIntersection(seg, facet)) { double* ptA = new double[3]; calcIntersectionPtSegmentFacet(seg, facet, ptA); _polygonA.push_back(ptA); - LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A"); + LOG(3,"Segment-facet (facet " << strTF(facet) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A"); if(facet == XYZ) { double* ptB = new double[3]; copyVector3(ptA, ptB); _polygonB.push_back(ptB); - LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B"); + LOG(3,"Segment-facet (facet " << strTF(facet) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B"); } } @@ -365,112 +379,113 @@ namespace INTERP_KERNEL double* ptA = new double[3]; calcIntersectionPtSegmentEdge(seg, edge, ptA); _polygonA.push_back(ptA); - LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A"); + LOG(3,"Segment-edge (edge " << strTE(edge) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A"); if(edge >= XY) { double* ptB = new double[3]; copyVector3(ptA, ptB); _polygonB.push_back(ptB); + LOG(3,"Segment-edge (edge " << strTE(edge) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to B"); } } } - + // segment - corner for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1)) { const bool doTest = - isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] && - isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] && - isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )]; + isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] && + isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] && + isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )]; if(doTest && testSegmentCornerIntersection(seg, corner)) { double* ptA = new double[3]; copyVector3(&COORDS_TET_CORNER[3 * corner], ptA); _polygonA.push_back(ptA); - LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A"); + LOG(3,"Segment-corner (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A"); if(corner != O) { double* ptB = new double[3]; _polygonB.push_back(ptB); copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); - LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B"); + LOG(3,"Segment-corner (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B"); } } } - // segment - ray - for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1)) + // segment - ray + for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1)) + { + if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner)) { - if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner)) - { - double* ptB = new double[3]; - copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); - _polygonB.push_back(ptB); - LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B"); - } + double* ptB = new double[3]; + copyVector3(&COORDS_TET_CORNER[3 * corner], ptB); + _polygonB.push_back(ptB); + LOG(3,"Segment-ray (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B"); } - - // segment - halfstrip - for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1)) - { + } + + // segment - halfstrip + for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1)) + { #if 0 - const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX - const bool doTest = - !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] && - !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]]; - + const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX + const bool doTest = + !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] && + !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]]; + - if(doTest && testSegmentHalfstripIntersection(seg, edge)) + if(doTest && testSegmentHalfstripIntersection(seg, edge)) #endif - if(testSegmentHalfstripIntersection(seg, edge)) - { - double* ptB = new double[3]; - calcIntersectionPtSegmentHalfstrip(seg, edge, ptB); - _polygonB.push_back(ptB); - LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); - } - } + if(testSegmentHalfstripIntersection(seg, edge)) + { + double* ptB = new double[3]; + calcIntersectionPtSegmentHalfstrip(seg, edge, ptB); + _polygonB.push_back(ptB); + LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B"); + } + } } - // inclusion tests - for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1)) + // inclusion tests + for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1)) + { + // { XYZ - inclusion only possible if in Tetrahedron? + // tetrahedron + if(testCornerInTetrahedron(corner)) { - // { XYZ - inclusion only possible if in Tetrahedron? - // tetrahedron - if(testCornerInTetrahedron(corner)) - { - double* ptA = new double[3]; - copyVector3(&_coords[5*corner], ptA); - _polygonA.push_back(ptA); - LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A"); - } - - // XYZ - plane - if(testCornerOnXYZFacet(corner)) - { - double* ptB = new double[3]; - copyVector3(&_coords[5*corner], ptB); - _polygonB.push_back(ptB); - LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B"); - } + double* ptA = new double[3]; + copyVector3(&_coords[5*corner], ptA); + _polygonA.push_back(ptA); + LOG(3,"Inclusion tetrahedron (corner " << strTriC(corner) << "): " << vToStr(ptA) << " added to A"); + } - // projection on XYZ - facet - if(testCornerAboveXYZFacet(corner)) - { - double* ptB = new double[3]; - copyVector3(&_coords[5*corner], ptB); - ptB[2] = 1 - ptB[0] - ptB[1]; - assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0)); - _polygonB.push_back(ptB); - LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B"); - } + // XYZ - plane + if(testCornerOnXYZFacet(corner)) + { + double* ptB = new double[3]; + copyVector3(&_coords[5*corner], ptB); + _polygonB.push_back(ptB); + LOG(3,"Inclusion XYZ-plane (corner " << strTriC(corner) << "): " << vToStr(ptB) << " added to B"); + } + // projection on XYZ - facet + if(testCornerAboveXYZFacet(corner)) + { + double* ptB = new double[3]; + copyVector3(&_coords[5*corner], ptB); + ptB[2] = 1 - ptB[0] - ptB[1]; + assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0)); + _polygonB.push_back(ptB); + LOG(3,"Projection XYZ-plane (corner " << strTriC(corner) << "): " << vToStr(ptB) << " added to B"); } } + } + /** * Calculates the intersection polygon A, performing the intersection tests * and storing the corresponding point in the vector _polygonA. @@ -504,10 +519,9 @@ namespace INTERP_KERNEL bool isZero[NO_DP]; // check beforehand which double-products are zero + // Test for "== 0.0" here is OK since doubleProduct has been fixed for rounding to zero already. for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1)) - { - isZero[dp] = (calcStableC(seg, dp) == 0.0); - } + isZero[dp] = (calcStableC(seg, dp) == 0.0); // segment - facet for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1)) @@ -751,18 +765,13 @@ namespace INTERP_KERNEL */ bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) const { - // coordinate to check const int coord = static_cast(facet); for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1)) - { - if(_coords[5*c + coord] != 0.0) - { - return false; - } - } - + if(_coords[5*c + coord] != 0.0) + return false; + return true; } @@ -801,9 +810,7 @@ namespace INTERP_KERNEL double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2]; if(epsilonEqual(sign, 0.)) - { - sign = 0.; - } + sign = 0.; return (sign < 0.) ? -1 : (sign > 0.) ? 1 : 0; } @@ -815,13 +822,10 @@ namespace INTERP_KERNEL bool TransformedTriangle::isTriangleBelowTetraeder() const { for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1)) - { - // check z-coords for all points - if(_coords[5*c + 2] >= 0.0) - { - return false; - } - } + // check z-coords for all points + if(_coords[5*c + 2] >= 0.0) + return false; + return true; } @@ -833,10 +837,10 @@ namespace INTERP_KERNEL { std::cout << "Coords : "; for(int i = 0 ; i < 3; ++i) - { - std::cout << vToStr(&_coords[5*i]) << ","; - } + std::cout << vToStr(&_coords[5*i]) << ","; + std::cout << std::endl; } - + + } // NAMESPACE diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 2904d95d8..841ec476f 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -380,9 +380,259 @@ namespace INTERP_KERNEL }; - // include definitions of inline methods + inline void TransformedTriangle::preCalculateTriangleSurroundsEdge() + { + for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1)) + { + _triangleSurroundsEdgeCache[edge] = testTriangleSurroundsEdge(edge); + } + } + + + // ---------------------------------------------------------------------------------- + // TransformedTriangle_math.cxx + // ---------------------------------------------------------------------------------- + + inline void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner) + { + // set the three corresponding double products to 0.0 + static const DoubleProduct DOUBLE_PRODUCTS[12] = + { + C_YZ, C_ZX, C_XY, // O + C_YZ, C_ZH, C_YH, // X + C_ZX, C_ZH, C_XH, // Y + C_XY, C_YH, C_XH // Z + }; + + for(int i = 0 ; i < 3 ; ++i) { + const DoubleProduct dp = DOUBLE_PRODUCTS[3*corner + i]; + + LOG(6, std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner); + _doubleProducts[8*seg + dp] = 0.0; + }; + } + + inline double TransformedTriangle::calcStableC(const TriSegment seg, const DoubleProduct dp) const + { + return _doubleProducts[8*seg + dp]; + } + + inline double TransformedTriangle::calcStableT(const TetraCorner corner) const + { + // assert(_isTripleProductsCalculated); + // assert(_validTP[corner]); + return _tripleProducts[corner]; + } + + inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp) const + { + + // 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]; + + return _coords[5*pt1 + off1] * _coords[5*pt2 + off2] - _coords[5*pt1 + off2] * _coords[5*pt2 + off1]; + } + + // ---------------------------------------------------------------------------------- + // TransformedTriangle_intersect.cxx + // ---------------------------------------------------------------------------------- + inline bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const + { + return _triangleSurroundsEdgeCache[edge] && testEdgeIntersectsTriangle(edge); + } + + inline bool TransformedTriangle::testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const + { + return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet); + } + + inline bool TransformedTriangle::testSurfaceRayIntersection(const TetraCorner corner) const + { + return testTriangleSurroundsRay( corner ) && testSurfaceAboveCorner( corner ); + } + + inline bool TransformedTriangle::testCornerInTetrahedron(const TriCorner corner) const + { + const double pt[4] = + { + _coords[5*corner], // x + _coords[5*corner + 1], // y + _coords[5*corner + 2], // z + _coords[5*corner + 3] // z + }; + + for(int i = 0 ; i < 4 ; ++i) + { + if(pt[i] < 0.0 || pt[i] > 1.0) + { + return false; + } + } + return true; + } + + inline bool TransformedTriangle::testCornerOnXYZFacet(const TriCorner corner) const + { +#if 0 + const double pt[4] = + { + _coords[5*corner], // x + _coords[5*corner + 1], // y + _coords[5*corner + 2], // z + _coords[5*corner + 3] // h + }; +#endif + const double* pt = &_coords[5*corner]; + + if(pt[3] != 0.0) + { + return false; + } + + for(int i = 0 ; i < 3 ; ++i) + { + if(pt[i] < 0.0 || pt[i] > 1.0) + { + return false; + } + } + return true; + } + + inline bool TransformedTriangle::testCornerAboveXYZFacet(const TriCorner corner) const + { + const double x = _coords[5*corner]; + const double y = _coords[5*corner + 1]; + const double h = _coords[5*corner + 3]; + const double H = _coords[5*corner + 4]; + + return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0; + + } + + inline bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const + { + + // assert(edge < H01); + + // correspondence edge - triple products + // for edges OX, ..., ZX (Grandy, table III) + static const TetraCorner TRIPLE_PRODUCTS[12] = + { + X, O, // OX + Y, O, // OY + Z, O, // OZ + X, Y, // XY + Y, Z, // YZ + Z, X, // ZX + }; + + // Grandy, [16] + 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 + } + + inline bool TransformedTriangle::testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const + { +#if 0 + const double signs[3] = + { + SIGN_FOR_SEG_FACET_INTERSECTION[3*facet], + SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 1], + SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 2] + }; +#endif + + const double* signs = &SIGN_FOR_SEG_FACET_INTERSECTION[3*facet]; + const double c1 = signs[0]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet]); + const double c2 = signs[1]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]); + const double c3 = signs[2]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]); + + return (c1*c3 > 0.0) && (c2*c3 > 0.0); + } + + inline bool TransformedTriangle::testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const + { + // use correspondence facet a = 0 <=> offset for coordinate a in _coords + // and also correspondence segment AB => corner A + const double coord1 = _coords[5*seg + facet]; + const double coord2 = _coords[5*( (seg + 1) % 3) + facet]; + + //? should we use epsilon-equality here in second test? + LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 ); + + return (coord1*coord2 <= 0.0) && (coord1 != coord2); + } + + inline bool TransformedTriangle::testSegmentIntersectsHPlane(const TriSegment seg) const + { + // get the H - coordinates + const double coord1 = _coords[5*seg + 4]; + const double coord2 = _coords[5*( (seg + 1) % 3) + 4]; + //? should we use epsilon-equality here in second test? + LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 ); + + return (coord1*coord2 <= 0.0) && (coord1 != coord2); + } + + 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; + } + } + + inline bool TransformedTriangle::testTriangleSurroundsRay(const TetraCorner corner) const + { + // assert(corner == X || corner == Y || corner == Z); + + // double products to use for the possible corners + static const DoubleProduct DP_FOR_RAY_INTERSECTION[4] = + { + DoubleProduct(0), // O - only here to fill out and make indices match + C_10, // X + C_01, // Y + C_XY // Z + }; + + const DoubleProduct dp = DP_FOR_RAY_INTERSECTION[corner]; + + const double cPQ = calcStableC(PQ, dp); + const double cQR = calcStableC(QR, dp); + const double cRP = calcStableC(RP, dp); + + //? NB here we have no correction for precision - is this good? + // Our authority Grandy says nothing + LOG(5, "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" ); + + return ( cPQ*cQR > 0.0 ) && ( cPQ*cRP > 0.0 ); + + } -#include "TransformedTriangleInline.hxx" } diff --git a/src/INTERP_KERNEL/TransformedTriangleInline.hxx b/src/INTERP_KERNEL/TransformedTriangleInline.hxx deleted file mode 100644 index 28be48fa9..000000000 --- a/src/INTERP_KERNEL/TransformedTriangleInline.hxx +++ /dev/null @@ -1,284 +0,0 @@ -// Copyright (C) 2007-2024 CEA, EDF -// -// This library is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com -// - -#ifndef __TRANSFORMEDTRIANGLEINLINE_HXX__ -#define __TRANSFORMEDTRIANGLEINLINE_HXX__ - -// This file contains inline versions of some of the methods in the TransformedTriangle*.cxx files. -// It replaces those methods if OPTIMIZE is defined. -// NB : most of these methods have documentation in their corresponding .cxx - file. - -// ---------------------------------------------------------------------------------- -// Optimization methods. These are only defined and used if OPTIMIZE is defined. -// ----------------------------------------------------------------------------------- - - -inline void TransformedTriangle::preCalculateTriangleSurroundsEdge() -{ - for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1)) - { - _triangleSurroundsEdgeCache[edge] = testTriangleSurroundsEdge(edge); - } -} - - -// ---------------------------------------------------------------------------------- -// TransformedTriangle_math.cxx -// ---------------------------------------------------------------------------------- - -inline void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner) -{ - // set the three corresponding double products to 0.0 - static const DoubleProduct DOUBLE_PRODUCTS[12] = - { - C_YZ, C_ZX, C_XY, // O - C_YZ, C_ZH, C_YH, // X - C_ZX, C_ZH, C_XH, // Y - C_XY, C_YH, C_XH // Z - }; - - for(int i = 0 ; i < 3 ; ++i) { - const DoubleProduct dp = DOUBLE_PRODUCTS[3*corner + i]; - - LOG(6, std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner); - _doubleProducts[8*seg + dp] = 0.0; - }; -} - -inline double TransformedTriangle::calcStableC(const TriSegment seg, const DoubleProduct dp) const -{ - return _doubleProducts[8*seg + dp]; -} - -inline double TransformedTriangle::calcStableT(const TetraCorner corner) const -{ - // assert(_isTripleProductsCalculated); - // assert(_validTP[corner]); - return _tripleProducts[corner]; -} - -inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp) const -{ - - // 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]; - - return _coords[5*pt1 + off1] * _coords[5*pt2 + off2] - _coords[5*pt1 + off2] * _coords[5*pt2 + off1]; -} - -// ---------------------------------------------------------------------------------- -// TransformedTriangle_intersect.cxx -// ---------------------------------------------------------------------------------- -inline bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const -{ - return _triangleSurroundsEdgeCache[edge] && testEdgeIntersectsTriangle(edge); -} - -inline bool TransformedTriangle::testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const -{ - return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet); -} - -inline bool TransformedTriangle::testSurfaceRayIntersection(const TetraCorner corner) const -{ - return testTriangleSurroundsRay( corner ) && testSurfaceAboveCorner( corner ); -} - -inline bool TransformedTriangle::testCornerInTetrahedron(const TriCorner corner) const -{ - const double pt[4] = - { - _coords[5*corner], // x - _coords[5*corner + 1], // y - _coords[5*corner + 2], // z - _coords[5*corner + 3] // z - }; - - for(int i = 0 ; i < 4 ; ++i) - { - if(pt[i] < 0.0 || pt[i] > 1.0) - { - return false; - } - } - return true; -} - -inline bool TransformedTriangle::testCornerOnXYZFacet(const TriCorner corner) const -{ -#if 0 - const double pt[4] = - { - _coords[5*corner], // x - _coords[5*corner + 1], // y - _coords[5*corner + 2], // z - _coords[5*corner + 3] // h - }; -#endif - const double* pt = &_coords[5*corner]; - - if(pt[3] != 0.0) - { - return false; - } - - for(int i = 0 ; i < 3 ; ++i) - { - if(pt[i] < 0.0 || pt[i] > 1.0) - { - return false; - } - } - return true; -} - -inline bool TransformedTriangle::testCornerAboveXYZFacet(const TriCorner corner) const -{ - const double x = _coords[5*corner]; - const double y = _coords[5*corner + 1]; - const double h = _coords[5*corner + 3]; - const double H = _coords[5*corner + 4]; - - return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0; - -} - -inline bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const -{ - - // assert(edge < H01); - - // correspondence edge - triple products - // for edges OX, ..., ZX (Grandy, table III) - static const TetraCorner TRIPLE_PRODUCTS[12] = - { - X, O, // OX - Y, O, // OY - Z, O, // OZ - X, Y, // XY - Y, Z, // YZ - Z, X, // ZX - }; - - // Grandy, [16] - 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 -} - -inline bool TransformedTriangle::testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const -{ -#if 0 - const double signs[3] = - { - SIGN_FOR_SEG_FACET_INTERSECTION[3*facet], - SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 1], - SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 2] - }; -#endif - - const double* signs = &SIGN_FOR_SEG_FACET_INTERSECTION[3*facet]; - const double c1 = signs[0]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet]); - const double c2 = signs[1]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]); - const double c3 = signs[2]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]); - - return (c1*c3 > 0.0) && (c2*c3 > 0.0); -} - -inline bool TransformedTriangle::testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const -{ - // use correspondence facet a = 0 <=> offset for coordinate a in _coords - // and also correspondence segment AB => corner A - const double coord1 = _coords[5*seg + facet]; - const double coord2 = _coords[5*( (seg + 1) % 3) + facet]; - - //? should we use epsilon-equality here in second test? - LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 ); - - return (coord1*coord2 <= 0.0) && (coord1 != coord2); -} - -inline bool TransformedTriangle::testSegmentIntersectsHPlane(const TriSegment seg) const -{ - // get the H - coordinates - const double coord1 = _coords[5*seg + 4]; - const double coord2 = _coords[5*( (seg + 1) % 3) + 4]; - //? should we use epsilon-equality here in second test? - LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 ); - - return (coord1*coord2 <= 0.0) && (coord1 != coord2); -} - -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; - } -} - -inline bool TransformedTriangle::testTriangleSurroundsRay(const TetraCorner corner) const -{ - // assert(corner == X || corner == Y || corner == Z); - - // double products to use for the possible corners - static const DoubleProduct DP_FOR_RAY_INTERSECTION[4] = - { - DoubleProduct(0), // O - only here to fill out and make indices match - C_10, // X - C_01, // Y - C_XY // Z - }; - - const DoubleProduct dp = DP_FOR_RAY_INTERSECTION[corner]; - - const double cPQ = calcStableC(PQ, dp); - const double cQR = calcStableC(QR, dp); - const double cRP = calcStableC(RP, dp); - - //? NB here we have no correction for precision - is this good? - // Our authority Grandy says nothing - LOG(5, "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" ); - - return ( cPQ*cQR > 0.0 ) && ( cPQ*cRP > 0.0 ); - -} -#endif diff --git a/src/INTERP_KERNEL/TransformedTriangleIntersect.cxx b/src/INTERP_KERNEL/TransformedTriangleIntersect.cxx index d72592592..9fc3a1674 100644 --- a/src/INTERP_KERNEL/TransformedTriangleIntersect.cxx +++ b/src/INTERP_KERNEL/TransformedTriangleIntersect.cxx @@ -168,13 +168,12 @@ namespace INTERP_KERNEL const double alpha = tA / (tA - tB); // calculate point - LOG(4, "corner A = " << corners[0] << " corner B = " << corners[1] ); + LOG(4, "corner A = " << strTC(corners[0]) << " corner B = " << strTC(corners[1]) ); LOG(4, "tA = " << tA << " tB = " << tB << " alpha= " << alpha ); for(int i = 0; i < 3; ++i) { - pt[i] = (1 - alpha) * COORDS_TET_CORNER[3*corners[0] + i] + - alpha * COORDS_TET_CORNER[3*corners[1] + i]; + pt[i] = (1 - alpha)*COORDS_TET_CORNER[3*corners[0] + i] + alpha*COORDS_TET_CORNER[3*corners[1] + i]; #if 0 pt[i] = (1 - alpha) * getCoordinateForTetCorner() + alpha * getCoordinateForTetCorner(); @@ -234,7 +233,7 @@ namespace INTERP_KERNEL /** * Tests if the given segment of the triangle intersects the given edge of the tetrahedron (Grandy, eq. [20] - * If the OPTIMIZE is defined, it does not do the test the double product that should be zero. + * * @param seg segment of the triangle * @param edge edge of tetrahedron * @return true if the segment intersects the edge @@ -249,13 +248,13 @@ namespace INTERP_KERNEL for(int i = 0 ; i < 2 ; ++i) { facet[i] = FACET_FOR_EDGE[2*edge + i]; - + // find the two c-values -> the two for the other edges of the facet int idx1 = 0 ; int idx2 = 1; DoubleProduct dp1 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1]; DoubleProduct dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2]; - + if(dp1 == DoubleProduct( edge )) { idx1 = 2; @@ -266,7 +265,7 @@ namespace INTERP_KERNEL idx2 = 2; dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2]; } - + const double c1 = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1]*calcStableC(seg, dp1); const double c2 = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2]*calcStableC(seg, dp2); @@ -347,24 +346,15 @@ namespace INTERP_KERNEL const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx]; c[j] = dpIdx < 0.0 ? 0.0 : sign * calcStableC(seg, dp); } - - // pt[i] = (c1*s1 + c2*s2) / (s1^2 + s2^2) pt[i] = (c[0] * s[0] + c[1] * s[1]) / denominator; - - // strange bug with -O2 enabled : assertion fails when we don't have the following - // trace - line - //std::cout << "pt[i] = " << pt[i] << std::endl; - //assert(pt[i] >= 0.0); // check we are in tetraeder - //assert(pt[i] <= 1.0); - } } /** * Tests if the given segment of the triangle intersects the given corner of the tetrahedron. - * (Grandy, eq. [21]). If OPTIMIZE is defined, the double products that should be zero are not verified. + * (Grandy, eq. [21]). * * @param seg segment of the triangle * @param corner corner of the tetrahedron @@ -372,8 +362,6 @@ namespace INTERP_KERNEL */ bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const { - - // facets meeting at a given corner static const TetraFacet FACETS_FOR_CORNER[12] = { @@ -388,9 +376,7 @@ namespace INTERP_KERNEL { const TetraFacet facet = FACETS_FOR_CORNER[3*corner + i]; if(testSegmentIntersectsFacet(seg, facet)) - { - return true; - } + return true; } return false; @@ -432,11 +418,10 @@ namespace INTERP_KERNEL }; const TetraFacet facet = FACET_FOR_HALFSTRIP_INTERSECTION[edgeIndex]; - // special case : facet H = 0 const bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet); - LOG(4, "Halfstrip tests (" << seg << ", " << edge << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) ); + LOG(4, "Halfstrip tests (" << strTriS(seg) << ", " << strTE(edge) << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) ); LOG(4, "c2 = " << cVals[2] << ", c3 = " << cVals[3] ); return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0); @@ -497,7 +482,6 @@ namespace INTERP_KERNEL /** * Tests if the given segment of triangle PQR intersects the ray pointing * in the upwards z - direction from the given corner of the tetrahedron. (Grandy eq. [29]) - * If OPTIMIZE is defined, the double product that should be zero is not verified. * * @param seg segment of the triangle PQR * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z) @@ -573,6 +557,7 @@ namespace INTERP_KERNEL // if two or more c-values are zero we disallow x-edge intersection // Grandy, p.446 + // test for == 0.0 is OK here, if you look at the way double products were pre-comuted: const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0); if(numZeros >= 2 ) diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx index efb5937ba..f3ce8f3ad 100644 --- a/src/INTERP_KERNEL/VectorUtils.hxx +++ b/src/INTERP_KERNEL/VectorUtils.hxx @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -78,7 +79,7 @@ namespace INTERP_KERNEL inline const std::string vToStr(const double* pt) { std::stringstream ss(std::ios::out); - ss << "[" << pt[0] << ", " << pt[1] << ", " << pt[2] << "]"; + ss << std::setprecision(16) << "[" << pt[0] << ", " << pt[1] << ", " << pt[2] << "]"; return ss.str(); } diff --git a/src/INTERP_KERNELTest/MeshTestToolkit.txx b/src/INTERP_KERNELTest/MeshTestToolkit.txx index 063885089..49097166a 100644 --- a/src/INTERP_KERNELTest/MeshTestToolkit.txx +++ b/src/INTERP_KERNELTest/MeshTestToolkit.txx @@ -166,7 +166,7 @@ namespace INTERP_TEST LOG(1, "Source volume inconsistent : vol of cell " << i << " = " << sVol[i] << " but the row sum is " << sum_row ); ok = false; } - LOG(1, "diff = " < eff = countNumberOfMatrixEntries(m); // LOG(1, eff.first << " of " << numTargetElems * numSrcElems << " intersections calculated : ratio = " // << double(eff.first) / double(numTargetElems * numSrcElems)); - LOG(1, eff.second << " non-zero elements of " << eff.first << " total : filter efficiency = " - << double(eff.second) / double(eff.first)); +// LOG(1, eff.second << " non-zero elements of " << eff.first << " total : filter efficiency = " +// << double(eff.second) / double(eff.first)); LOG(1, "Intersection calculation done. " << std::endl ); diff --git a/src/INTERP_KERNELTest/TestInterpKernel.cxx b/src/INTERP_KERNELTest/TestInterpKernel.cxx index 6f3462721..1a885ad51 100644 --- a/src/INTERP_KERNELTest/TestInterpKernel.cxx +++ b/src/INTERP_KERNELTest/TestInterpKernel.cxx @@ -52,10 +52,10 @@ CPPUNIT_TEST_SUITE_REGISTRATION( UnitTetra3D2DIntersectionTest ); #ifndef MEDCOUPLING_MICROMED // These test suites need MEDLoader to load some test files: CPPUNIT_TEST_SUITE_REGISTRATION( InterpolationOptionsTest ); +CPPUNIT_TEST_SUITE_REGISTRATION( SingleElementTetraTests ); CPPUNIT_TEST_SUITE_REGISTRATION( HexaTests ); CPPUNIT_TEST_SUITE_REGISTRATION( MultiElement2DTests ); CPPUNIT_TEST_SUITE_REGISTRATION( MultiElementTetraTests ); -CPPUNIT_TEST_SUITE_REGISTRATION( SingleElementTetraTests ); CPPUNIT_TEST_SUITE_REGISTRATION( ThreeDSurfProjectionTest ); #endif // --- generic Main program from KERNEL_SRC/src/Basics/Test