From b719896a4b47d4d333da259dc2bf033a616dff2d Mon Sep 17 00:00:00 2001 From: vbd Date: Thu, 5 Jul 2007 14:35:18 +0000 Subject: [PATCH] ronnas: Adding the other two implementation files for the class TransformedTriangle --- .../TransformedTriangle_intersect.cxx | 457 +++++++++++++++++ .../TransformedTriangle_math.cxx | 467 ++++++++++++++++++ 2 files changed, 924 insertions(+) create mode 100644 src/INTERP_KERNEL/TransformedTriangle_intersect.cxx create mode 100644 src/INTERP_KERNEL/TransformedTriangle_math.cxx diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx new file mode 100644 index 000000000..007ec67c6 --- /dev/null +++ b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx @@ -0,0 +1,457 @@ +#include "TransformedTriangle.hxx" +#include +#include +#include +#include + +namespace INTERP_UTILS +{ + + //////////////////////////////////////////////////////////////////////////////////// + /// Constants ///////////////// + //////////////////////////////////////////////////////////////////////////////////// + + // correspondance facet - double product + // Grandy, table IV + const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_SEG_FACET_INTERSECTION[12] = + { + C_XH, C_XY, C_ZX, // OYZ + C_YH, C_YZ, C_XY, // OZX + C_ZH, C_ZX, C_YZ, // OXY + C_XH, C_YH, C_ZH // XYZ + }; + // signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION + const double TransformedTriangle::SIGN_FOR_SEG_FACET_INTERSECTION[12] = + { + 1.0, 1.0, -1.0, + 1.0, 1.0, -1.0, + 1.0, 1.0, -1.0, + 1.0, 1.0, 1.0 + }; + + + //////////////////////////////////////////////////////////////////////////////////// + /// Intersection test methods and intersection point calculations //////// + //////////////////////////////////////////////////////////////////////////////////// + /** + * Tests if the given edge of the tetrahedron intersects the triangle PQR. (Grandy, eq [17]) + * + * @param edge edge of tetrahedron + * @returns true if PQR intersects the edge, and the edge is not in the plane of the triangle. + */ + bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const + { + return testTriangleSurroundsEdge(edge) && testEdgeIntersectsTriangle(edge); + } + + /** + * Calculates the point of intersection between the given edge of the tetrahedron and the + * triangle PQR. (Grandy, eq [22]) + * + * @pre testSurfaceEdgeIntersection(edge) returns true + * @param edge edge of tetrahedron + * @param pt array of three doubles in which to store the coordinates of the intersection point + */ + void TransformedTriangle::calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const + { + // not implemented + } + + /** + * Tests if the given segment of the triangle intersects the given facet of the tetrahedron. + * (Grandy, eq. [19]) + * + * @param seg segment of the triangle + * @param facet facet of the tetrahedron + * @returns true if the segment intersects the facet + */ + bool TransformedTriangle::testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const + { + return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet); + } + + /** + * Calculates the point of intersection between the given segment of the triangle + * and the given facet of the tetrahedron. (Grandy, eq. [23]) + * + * @pre testSurfaceEdgeIntersection(seg, facet) returns true + * + * @param seg segment of the triangle + * @param facet facet of the tetrahedron + * @param pt array of three doubles in which to store the coordinates of the intersection point + */ + void TransformedTriangle::calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const + { + // not implemented + } + + /** + * Tests if the given segment of the triangle intersects the given edge of the tetrahedron (Grandy, eq. [20] + * + * @param seg segment of the triangle + * @param edge edge of tetrahedron + * @returns true if the segment intersects the edge + */ + bool TransformedTriangle::testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const + { + // facets shared by each edge + static const TetraFacet FACET_FOR_EDGE[12] = + { + OXY, OZX, // OX + OXY, OYZ, // OY + OZX, OYZ, // OZ + OXY, XYZ, // XY + OYZ, XYZ, // YZ + OZX, XYZ // ZX + }; + + if(calcStableC(seg,DoubleProduct( edge )) != 0.0) + { + return false; + } + else + { + // check condition that the double products for one of the two + // facets adjacent to the edge has a positive product + bool isFacetCondVerified = false; + TetraFacet facet[2]; + 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; + dp1 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1]; + } + else if(dp2 == DoubleProduct( edge )) + { + 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); + + isFacetCondVerified = isFacetCondVerified || c1*c2 > 0.0; + } + if(!isFacetCondVerified) + { + return false; + } + else + { + return testSegmentIntersectsFacet(seg, facet[0]) || testSegmentIntersectsFacet(seg, facet[1]); + } + } + } + + /** + * Calculates the point of intersection between the given segment of the triangle + * and the given edge of the tetrahedron. (Grandy, eq. [25]) + * + * @pre testSegmentEdgeIntersection(seg, edge) returns true + * + * @param seg segment of the triangle + * @param edge edge of the tetrahedron + * @param pt array of three doubles in which to store the coordinates of the intersection point + */ + void TransformedTriangle::calcIntersectionPtSegmentEdge(const TriSegment seg, const TetraEdge edge, double* pt) const + { + // not implemented + } + + /** + * Tests if the given segment of the triangle intersects the given corner of the tetrahedron. + * (Grandy, eq. [21]) + * + * @param seg segment of the triangle + * @param corner corner of the tetrahedron + * @returns true if the segment intersects the corner + */ + bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const + { + // edges meeting at a given corner + static const TetraEdge EDGES_FOR_CORNER[12] = + { + OX, OY, OZ, // O + OX, XY, ZX, // X + OY, XY, YZ, // Y + OZ, ZX, YZ // Z + }; + + // facets meeting at a given corner + static const TetraFacet FACETS_FOR_CORNER[12] = + { + OXY, OYZ, OZX, // O + OZX, OXY, XYZ, // X + OYZ, XYZ, OXY, // Y + OZX, XYZ, OYZ // Z + }; + + // check double products are zero + for(int i = 0 ; i < 3 ; ++i) + { + const TetraEdge edge = EDGES_FOR_CORNER[3*corner + i]; + const DoubleProduct dp = DoubleProduct( edge ); + const double c = calcStableC(seg, dp); + if(c != 0.0) + { + return false; + } + } + + // check segment intersect a facet + for(int i = 0 ; i < 3 ; ++i) + { + const TetraFacet facet = FACETS_FOR_CORNER[3*corner + i]; + if(testSegmentIntersectsFacet(seg, facet)) + { + return true; + } + } + + return false; + } + + + /** + * Tests if the triangle PQR intersects the ray pointing in the upwards z - direction + * from the given corner of the tetrahedron. (Grandy eq. [29]) + * + * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z) + * @returns true if the upwards ray from the corner intersects the triangle. + */ + bool TransformedTriangle::testSurfaceRayIntersection(const TetraCorner corner) const + { + // not implemented + //return testTriangleSurroundsEdge(H_10,etc) && testSurfaceAboveCorner(corner); + return false; + } + + /** + * Tests if the given segment of the triangle intersects the half-strip above the + * given edge of the h = 0 plane. (Grandy, eq. [30]) + * + * @param seg segment of the triangle + * @param edge edge of the h = 0 plane of the tetrahedron (XY, YZ, ZX) + * @returns true if the upwards ray from the corner intersects the triangle. + */ + bool TransformedTriangle::testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edg) + { + // not implemented + return false; + } + + /** + * Calculates the point of intersection between the given segment of the triangle + * and the halfstrip above the given edge of the tetrahedron. (Grandy, eq. [31]) + * + * @pre testSegmentHalfstripIntersection(seg, edge) returns true + * + * @param seg segment of the triangle + * @param edge edge of the tetrahedron defining the halfstrip + * @param pt array of three doubles in which to store the coordinates of the intersection point + */ + void TransformedTriangle::calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const + { + // not implemented + } + + /** + * 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]) + * + * @param seg segment of the triangle PQR + * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z) + * @returns true if the upwards ray from the corner intersects the segment. + */ + bool TransformedTriangle::testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const + { + // not implemented + return false; + } + + /** + * Tests if the given corner of triangle PQR lies in the interior of the unit tetrahedron + * (0 <= x_p, y_p, z_p, h_p <= 1) + * + * @param corner corner of the triangle PQR + * @returns true if the corner lies in the interior of the unit tetrahedron. + */ + 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; + } + + /** + * Tests if the given corner of triangle PQR lies on the facet h = 0 (the XYZ facet) + * (0 <= x_p, y_p, z_p <= 1 && h_p = 0) + * + * @param corner corner of the triangle PQR + * @returns true if the corner lies on the facet h = 0 + */ + bool TransformedTriangle::testCornerOnXYZFacet(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] // h + }; + + 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; + } + + + //////////////////////////////////////////////////////////////////////////////////// + /// Utility methods used in intersection tests /////////////// + //////////////////////////////////////////////////////////////////////////////////// + /** + * Tests if the triangle PQR surrounds the axis on which the + * given edge of the tetrahedron lies. + * + * @param edge edge of tetrahedron + * @returns true if PQR surrounds edge, false if not (see Grandy, eq. [53]) + */ + bool TransformedTriangle::testTriangleSurroundsEdge(const TetraEdge edge) const + { + // NB DoubleProduct enum corresponds to TetraEdge enum according to Grandy, table III + // so we can use the edge directly + + // optimization : we could use _doubleProducts directly here + const double cPQ = calcStableC(TransformedTriangle::PQ, DoubleProduct(edge)); + const double cQR = calcStableC(TransformedTriangle::QR, DoubleProduct(edge)); + const double cRP = calcStableC(TransformedTriangle::RP, DoubleProduct(edge)); + + // 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); + + return (cPQ*cQR >= 0.0) && (cQR*cRP >=0.0) && (cRP*cPQ >= 0.0) && numZeros < 2; + } + + /** + * Tests if the corners of the given edge lie on different sides of the triangle PQR. + * + * @param edge edge of the tetrahedron + * @returns true if the corners of the given edge lie on different sides of the triangle PQR + * or if one (but not both) lies in the plane of the triangle. + */ + bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const + { + + assert(edge < H01); + + // correspondance 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 + Y, Z, // YZ + Z, X, // ZX + X, Y // XY + }; + + // 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? + return (t1*t2 <= 0.0) && (t1 - t2 != 0.0); + } + + /** + * Tests if the given facet of the tetrahedron surrounds the axis on which the + * given segment of the triangle lies. + * + * @param seg segment of the triangle + * @param facet facet of the tetrahedron + * @returns true if the facet surrounds the segment, false if not (see Grandy, eq. [18]) + */ + bool TransformedTriangle::testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const + { + + 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] + }; + + 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); + } + + /** + * Tests if the corners of the given segment lie on different sides of the given facet. + * + * @param seg segment of the triangle + * @param facet facet of the tetrahedron + * @returns true if the corners of the given segment lie on different sides of the facet + * or if one (but not both) lies in the plane of the facet. (see Grandy, eq. [18]) + */ + bool TransformedTriangle::testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const + { + // use correspondance facet a = 0 <=> offset for coordinate a in _coords + // and also coorespondance 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? + return (coord1*coord2 <= 0.0) && (coord1 != coord2); + } + + /** + * Tests if the triangle PQR lies above a given corner in the z-direction (implying that the + * ray pointing upward in the z-direction from the corner can intersect the triangle) + * (Grandy eq.[28]) + * + * @param corner corner of the tetrahedron + * @returns true if the triangle lies above the corner in the z-direction + */ + bool TransformedTriangle::testSurfaceAboveCorner(const TetraCorner corner) + { + // not implemented + return false; + } + +}; // NAMESPACE diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx new file mode 100644 index 000000000..145cf93ed --- /dev/null +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -0,0 +1,467 @@ +#include "TransformedTriangle.hxx" +#include +#include +#include +#include + +namespace INTERP_UTILS +{ + + //////////////////////////////////////////////////////////////////////////////////// + /// Constants ///////////////// + //////////////////////////////////////////////////////////////////////////////////// + const int TransformedTriangle::DP_OFFSET_1[8] = {1, 2, 0, 2, 0, 1, 4, 1}; + const int TransformedTriangle::DP_OFFSET_2[8] = {2, 0, 1, 3, 3, 3, 0, 4}; + + const int TransformedTriangle::COORDINATE_FOR_DETERMINANT_EXPANSION[12] = + { + 0, 1, 2, // O + 3, 1, 2, // X + 0, 3, 2, // Y + 0, 1, 3 // Z + }; + + const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_DETERMINANT_EXPANSION[12] = + { + C_YZ, C_ZX, C_XY, // O + C_YZ, C_ZH, C_YH, // X + C_ZH, C_ZX, C_XH, // Y + C_YH, C_XH, C_XY // Z + }; + + const double TransformedTriangle::MACH_EPS = 1.0e-15; + const double TransformedTriangle::MULT_PREC_F = 4.0*TransformedTriangle::MACH_EPS; + const double TransformedTriangle::THRESHOLD_F = 20.0; + + const double TransformedTriangle::TRIPLE_PRODUCT_ANGLE_THRESHOLD = 0.1; + + //////////////////////////////////////////////////////////////////////////////////// + /// Double and triple product calculations /////////////// + //////////////////////////////////////////////////////////////////////////////////// + + /** + * Pre-calculates all double products for this triangle, and stores + * them internally. This method makes compensation for precision errors, + * and it is thus the "stable" double products that are stored. + * + */ + void TransformedTriangle::preCalculateDoubleProducts(void) + { + if(_isDoubleProductsCalculated) + return; + + // aliases to shorten code + typedef TransformedTriangle::DoubleProduct DoubleProduct; + typedef TransformedTriangle::TetraCorner TetraCorner; + typedef TransformedTriangle::TriSegment TriSegment; + typedef TransformedTriangle TT; + + // -- calculate all unstable double products -- store in _doubleProducts + for(TriSegment seg = TT::PQ ; seg <= TT::RP ; seg = TriSegment(seg + 1)) + { + for(DoubleProduct dp = TT::C_YZ ; dp <= TT::C_10 ; dp = DoubleProduct(dp + 1)) + { + _doubleProducts[8*seg + dp] = calcUnstableC(seg, dp); + } + } + + + // -- (1) for each segment : check that double products satisfy Grandy, [46] + // -- and make corrections if not + for(TriSegment seg = TT::PQ ; seg <= TT::RP ; seg = TriSegment(seg + 1)) + { + const double term1 = _doubleProducts[8*seg + TT::C_YZ] * _doubleProducts[8*seg + TT::C_XH]; + const double term2 = _doubleProducts[8*seg + TT::C_ZX] * _doubleProducts[8*seg + TT::C_YH]; + const double term3 = _doubleProducts[8*seg + TT::C_XY] * _doubleProducts[8*seg + TT::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; + + // 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); + const int num_neg = (term1 < 0.0 ? 1 : 0) + (term2 < 0.0 ? 1 : 0) + (term3 < 0.0 ? 1 : 0); + + + if(num_zero == 2 || (num_zero !=3 && num_neg == 0) || (num_zero !=3 && num_neg == 3)) + { + std::cout << "inconsistent! "<< std::endl; + + // find TetraCorner nearest to segment + double min_dist; + TetraCorner min_corner; + + for(TetraCorner corner = TT::O ; corner <= TT::Z ; corner = TetraCorner(corner + 1)) + { + // calculate distance corner - segment axis + // NB uses fact that TriSegment <=> TriCorner that is first point of segment (PQ <=> P) + const TriCorner ptP_idx = TriCorner(seg); + const TriCorner ptQ_idx = TriCorner(seg + 1 % 3); + + const double ptP[3] = { _coords[5*ptP_idx], _coords[5*ptP_idx + 1], _coords[5*ptP_idx + 2] }; + const double ptQ[3] = { _coords[5*ptQ_idx], _coords[5*ptQ_idx + 1], _coords[5*ptQ_idx + 2] }; + + // coordinates of corner + const double ptTetCorner[3] = + { + corner == TT::X ? 1.0 : 0.0, + corner == TT::Y ? 1.0 : 0.0, + corner == TT::Z ? 1.0 : 0.0 + }; + + // dist^2 = ( PQ x CP )^2 / |PQ|^2 where C is the corner point + + // difference vectors + const double diffPQ[3] = { ptQ[0] - ptP[0], ptQ[1] - ptP[1], ptQ[2] - ptP[2] }; + const double diffCornerP[3] = { ptP[0] - ptTetCorner[0], ptP[1] - ptTetCorner[1], ptP[2] - ptTetCorner[2] }; + + // 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], + }; + + + 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 dist = cross_squared / norm_diffPQ_squared; + + // update minimum (initializs with first corner) + if(corner == TT::O || dist < min_dist) + { + min_dist = dist; + min_corner = corner; + } + } + + // set the three corresponding double products to 0.0 + static const DoubleProduct DOUBLE_PRODUCTS[12] = + { + TT::C_YZ, TT::C_ZX, TT::C_XY, // O + TT::C_YZ, TT::C_ZH, TT::C_YH, // X + TT::C_ZX, TT::C_ZH, TT::C_XH, // Y + TT::C_XY, TT::C_YH, TT::C_XH // Z + }; + + for(int i = 0 ; i < 3 ; ++i) { + DoubleProduct dp = DOUBLE_PRODUCTS[3*min_corner + i]; + + // std::cout << std::endl << "in code inconsistent dp :" << dp << std::endl; + + _doubleProducts[8*seg + dp] = 0.0; + } + + } + } + + + // -- (2) check that each double product statisfies Grandy, [47], else set to 0 + + for(TriSegment seg = TT::PQ ; seg <= TT::RP ; seg = TriSegment(seg + 1)) + { + for(DoubleProduct dp = TT::C_YZ ; dp <= TT::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 double term1 = _coords[5*pt1 + off1] * _coords[5*pt2 + off2]; + const double term2 = _coords[5*pt1 + off2] * _coords[5*pt2 + off1]; + const double absTerm1 = (term1 > 0.0 ? term1 : -term1); + const double absTerm2 = (term2 > 0.0 ? term2 : -term2); + + const double delta = MULT_PREC_F*(absTerm1 + absTerm2); + const double absDoubleProduct = _doubleProducts[8*seg + dp] > 0.0 ? _doubleProducts[8*seg + dp] : -_doubleProducts[8*seg + dp]; + + if(absDoubleProduct < THRESHOLD_F*delta) + { + std::cout << "Double product " << 8*seg+dp << " = " << absDoubleProduct << " is imprecise, reset to 0.0" << std::endl; + _doubleProducts[8*seg + dp] = 0.0; + } + } + } + + _isDoubleProductsCalculated = true; + } + + /** + * Pre-calculates all triple products for the tetrahedron with respect to + * this triangle, and stores them internally. This method takes into account + * the problem of errors due to cancellation. + * + */ + void TransformedTriangle::preCalculateTripleProducts(void) + { + if(_isTripleProductsCalculated) + return; + + for(TetraCorner corner = O ; corner <= Z ; corner = TetraCorner(corner + 1)) + { + bool isGoodRowFound = false; + + // find edge / row to use + int minRow; + double minAngle; + bool isMinInitialised = false; + + for(int row = 1 ; row < 4 ; ++row) + { + DoubleProduct dp = DP_FOR_DETERMINANT_EXPANSION[3*corner + (row - 1)]; + + // get edge by using correspondance between Double Product and Edge + TetraEdge edge = TetraEdge(dp); + + // use edge only if it is surrounded by the surface + if( testTriangleSurroundsEdge(edge) ) + { + isGoodRowFound = true; + + // -- calculate angle between edge and PQR + // find normal to PQR - cross PQ and PR + const double pq[3] = + { + _coords[5*Q] - _coords[5*P], + _coords[5*Q + 1] - _coords[5*P + 1], + _coords[5*Q + 2] - _coords[5*P + 2] + }; + + const double pr[3] = + { + _coords[5*R] - _coords[5*P], + _coords[5*R + 1] - _coords[5*P + 1], + _coords[5*R + 2] - _coords[5*P + 2] + }; + + const double normal[3] = + { + pq[1]*pr[2] - pq[2]*pr[1], + pq[2]*pr[0] - pq[0]*pr[2], + pq[0]*pr[1] - pq[1]*pr[0] + }; + + static const double EDGE_VECTORS[18] = + { + 1.0, 0.0, 0.0, // OX + 0.0, 1.0, 0.0, // OY + 0.0, 0.0, 1.0, // OZ + 0.0,-1.0, 1.0, // YZ + 1.0, 0.0,-1.0, // ZX + -1.0,-1.0, 1.0 // XY + }; + + const double edgeVec[3] = { + EDGE_VECTORS[3*edge], + EDGE_VECTORS[3*edge + 1], + EDGE_VECTORS[3*edge + 2], + }; + + const double lenNormal = sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); + const double lenEdgeVec = sqrt(edgeVec[0]*edgeVec[0] + edgeVec[1]*edgeVec[1] + edgeVec[2]*edgeVec[2]); + const double dotProd = normal[0]*edgeVec[0] + normal[1]*edgeVec[1] + normal[2]*edgeVec[2]; + + const double angle = 3.14159262535 - acos( dotProd / ( lenNormal * lenEdgeVec ) ); + + if(!isMinInitialised || angle < minAngle) + { + minAngle = angle; + minRow = row; + isMinInitialised = true; + } + + } + } + + if(isGoodRowFound) + { + if(minAngle < TRIPLE_PRODUCT_ANGLE_THRESHOLD) + { + _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, true); + } + else + { + _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, false); + } + } + else + { + // this value will not be used + // we set it to whatever + _tripleProducts[corner] = -3.14159265; + } + + } + + _isTripleProductsCalculated = true; + } + + /** + * Returns the stable double product c_{xy}^{pq} + * + * @pre The stable double products have been calculated with preCalculateDoubleProducts. + * @param seg segment of triangle + * @param dp double product sought + * + * @returns stabilised double product c_{xy}^{pq} + * + */ + double TransformedTriangle::calcStableC(const TriSegment seg, const DoubleProduct dp) const + { + assert(_isDoubleProductsCalculated); + return _doubleProducts[8*seg + dp]; + } + + + /** + * Returns the stable triple product t_X for a given corner + * The triple product gives the signed volume of the tetrahedron between + * this corner and the triangle PQR. These triple products have been calculated + * in a way to avoid problems with cancellation. + * + * @pre double products have already been calculated + * @pre triple products have already been calculated + * @param corner corner for which the triple product is calculated + * @returns triple product associated with corner (see Grandy, eqs. [50]-[52]) + */ + double TransformedTriangle::calcStableT(const TetraCorner corner) const + { + assert(_isTripleProductsCalculated); + return _tripleProducts[corner]; + } + + /** + * Calculates the given double product c_{xy}^{pq} = x_p*y_q - y_p*x_q for a + * a segment PQ of the triangle. This method does not compensate for + * precision errors. + * + * @param seg segment of triangle + * @param dp double product sought + * + * @returns double product c_{xy}^{pq} + * + */ + 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]; + } + + /** + * 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. + * + * @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 + * @param project indicates whether or not to perform projection as inidicated in Grandy, p.446 + * @returns triple product associated with corner (see Grandy, [50]-[52]) + */ + + double TransformedTriangle::calcTByDevelopingRow(const TetraCorner corner, const int row, const bool project) const + { + + // OVERVIEW OF CALCULATION + // --- sign before the determinant + // the sign used depends on the sign in front of the triple product (Grandy, [15]), + // and the convention used in the definition of the double products + + // the sign in front of the determinant gives the following schema for the three terms (I): + // corner/row 1 2 3 + // O (sign:+) + - + + // X (sign:-) - + - + // Y (sign:-) - + - + // Z (sign:-) - + - + + // the 2x2 determinants are the following (C_AB <=> A_p*B_q - B_p*A_q, etc) + // corner/row 1 2 3 + // O (sign:+) C_YZ C_XZ C_XY + // X (sign:-) C_YZ C_HZ C_HY + // Y (sign:-) C_HZ C_XZ C_XH + // Z (sign:-) C_YH C_XH C_XY + + // these are represented in DP_FOR_DETERMINANT_EXPANSION, + // except for the fact that certain double products are inversed (C_AB <-> C_BA) + + // comparing with the DOUBLE_PRODUCTS and using the fact that C_AB = -C_BA + // we deduce the following schema (II) : + // corner/row 1 2 3 + // O (sign:+) + - + + // X (sign:-) + - - + // Y (sign:-) - - + + // Z (sign:-) + + + + + // comparing the two schemas (I) and (II) gives us the following matrix of signs, + // putting 1 when the signs in (I) and (II) are equal and -1 when they are different : + + static const int SIGNS[12] = + { + 1, 1, 1, + -1,-1, 1, + 1,-1,-1, + -1, 1,-1 + }; + + // find the offsets of the rows of the determinant + const int offset = COORDINATE_FOR_DETERMINANT_EXPANSION[3 * corner + (row - 1)]; + + const DoubleProduct dp = DP_FOR_DETERMINANT_EXPANSION[3 * corner + (row - 1)]; + + const int sign = SIGNS[3 * corner + (row - 1)]; + + const double cQR = calcStableC(QR, dp); + const double cRP = calcStableC(RP, dp); + const double cPQ = calcStableC(PQ, dp); + + // coordinate to use for projection (Grandy, [57]) with edges + // OX, OY, OZ, YZ, ZX, XY 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 } ; + + double alpha = 0.0; + + 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[0] * cPQ }; + + const double sumDPProd = coordDPProd[0] + coordDPProd[1] + coordDPProd[2]; + const double sumDPProdSq = coordDPProd[0]*coordDPProd[0] + coordDPProd[1]*coordDPProd[1] + coordDPProd[2]*coordDPProd[2]; + alpha = sumDPProd / sumDPProdSq; + } + + 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 ); + + } + +}; // NAMESPACE -- 2.39.2