--- /dev/null
+#include "TransformedTriangle.hxx"
+#include <iostream>
+#include <fstream>
+#include <cassert>
+#include <cmath>
+
+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
--- /dev/null
+#include "TransformedTriangle.hxx"
+#include <iostream>
+#include <fstream>
+#include <cassert>
+#include <cmath>
+
+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
