From 0340465720812af544a1a60882b932c852279ac3 Mon Sep 17 00:00:00 2001 From: vbd Date: Wed, 22 Aug 2007 07:59:32 +0000 Subject: [PATCH] staffan : * added compile-time switches for different experimental pieces of code * fixed bug consistency test for double products which caused many tests to fail -> this seems to have made some of the experimental code unnecessary * removed log messages as these are encumbering when running tests with large meshes status : all pertinent tests pass ! --- src/INTERP_KERNEL/BoundingBox.cxx | 7 +- src/INTERP_KERNEL/Interpolation3D.cxx | 19 ++-- src/INTERP_KERNEL/MeshElement.cxx | 8 +- src/INTERP_KERNEL/TetraAffineTransform.hxx | 26 ++++- src/INTERP_KERNEL/TransformedTriangle.cxx | 98 ++++++++++++++++-- src/INTERP_KERNEL/TransformedTriangle.hxx | 5 + .../TransformedTriangle_intersect.cxx | 77 ++++++++++----- .../TransformedTriangle_math.cxx | 99 +++++++++++++------ 8 files changed, 261 insertions(+), 78 deletions(-) diff --git a/src/INTERP_KERNEL/BoundingBox.cxx b/src/INTERP_KERNEL/BoundingBox.cxx index 2a43dceb5..3817afb3a 100644 --- a/src/INTERP_KERNEL/BoundingBox.cxx +++ b/src/INTERP_KERNEL/BoundingBox.cxx @@ -84,15 +84,18 @@ namespace INTERP_UTILS // minimum coordinate of one is greater than the maximum coordinate of the other //? stable version - /* const double tol = 1.0e-2; + /*const double tol = 1.0e-2*_coords[c]; if(_coords[c] > otherMaxCoord + tol || _coords[c + 3] < otherMinCoord - tol) */ + if(_coords[c] > otherMaxCoord || _coords[c + 3] < otherMinCoord) - { + + { return true; } + } return false; } diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx index 84aeac44f..8d96f5d45 100644 --- a/src/INTERP_KERNEL/Interpolation3D.cxx +++ b/src/INTERP_KERNEL/Interpolation3D.cxx @@ -7,6 +7,7 @@ #include "RegionNode.hxx" #include "TetraAffineTransform.hxx" #include "TransformedTriangle.hxx" +#include "VectorUtils.hxx" using namespace INTERP_UTILS; using namespace MEDMEM; @@ -157,7 +158,8 @@ namespace MEDMEM iter != currNode->getTargetRegion().getEndElements() ; ++iter) { const double vol = calculateIntersectionVolume(*srcElement, **iter); - if(vol != 0.0) + if(!epsilonEqual(vol, 0.0, 1.0e-10)) + // if(vol != 0.0) { const int targetIdx = (*iter)->getIndex(); @@ -281,11 +283,11 @@ namespace MEDMEM // (a), (b) and (c) not yet implemented // (d) : without fine-level filtering (a) - (c) for the time being - - // std::cout << "Source : "; - // srcElement.dumpCoords(); - // std::cout << "Target : "; - // targetElement.dumpCoords(); + // std::cout << "------------------" << std::endl; + // std::cout << "Source : "; + srcElement.dumpCoords(); + // std::cout << "Target : "; + targetElement.dumpCoords(); // get array of points of target tetraeder const double* tetraCorners[4]; @@ -298,7 +300,7 @@ namespace MEDMEM TetraAffineTransform T( tetraCorners ); // check if we have planar tetra element - if(T.determinant() == 0.0) + if(epsilonEqual(T.determinant(), 0.0, 1.0e-16)) { // tetra is planar // std::cout << "Planar tetra -- volume 0" << std::endl; @@ -321,7 +323,7 @@ namespace MEDMEM for(vector::iterator iter = triangles.begin() ; iter != triangles.end(); ++iter) { // std::cout << std::endl << "= > Triangle " << ++i << std::endl; - // iter->dumpCoords(); + iter->dumpCoords(); volume += iter->calculateIntersectionVolume(); } @@ -332,6 +334,7 @@ namespace MEDMEM //? fault in article, Grandy, [8] : it is the determinant of the inverse transformation that // should be used + return std::abs(1.0 / T.determinant() * volume) ; } diff --git a/src/INTERP_KERNEL/MeshElement.cxx b/src/INTERP_KERNEL/MeshElement.cxx index 072c4ce75..32077a5c7 100644 --- a/src/INTERP_KERNEL/MeshElement.cxx +++ b/src/INTERP_KERNEL/MeshElement.cxx @@ -162,9 +162,11 @@ namespace INTERP_UTILS assert(faceType == MED_TRIA3); // create transformed triangles from face + switch(faceType) { case MED_TRIA3: + // std::cout << std::endl<< "** Adding triangle " << i << std::endl; triangles.push_back(TransformedTriangle(&transformedNodes[0], &transformedNodes[3], &transformedNodes[6])); break; @@ -186,12 +188,12 @@ namespace INTERP_UTILS void MeshElement::dumpCoords() const { - std::cout << "Element " << _index + 1 << " has nodes " << std::endl; + // std::cout << "Element " << _index + 1 << " has nodes " << std::endl; for(int i = 1 ; i <= getNumberNodes() ; ++i) { - std::cout << vToStr(getCoordsOfNode(i)) << ", "; + // std::cout << vToStr(getCoordsOfNode(i)) << ", "; } - std::cout << std::endl; + // std::cout << std::endl; } diff --git a/src/INTERP_KERNEL/TetraAffineTransform.hxx b/src/INTERP_KERNEL/TetraAffineTransform.hxx index c7263e10e..6b781bf63 100644 --- a/src/INTERP_KERNEL/TetraAffineTransform.hxx +++ b/src/INTERP_KERNEL/TetraAffineTransform.hxx @@ -6,6 +6,8 @@ #include "VectorUtils.hxx" +#undef NULL_COORD_CORRECTION + namespace INTERP_UTILS { @@ -69,7 +71,7 @@ namespace INTERP_UTILS // and O is the position vector of the point that is mapped onto the origin for(int i = 0 ; i < 3 ; ++i) { - _translation[i] = -_linearTransform[3*i]*(pts[0])[0] - _linearTransform[3*i+1]*(pts[0])[1] - _linearTransform[3*i+2]*(pts[0])[2] ; + _translation[i] = -(_linearTransform[3*i]*(pts[0])[0] + _linearTransform[3*i+1]*(pts[0])[1] + _linearTransform[3*i+2]*(pts[0])[2]) ; } // precalculate determinant (again after inversion of transform) @@ -115,6 +117,14 @@ namespace INTERP_UTILS // translation dest[i] += _translation[i]; + +#ifdef NULL_COORD_CORRECTION + if(epsilonEqual(dest[i], 0.0, 1.0e-15)) + { + dest[i] = 0.0; + } +#endif + } if(selfAllocation) @@ -154,6 +164,7 @@ namespace INTERP_UTILS { //{ we copy the matrix for the lu-factorization // maybe inefficient + double lu[9]; for(int i = 0 ; i < 9; ++i) { @@ -162,7 +173,11 @@ namespace INTERP_UTILS // calculate LU factorization int idx[3]; - factorizeLU(lu, idx); + // double parity = 1.0; + factorizeLU(lu, idx/*, &parity*/); + + //_determinant = parity / (lu[3*idx[0]] * lu[3*idx[1] + 1] * lu[3*idx[2] + 2]); + //std::cout << "lu-determinant = " << _determinant << std::endl; // calculate inverse by forward and backward substitution // store in _linearTransform @@ -212,7 +227,7 @@ namespace INTERP_UTILS ///////////////////////////////////////////////// /// Auxiliary methods for inverse calculation /// ///////////////////////////////////////////////// - void factorizeLU(double* lu, int* idx) const + void factorizeLU(double* lu, int* idx/*, double* parity*/) const { // 3 x 3 LU factorization // initialise idx @@ -242,6 +257,9 @@ namespace INTERP_UTILS int tmp = idx[k]; idx[k] = idx[row]; idx[row] = tmp; + /* if(k != row) + *parity = -(*parity); + */ // calculate row for(int j = k + 1 ; j < 3 ; ++j) @@ -250,7 +268,7 @@ namespace INTERP_UTILS lu[3*idx[j] + k] /= lu[3*idx[k] + k]; for(int s = k + 1; s < 3 ; ++s) { - // case s = k will always become zero, and then replaced by + // case s = k will always become zero, and then be replaced by // the l-value // there's no need to calculate this explicitly diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index 0e73c592d..6f2a9650e 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -4,8 +4,13 @@ #include #include #include +#include #include "VectorUtils.hxx" #include +#include + +#undef MERGE_CALC +#undef COORDINATE_CORRECTION 1.0e-15 class CircularSortOrder { @@ -84,6 +89,16 @@ private: const double _a, _b; }; +class Vector3Cmp +{ + public: + bool operator()(double* const& pt1, double* const& pt2) + { + // std::cout << "points are equal ? : " << int((pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2])) << std::endl; + return (pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2]); + } +}; + namespace INTERP_UTILS { //////////////////////////////////////////////////////////////////////////// @@ -113,15 +128,34 @@ namespace INTERP_UTILS } // h coordinate + _coords[5*P + 3] = 1 - p[0] - p[1] - p[2]; _coords[5*Q + 3] = 1 - q[0] - q[1] - q[2]; _coords[5*R + 3] = 1 - r[0] - r[1] - r[2]; - + + // order of substractions give different results ? + /* + _coords[5*P + 3] = 1 - (p[0] - (p[1] - p[2])); + _coords[5*Q + 3] = 1 - (q[0] - (q[1] - q[2])); + std::cout << "old = " << 1 - q[0] - q[1] - q[2] << " calculated = " << 1 - (q[0] - (q[1] - q[2])) << " stored : " << _coords[5*Q + 3] << std::endl; + _coords[5*R + 3] = 1 - (r[0] -(r[1] - r[2])); + */ + + // 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]; +#ifdef COORDINATE_CORRECTION + + // correction of small (imprecise) coordinate values + for(int i = 0 ; i < 15 ; ++i) + { + _coords[i] = epsilonEqual(_coords[i], 0.0, COORDINATE_CORRECTION) ? 0.0 : _coords[i]; + } +#endif + // initialise rest of data preCalculateDoubleProducts(); preCalculateTripleProducts(); @@ -170,12 +204,14 @@ namespace INTERP_UTILS double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2]; + if(sign == 0.0) { // std::cout << std::endl << "Triangle is perpendicular to z-plane - V = 0.0" << std::endl << std::endl; return 0.0; } + // normalize sign = sign > 0.0 ? 1.0 : -1.0; @@ -183,6 +219,17 @@ namespace INTERP_UTILS // std::cout << std::endl << "-- Calculating intersection polygons ... " << std::endl; calculateIntersectionPolygons(); +#ifdef MERGE_CALC + const bool mergeCalculation = isPolygonAOnHFacet(); + if(mergeCalculation) + { + // move points in B to A to avoid missing points + // NB : need to remove elements from B in order to handle deletion properly + _polygonA.insert(_polygonA.end(), _polygonB.begin(), _polygonB.end()); + _polygonB.clear(); + } +#endif + double barycenter[3]; // calculate volume under A @@ -197,11 +244,16 @@ namespace INTERP_UTILS } double volB = 0.0; - // if triangle is not in h = 0 plane, calculate volume under B +#ifdef MERGE_CALC + if((!mergeCalculation) && _polygonB.size() > 2) +#else if(!isTriangleInPlaneOfFacet(XYZ) && _polygonB.size() > 2) +#endif { // std::cout << std::endl << "-- Treating polygon B ... " << std::endl; + // std::cout << _coords[5*P + 3] << ", " << _coords[5*Q + 3] << ", " << _coords[5*R+ 3] << std::endl; + calculatePolygonBarycenter(B, barycenter); sortIntersectionPolygon(B, barycenter); volB = calculateVolumeUnderPolygon(B, barycenter); @@ -519,7 +571,7 @@ namespace INTERP_UTILS vol += (factor1 * factor2) / 6.0; } - // // std::cout << "Abs. Volume is " << vol << std::endl; + // std::cout << "Abs. Volume is " << vol << std::endl; return vol; } @@ -543,6 +595,7 @@ namespace INTERP_UTILS for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1)) { // ? should have epsilon-equality here? + //if(!epsilonEqual(_coords[5*c + coord], 0.0, 1.0e-15)) if(_coords[5*c + coord] != 0.0) { return false; @@ -552,6 +605,39 @@ namespace INTERP_UTILS return true; } + bool TransformedTriangle::isPolygonAOnHFacet() const + { + // need to have vector of unique points in order to determine the "real" number of + // of points in the polygon, to avoid problems when polygon A has less than 3 points + + using ::Vector3Cmp; + std::vector pAUnique; + pAUnique.reserve(_polygonA.size()); + Vector3Cmp cmp; + unique_copy(_polygonA.begin(), _polygonA.end(), back_inserter(pAUnique), cmp); + //for(std::vector::const_iterator iter = pAUnique.begin() ; iter != pAUnique.end() ; ++iter) + //std::cout << "next : " << vToStr(*iter) << std::endl; + + // std::cout << "paunique size = " << pAUnique.size() << std::endl; + if(pAUnique.size() < 3) + { + return false; + } + for(std::vector::const_iterator iter = _polygonA.begin() ; iter != _polygonA.end() ; ++iter) + { + const double* pt = *iter; + const double h = 1.0 - pt[0] - pt[1] - pt[2]; + //// std::cout << "h = " << h << std::endl; + + //if(h != 0.0) + if(!epsilonEqual(h, 0.0)) + { + return false; + } + } + // std::cout << "Polygon A is on h = 0 facet" << std::endl; + return true; + } bool TransformedTriangle::isTriangleBelowTetraeder() { @@ -568,12 +654,12 @@ namespace INTERP_UTILS void TransformedTriangle::dumpCoords() { - std::cout << "Coords : "; + // 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; + // std::cout << std::endl; } }; // NAMESPACE diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 771da3ff4..6797310eb 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -89,6 +89,9 @@ namespace INTERP_UTILS bool isTriangleBelowTetraeder(); + bool isPolygonAOnHFacet() const; + + //////////////////////////////////////////////////////////////////////////////////// /// Intersection test methods and intersection point calculations //////// //////////////////////////////////////////////////////////////////////////////////// @@ -145,6 +148,8 @@ namespace INTERP_UTILS void preCalculateDoubleProducts(void); + void resetDoubleProducts(const TriSegment seg, const TetraCorner corner); + double calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const; void preCalculateTripleProducts(void); diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx index 9b4ef0e72..300ee72c0 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx @@ -5,6 +5,9 @@ #include #include "VectorUtils.hxx" +#define SEG_RAY_TABLE 1 // seems correct +#undef TEST_EPS 0.0//1.0e-14 + namespace INTERP_UTILS { @@ -94,8 +97,6 @@ namespace INTERP_UTILS // : (x,y,z)* = (1-alpha)*A + alpha*B where // alpha = t_A / (t_A - t_B) - - const TetraCorner corners[2] = { CORNERS_FOR_EDGE[2*edge], @@ -108,12 +109,13 @@ namespace INTERP_UTILS const double alpha = tA / (tA - tB); // calculate point + // std::cout << "corner A = " << corners[0] << " corner B = " << corners[1] << std::endl; + // std::cout << "tA = " << tA << " tB = " << tB << " alpha= " << alpha << std::endl; for(int i = 0; i < 3; ++i) { - // std::cout << "tA = " << tA << " tB = " << tB << " alpha= " << alpha << std::endl; pt[i] = (1 - alpha) * COORDS_TET_CORNER[3*corners[0] + i] + alpha * COORDS_TET_CORNER[3*corners[1] + i]; - // std::cout << pt[i] << std::endl; + // // std::cout << pt[i] << std::endl; assert(pt[i] >= 0.0); assert(pt[i] <= 1.0); } @@ -169,7 +171,8 @@ 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; - // std::cout << "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i] << std::endl; + // std::cout << "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i] << std::endl; + // std::cout << "c(" << seg << ", " << dp << ") = " << sign * calcStableC(seg, dp) << std::endl; assert(pt[i] >= 0.0); assert(pt[i] <= 1.0); } @@ -197,7 +200,8 @@ namespace INTERP_UTILS OZX, XYZ // ZX }; - if(calcStableC(seg,DoubleProduct( edge )) != 0.0) + if(calcStableC(seg,DoubleProduct( edge )) != 0.0) + // if(!epsilonEqual(calcStableC(seg,DoubleProduct( edge )), 0.0, TEST_EPS)) { return false; } @@ -216,6 +220,7 @@ namespace INTERP_UTILS 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; @@ -227,11 +232,16 @@ namespace INTERP_UTILS 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); + 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); - isFacetCondVerified = isFacetCondVerified || c1*c2 > 0.0; + //isFacetCondVerified = isFacetCondVerified || c1*c2 > 0.0; + if(c1*c2 > 0.0) + { + isFacetCondVerified = true; + } } + if(!isFacetCondVerified) { return false; @@ -347,6 +357,7 @@ namespace INTERP_UTILS const DoubleProduct dp = DoubleProduct( edge ); const double c = calcStableC(seg, dp); if(c != 0.0) + // if(!epsilonEqual(c, 0.0, TEST_EPS)) { return false; } @@ -522,25 +533,23 @@ 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] = +#if SEG_RAY_TABLE==1 + 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 - };*/ + }; +#else + 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 - }; - */ +#endif + // facets to use //? not sure this is correct static const TetraFacet FIRST_FACET_SEGMENT_RAY_INTERSECTION[3] = @@ -688,7 +697,8 @@ namespace INTERP_UTILS // if two or more c-values are zero we disallow x-edge intersection // Grandy, p.446 - const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0); + const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0); + //? const int numZeros = (epsilonEqual(cPQ,0.0) ? 1 : 0) + (epsilonEqual(cQR,0.0) ? 1 : 0) + (epsilonEqual(cRP, 0.0) ? 1 : 0); if(numZeros >= 2 ) { @@ -717,7 +727,7 @@ namespace INTERP_UTILS X, O, // OX Y, O, // OY Z, O, // OZ - X, Y, // XY + X, Y, // XY Y, Z, // YZ Z, X, // ZX }; @@ -753,7 +763,16 @@ namespace INTERP_UTILS 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); + if(false) + //if(epsilonEqual(c1, 0.0, TEST_EPS) || epsilonEqual(c2, 0.0, TEST_EPS) || epsilonEqual(c3, 0.0, TEST_EPS)) + { + return false; + } + else + { + return (c1*c3 > 0.0) && (c2*c3 > 0.0); + } + } /** @@ -773,6 +792,7 @@ namespace INTERP_UTILS //? should we use epsilon-equality here in second test? // std::cout << "coord1 : " << coord1 << " coord2 : " << coord2 << std::endl; + //return (coord1*coord2 <= 0.0) && epsilonEqual(coord1,coord2); return (coord1*coord2 <= 0.0) && (coord1 != coord2); } @@ -799,11 +819,20 @@ namespace INTERP_UTILS //? is it always YZ here ? //? changed to XY ! const double normal = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY); - // 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; + //std::cout << "surface above corner " << corner << " : " << "n = " << normal << ", t = [" << calcTByDevelopingRow(corner, 1, false) << ", " << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) << std::endl; + //std::cout << "] - 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, 1, false) * normal ) >= 0.0; - //return ( calcStableT(corner) * normal ) >= 0.0; + //return ( calcTByDevelopingRow(corner, 1, false) * normal ) >= 0.0; + if(!_validTP[corner]) + { + return ( calcTByDevelopingRow(corner, 1, false) * normal ) >= 0.0; + } + else + { + return ( calcStableT(corner) * normal ) >= 0.0; + } } /** diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index aab69e295..e003671f6 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -9,6 +9,10 @@ #include "VectorUtils.hxx" +#undef SECOND_CORNER_RESET +#undef FIXED_DELTA 1.0e-15 + + namespace INTERP_UTILS { @@ -77,8 +81,8 @@ namespace INTERP_UTILS const double term3 = _doubleProducts[8*seg + C_XY] * _doubleProducts[8*seg + C_ZH]; // std::cout << std::endl; - //std::cout << "for seg " << seg << " consistency " << term1 + term2 + term3 << std::endl; - //std::cout << "term1 :" << term1 << " term2 :" << term2 << " term3: " << term3 << std::endl; + // std::cout << "for seg " << seg << " consistency " << term1 + term2 + term3 << std::endl; + // std::cout << "term1 :" << term1 << " term2 :" << term2 << " term3: " << term3 << std::endl; // if(term1 + term2 + term3 != 0.0) const int num_zero = (term1 == 0.0 ? 1 : 0) + (term2 == 0.0 ? 1 : 0) + (term3 == 0.0 ? 1 : 0); @@ -89,7 +93,7 @@ namespace INTERP_UTILS // * two terms zero // * all terms positive // * all terms negative - if((num_zero == 1 && num_neg != 1) || num_zero == 2 || num_neg == 0 || num_neg == 3 ) + if((num_zero == 1 && num_neg != 1) || num_zero == 2 || (num_neg == 0 && num_zero != 3) || num_neg == 3 ) { std::map distances; //std::cout << "inconsistent! "<< std::endl; @@ -104,23 +108,18 @@ namespace INTERP_UTILS // first element -> minimum distance const TetraCorner minCorner = distances.begin()->second; - // set the three corresponding double products to 0.0 - static const DoubleProduct DOUBLE_PRODUCTS[12] = + resetDoubleProducts(seg, minCorner); + // test : reset also for second corner if distance is small enough +#ifdef SECOND_CORNER_RESET + std::map::const_iterator iter = distances.begin(); + ++iter; + if(iter->first < 1.0e-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) { - DoubleProduct dp = DOUBLE_PRODUCTS[3*minCorner + i]; - - // std::cout << std::endl << "inconsistent dp :" << dp << std::endl; - _doubleProducts[8*seg + dp] = 0.0; - } - + resetDoubleProducts(seg, iter->second); + } +#endif } + } @@ -128,7 +127,7 @@ namespace INTERP_UTILS for(TriSegment seg = PQ ; seg <= RP ; seg = TriSegment(seg + 1)) { - for(DoubleProduct dp = C_YZ ; dp <= C_10 ; dp = DoubleProduct(dp + 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 @@ -141,15 +140,18 @@ namespace INTERP_UTILS const double term1 = _coords[5*pt1 + off1] * _coords[5*pt2 + off2]; const double term2 = _coords[5*pt1 + off2] * _coords[5*pt2 + off1]; - +#ifdef FIXED_DELTA + const double delta = FIXED_DELTA; +#else const long double delta = MULT_PREC_F * ( std::abs(term1) + std::abs(term2) ); +#endif - if( std::abs(_doubleProducts[8*seg + dp]) < THRESHOLD_F*delta ) + if( epsilonEqual(_doubleProducts[8*seg + dp], 0.0, THRESHOLD_F * delta)) { if(_doubleProducts[8*seg + dp] != 0.0) { - // std::cout << "Double product for (seg,dp) = (" << seg << ", " << dp << ") = " - // << std::abs(_doubleProducts[8*seg + dp]) << " is imprecise, reset to 0.0" << std::endl; + // std::cout << "Double product for (seg,dp) = (" << seg << ", " << dp << ") = " + // << std::abs(_doubleProducts[8*seg + dp]) << " is imprecise, reset to 0.0" << std::endl; } _doubleProducts[8*seg + dp] = 0.0; @@ -160,6 +162,25 @@ namespace INTERP_UTILS _isDoubleProductsCalculated = true; } + 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]; + + // std::cout << std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner << std::endl; + _doubleProducts[8*seg + dp] = 0.0; + }; + } + double TransformedTriangle::calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const { // NB uses fact that TriSegment <=> TriCorner that is first point of segment (PQ <=> P) @@ -239,10 +260,12 @@ namespace INTERP_UTILS if(minAngle < TRIPLE_PRODUCT_ANGLE_THRESHOLD) { _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, true); + //_tripleProducts[corner] = calcTByDevelopingRow(corner, 1, false); } else { - _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, false); + _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, false); + //_tripleProducts[corner] = calcTByDevelopingRow(corner, 1, false); } _validTP[corner] = true; } @@ -287,7 +310,7 @@ namespace INTERP_UTILS 1.0, 0.0, 0.0, // OX 0.0, 1.0, 0.0, // OY 0.0, 0.0, 1.0, // OZ - -1.0,-1.0, 1.0, // XY + -1.0, 1.0, 0.0, // XY 0.0,-1.0, 1.0, // YZ 1.0, 0.0,-1.0 // ZX }; @@ -297,12 +320,20 @@ namespace INTERP_UTILS EDGE_VECTORS[3*edge + 1], EDGE_VECTORS[3*edge + 2], }; - + + //return angleBetweenVectors(normal, edgeVec); + const double lenNormal = norm(normal); const double lenEdgeVec = norm(edgeVec); const double dotProd = dot(normal, edgeVec); + // return asin( dotProd / ( lenNormal * lenEdgeVec ) ); + //# assert(dotProd / ( lenNormal * lenEdgeVec ) + 1.0 >= 0.0); + //# assert(dotProd / ( lenNormal * lenEdgeVec ) - 1.0 <= 0.0); + //? is this more stable? -> no subtraction + // return asin( dotProd / ( lenNormal * lenEdgeVec ) ) + 3.141592625358979 / 2.0; return 3.14159262535 - acos( dotProd / ( lenNormal * lenEdgeVec ) ); + } /** @@ -456,7 +487,8 @@ namespace INTERP_UTILS const double sumDPProd = coordDPProd[0] + coordDPProd[1] + coordDPProd[2]; const double sumDPProdSq = dot(coordDPProd, coordDPProd); - alpha = sumDPProd / sumDPProdSq; + // alpha = sumDPProd / sumDPProdSq; + alpha = (sumDPProdSq != 0.0) ? sumDPProd / sumDPProdSq : 0.0; } const double cQRbar = cQR * (1.0 - alpha * coordValues[0] * cQR); @@ -464,9 +496,9 @@ namespace INTERP_UTILS const double cPQbar = cPQ * (1.0 - alpha * coordValues[2] * cPQ); // check sign of barred products - should not change - assert(cQRbar * cQR >= 0.0); - assert(cRPbar * cRP >= 0.0); - assert(cPQbar * cPQ >= 0.0); + // assert(cQRbar * cQR >= 0.0); + //assert(cRPbar * cRP >= 0.0); + //assert(cPQbar * cPQ >= 0.0); const double p_term = _coords[5*P + offset] * cQRbar; const double q_term = _coords[5*Q + offset] * cRPbar; @@ -474,10 +506,15 @@ namespace INTERP_UTILS // 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::abs(p_term) + std::abs(q_term) + std::abs(r_term)); - +#endif + if( epsilonEqual( p_term + q_term + r_term, 0.0, THRESHOLD_F * delta) ) { + // std::cout << "Reset imprecise triple product for corner " << corner << " to zero" << std::endl; return 0.0; } else -- 2.39.2