From 03053f34c6fb4f471d70ce0a011970f3c61ccfc2 Mon Sep 17 00:00:00 2001 From: vbd Date: Tue, 7 Aug 2007 14:46:30 +0000 Subject: [PATCH] staffan : * corrections in Segment-Ray intersection and segment - halfstrip intersection --- src/INTERP_KERNEL/Interpolation3D.cxx | 2 +- src/INTERP_KERNEL/MeshElement.cxx | 54 +--------- src/INTERP_KERNEL/TransformedTriangle.cxx | 7 +- .../TransformedTriangle_intersect.cxx | 99 ++++++++++++------- .../TransformedTriangle_math.cxx | 39 +++++--- src/INTERP_KERNEL/VectorUtils.hxx | 1 + 6 files changed, 94 insertions(+), 108 deletions(-) diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx index 9355fbd27..59ae237a1 100644 --- a/src/INTERP_KERNEL/Interpolation3D.cxx +++ b/src/INTERP_KERNEL/Interpolation3D.cxx @@ -157,7 +157,7 @@ namespace MEDMEM iter != currNode->getTargetRegion().getEndElements() ; ++iter) { const double vol = calculateIntersectionVolume(*srcElement, **iter); - // if(vol != 0.0) + if(vol != 0.0) { const int targetIdx = (*iter)->getIndex(); diff --git a/src/INTERP_KERNEL/MeshElement.cxx b/src/INTERP_KERNEL/MeshElement.cxx index 1501b5322..f03ffbfbd 100644 --- a/src/INTERP_KERNEL/MeshElement.cxx +++ b/src/INTERP_KERNEL/MeshElement.cxx @@ -157,6 +157,7 @@ namespace INTERP_UTILS } + // to be removed assert(faceType == MED_TRIA3); // create transformed triangles from face @@ -165,6 +166,8 @@ namespace INTERP_UTILS case MED_TRIA3: triangles.push_back(TransformedTriangle(&transformedNodes[0], &transformedNodes[3], &transformedNodes[6])); break; + + // add other cases here to treat hexahedra, pyramides, etc default: std::cout << "Only elements with triangular faces are supported at the moment." << std::endl; @@ -174,57 +177,6 @@ namespace INTERP_UTILS } - -#if 0 - void MeshElement::triangulate(std::vector& triangles, const TetraAffineTransform& T) const - { - std::cout << "Triangulating element " << getIndex() << std::endl; - switch(_type) - { - case MED_TETRA4 : - // triangles == faces - const int beginFaceIdx = _mesh->getConnectivityIndex(MED_DESCENDING, MED_CELL)[_index]; - const int endFaceIdx = _mesh->getConnectivityIndex(MED_DESCENDING, MED_CELL)[_index + 1]; - - std::cout << "elements has faces at indices " << beginFaceIdx << " to " << endFaceIdx << std::endl; - assert(endFaceIdx - beginFaceIdx == 4); - - for(int i = beginFaceIdx ; i < endFaceIdx ; ++i) // loop over faces of element - { - const int faceIdx = _mesh->getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL, MED_ALL_ELEMENTS)[i - 1]; - const int beginNodeIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_FACE)[faceIdx - 1]; - const int endNodeIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_FACE)[faceIdx]; - std::cout << "Face " << faceIdx << " with nodes in [" << beginNodeIdx << "," << endNodeIdx << "[" << std::endl; - assert(endNodeIdx - beginNodeIdx == 3); - - double transformedPts[9]; - - for(int j = 0 ; j < 3 ; ++j) // loop over nodes of face - { - // { optimise here using getCoordinatesForNode ? - // could maybe use the connNodeIdx directly and only transform each node once - // instead of once per face - const int connNodeIdx = - _mesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL, MED_FACE, MED_ALL_ELEMENTS)[beginNodeIdx + j - 1] - 1; - const double* pt = &(_mesh->getCoordinates(MED_FULL_INTERLACE)[3*connNodeIdx]); - - //const double* pt = getCoordsOfNode(j + 1); - - // transform - T.apply(&transformedPts[3*j], pt); - std::cout << "Transforming : " << vToStr(pt) << " -> " << vToStr(&transformedPts[3*j]) << std::endl; - } - - triangles.push_back(TransformedTriangle(&transformedPts[0], &transformedPts[3], &transformedPts[6])); - } - - break; - default: - break; - } - } -#endif - int MeshElement::getIndex() const { return _index + 1; diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index d2807adc7..da0b401bb 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -183,14 +183,13 @@ namespace INTERP_UTILS std::cout << std::endl << "-- Calculating intersection polygons ... " << std::endl; calculateIntersectionPolygons(); - // calculate volume under A - std::cout << std::endl << "-- Treating polygon A ... " << std::endl; double barycenter[3]; - - double volA = 0.0; + // calculate volume under A + double volA = 0.0; if(_polygonA.size() > 2) { + std::cout << std::endl << "-- Treating polygon A ... " << std::endl; calculatePolygonBarycenter(A, barycenter); sortIntersectionPolygon(A, barycenter); volA = calculateVolumeUnderPolygon(A, barycenter); diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx index b601b7fe2..a0d8adef2 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx @@ -3,6 +3,7 @@ #include #include #include +#include "VectorUtils.hxx" namespace INTERP_UTILS { @@ -168,8 +169,9 @@ namespace INTERP_UTILS const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx]; const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx]; pt[i] = -( sign * calcStableC(seg, dp) ) / s; - assert(pt[i] > 0.0); - assert(pt[i] < 1.0); + // std::cout << "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i] << std::endl; + assert(pt[i] >= 0.0); + assert(pt[i] <= 1.0); } } @@ -394,12 +396,18 @@ namespace INTERP_UTILS // products 3 and 4 for each edge -> third condition // NB : some uncertainty whether these last are correct static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12] = + { + C_10, C_01, C_ZH, C_10, // XY + C_01, C_XY, C_XH, C_01, // YZ + C_XY, C_10, C_YH, C_XY // ZX + }; + /*static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12] = { C_10, C_01, C_ZH, C_10, // XY C_01, C_XY, C_XH, C_ZX, // YZ C_XY, C_10, C_YH, C_XY // ZX }; - + */ /*static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12] = { C_10, C_01, C_ZH, C_XY, // XY @@ -491,7 +499,7 @@ namespace INTERP_UTILS assert(pt[i] >= 0.0); assert(pt[i] <= 1.0); } - assert(std::abs(pt[0] + pt[1] + pt[2] - 1.0) < 1.0e-9); + assert(epsilonEqual(pt[0] + pt[1] + pt[2] - 1.0, 0.0, 1.0e-9)); } /** @@ -505,6 +513,7 @@ namespace INTERP_UTILS bool TransformedTriangle::testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const { assert(corner == X || corner == Y || corner == Z); + // std::cout << "Testing seg - ray intersection for seg = " << seg << ", corner = " << corner << std::endl; // readjust index since O is not used const int cornerIdx = static_cast(corner) - 1; @@ -513,56 +522,74 @@ namespace INTERP_UTILS // dp 1 -> cond 1 // dp 2-7 -> cond 3 //? NB : last two rows are not completely understood and may contain errors + /*static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21] = + { + C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X + C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_10, // Y + C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_01 // Z + };*/ static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21] = + { + C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X + C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_YZ, // Y + C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_ZX // Z + }; + /*static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21] = { C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X C_01, C_XH, C_ZH, C_XY, C_10, C_XH, C_ZX, // Y C_XY, C_YH, C_XH, C_10, C_01, C_ZH, C_YZ // Z }; - + */ // facets to use //? not sure this is correct - static const TetraFacet FACET_SEGMENT_RAY_INTERSECTION[3] = + static const TetraFacet FIRST_FACET_SEGMENT_RAY_INTERSECTION[3] = { OZX, // X - OXY, // Y - OYZ, // Z + OYZ, // Y + OZX, // Z }; const DoubleProduct dp0 = DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx]; + const double cVal0 = calcStableC(seg, dp0); //? epsilon-equality here? - if(dp0 == 0.0) // cond. 1 + // cond. 1 + if(cVal0 != 0.0) { + // std::cout << "SR fails at cond 1 cVal0 = " << cVal0 << std::endl; + return false; + } - if(testSegmentIntersectsFacet(seg, FACET_SEGMENT_RAY_INTERSECTION[cornerIdx])) // cond. 2.1 - { - const double H1 = _coords[5*seg + 4]; - const double H2 = _coords[5*( (seg + 1) % 3) + 4]; - - // S_H -> cond. 2.2 - // std::cout << "H1 = " << H1 << ", H2= " << H2 << std::endl; - if(H1*H2 <= 0.0 && H1 != H2) // should equality be in "epsilon-sense" ? - { + // cond 2 + const bool cond21 = testSegmentIntersectsFacet(seg, FIRST_FACET_SEGMENT_RAY_INTERSECTION[cornerIdx]); + const bool cond22 = (corner == Z) ? testSegmentIntersectsFacet(seg, OYZ) : testSegmentIntersectsHPlane(seg); - const double cVals[6] = - { - calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 1]), - calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 2]), - calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 3]), - calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 4]), - calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 5]), - calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 6]), - }; - - // cond. 3 - return ( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) < 0.0; - } - } + if(!(cond21 || cond22)) + { + // std::cout << "SR fails at cond 2 : cond21 = " << cond21 << ", cond22 = " << cond22 << std::endl; + return false; } - return false; + // cond 3 + const double cVals[6] = + { + calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 1]), + calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 2]), + calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 3]), + calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 4]), + calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 5]), + calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 6]), + }; + + // cond. 3 + if(( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) >= 0.0) + { + std::cout << "SR fails at cond 3 : " << (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] << std::endl; + } + return ( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) < 0.0; + } /** @@ -631,7 +658,7 @@ namespace INTERP_UTILS 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; + return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0; } @@ -775,8 +802,8 @@ namespace INTERP_UTILS // std::cout << "surface above corner " << corner << " : " << "n = " << normal << ", t = [" << calcTByDevelopingRow(corner, 1, false) << ", " << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) << "] - stable : " << calcStableT(corner) << std::endl; //? 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]) - // return ( calcTByDevelopingRow(corner, 2, false) * normal ) >= 0.0; - return ( calcStableT(corner) * normal ) >= 0.0; + return ( calcTByDevelopingRow(corner, 2, false) * normal ) >= 0.0; + //return ( calcStableT(corner) * normal ) >= 0.0; } /** diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index c513a2432..6a5b4ba46 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -4,6 +4,7 @@ #include #include #include +#include "VectorUtils.hxx" namespace INTERP_UTILS { @@ -129,16 +130,12 @@ namespace INTERP_UTILS //{ use functions in VectorUtils for this // cross product of difference vectors - const double cross[3] = - { - diffPQ[1]*diffCornerP[2] - diffPQ[2]*diffCornerP[1], - diffPQ[2]*diffCornerP[0] - diffPQ[0]*diffCornerP[2], - diffPQ[0]*diffCornerP[1] - diffPQ[1]*diffCornerP[0], - }; - + + double crossProd[3]; + cross(diffPQ, diffCornerP, crossProd); - const double cross_squared = cross[0]*cross[0] + cross[1]*cross[1] + cross[2]*cross[2]; - const double norm_diffPQ_squared = diffPQ[0]*diffPQ[0] + diffPQ[1]*diffPQ[1] + diffPQ[2]*diffPQ[2]; + const double cross_squared = dot(crossProd, crossProd); + const double norm_diffPQ_squared = dot(diffPQ, diffPQ); const double dist = cross_squared / norm_diffPQ_squared; // update minimum (initializs with first corner) @@ -351,7 +348,7 @@ namespace INTERP_UTILS double TransformedTriangle::calcStableT(const TetraCorner corner) const { assert(_isTripleProductsCalculated); - // assert(_validTP[corner]); + assert(_validTP[corner]); return _tripleProducts[corner]; } @@ -451,7 +448,7 @@ namespace INTERP_UTILS const double cPQ = calcStableC(PQ, dp); // coordinate to use for projection (Grandy, [57]) with edges - // OX, OY, OZ, YZ, ZX, XY in order : + // 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 } ; @@ -476,11 +473,21 @@ namespace INTERP_UTILS const double p_term = _coords[5*P + offset] * cQR * (1.0 - alpha * coordValues[0] * cQR) ; const double q_term = _coords[5*Q + offset] * cRP * (1.0 - alpha * coordValues[1] * cRP) ; const double r_term = _coords[5*R + offset] * cPQ * (1.0 - alpha * coordValues[2] * cPQ) ; - - // NB : using plus also for the middle term compensates for a double product - // which is inversely ordered - // std::cout << "Triple product for corner " << corner << ", row " << row << " = " << sign*( p_term + q_term + r_term ) << std::endl; - return sign*( p_term + q_term + r_term ); + + // check that we are not so close to zero that numerical errors could take over, + // otherwise reset to zero (cf Grandy, p. 446) + const long double delta = MULT_PREC_F * (std::abs(p_term) + std::abs(q_term) + std::abs(r_term)); + if( epsilonEqual( p_term + q_term + r_term, 0.0, THRESHOLD_F * delta) ) + { + return 0.0; + } + else + { + // NB : using plus also for the middle term compensates for a double product + // which is inversely ordered + // std::cout << "Triple product for corner " << corner << ", row " << row << " = " << sign*( p_term + q_term + r_term ) << std::endl; + return sign*( p_term + q_term + r_term ); + } } diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx index af5159671..f8e58d9d8 100644 --- a/src/INTERP_KERNEL/VectorUtils.hxx +++ b/src/INTERP_KERNEL/VectorUtils.hxx @@ -5,6 +5,7 @@ #include #include #include +#include namespace INTERP_UTILS { -- 2.39.2