-// Copyright (C) 2007-2016 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2024 CEA, EDF
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
#define __TRANSFORMED_TRIANGLE_HXX__
#include "INTERPKERNELDefines.hxx"
+#include "VectorUtils.hxx"
+#include "assert.h"
#include <vector>
*
*/
- /** \file TransformedTriangle.hxx
- *
+ /**
* OVERVIEW of how the class works : (details can be found in the documentation of each method)
*
* Constructor :
* coordinates of the corners of the triangle. It copies their coordinates and then proceeds to pre-calculating certain
* entities used in the intersection calculation : the double products, triple products and the values of the function E
* (Grandy, [53]).
+ * It is also at this point in constructor that:
+ * - the special case of PQR included in the XYZ plane is treated
+ * - the inconsistencies between double products/triple products computation is handled
*
* calculateIntersectionVolume() :
* This is the only method in the public interface. It calculates the volume under the intersection polygons
* When an intersection point has been detected it is calculated with a corresponding calc* - method in the cases where it
* is not known directly. It is then added to the polygon A and/or B as necessary.
*
- * OPTIMIZE :
- * If OPTIMIZE is defined, a large number of methods will be prefixed with inline and some optimizations concerning the tests
- * with zero double products will be used.
*/
class INTERPKERNEL_EXPORT TransformedTriangle
{
public:
-
- friend class INTERP_TEST::TransformedTriangleTest;
friend class INTERP_TEST::TransformedTriangleIntersectTest;
+ friend class INTERP_TEST::TransformedTriangleTest;
/*
* Enumerations representing the different geometric elements of the unit tetrahedron
* and the triangle. The end element, NO_* gives the number of elements in the enumeration
double calculateIntersectionVolume();
double calculateIntersectionSurface(TetraAffineTransform* tat);
-
void dumpCoords() const;
// Queries of member values used by UnitTetraIntersectionBary
-
const double* getCorner(TriCorner corner) const { return _coords + 5*corner; }
-
const std::vector<double*>& getPolygonA() const { return _polygonA; }
-
double getVolume() const { return _volume; }
protected:
-
TransformedTriangle() { }
// ----------------------------------------------------------------------------------
// High-level methods called directly by calculateIntersectionVolume()
// ----------------------------------------------------------------------------------
void calculateIntersectionAndProjectionPolygons();
-
void calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter);
-
void sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter);
-
double calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter);
// ----------------------------------------------------------------------------------
// High-level methods called directly by calculateIntersectionSurface()
// ----------------------------------------------------------------------------------
void calculateIntersectionPolygon();
-
double calculateSurfacePolygon();
// ----------------------------------------------------------------------------------
// Detection of degenerate triangles
// ----------------------------------------------------------------------------------
-
bool isTriangleInPlaneOfFacet(const TetraFacet facet) const;
-
bool isTriangleParallelToFacet(const TetraFacet facet) const;
-
int isTriangleInclinedToFacet(const TetraFacet facet) const;
-
bool isTriangleBelowTetraeder() const;
// ----------------------------------------------------------------------------------
// Intersection test methods and intersection point calculations
// ----------------------------------------------------------------------------------
-
inline bool testSurfaceEdgeIntersection(const TetraEdge edge) const;
-
void calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const;
-
inline bool testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const;
-
void calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const;
-
bool testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const;
-
void calcIntersectionPtSegmentEdge(const TriSegment seg, const TetraEdge edge, double* pt) const ;
-
bool testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const ;
-
inline bool testSurfaceRayIntersection(const TetraCorner corner) const;
-
bool testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edg);
-
void calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const;
-
bool testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const;
-
inline bool testCornerInTetrahedron(const TriCorner corner) const;
-
inline bool testCornerOnXYZFacet(const TriCorner corner) const;
-
inline bool testCornerAboveXYZFacet(const TriCorner corner) const;
// ----------------------------------------------------------------------------------
// Utility methods used in intersection tests
// ----------------------------------------------------------------------------------
-
bool testTriangleSurroundsEdge(const TetraEdge edge) const;
-
inline bool testEdgeIntersectsTriangle(const TetraEdge edge) const;
-
inline bool testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const;
-
inline bool testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const;
-
bool testSegmentIntersectsHPlane(const TriSegment seg) const;
-
bool testSurfaceAboveCorner(const TetraCorner corner) const;
-
bool testTriangleSurroundsRay(const TetraCorner corner) const;
// ----------------------------------------------------------------------------------
// Double and triple product calculations
// ----------------------------------------------------------------------------------
-
- void resetNearZeroCoordinates();
-
+ void handleDegenerateCases();
bool areDoubleProductsConsistent(const TriSegment seg) const;
-
- void preCalculateDoubleProducts(void);
-
+ void preCalculateDoubleProducts();
inline void resetDoubleProducts(const TriSegment seg, const TetraCorner corner);
-
double calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const;
-
- void preCalculateTripleProducts(void);
-
+ void preCalculateTripleProducts();
double calculateAngleEdgeTriangle(const TetraEdge edge) const;
-
inline double calcStableC(const TriSegment seg, const DoubleProduct dp) const;
-
inline double calcStableT(const TetraCorner corner) const;
+ inline double calcUnstableC(const TriSegment seg, const DoubleProduct dp, double & delta) const;
+ double calcTByDevelopingRow(const TetraCorner corner, const int row, const bool project) const;
- inline double calcUnstableC(const TriSegment seg, const DoubleProduct dp) const;
-
- double calcTByDevelopingRow(const TetraCorner corner, const int row = 1, const bool project = false) const;
+ // ----------------------------------------------------------------------------------
+ // Debug
+ // ----------------------------------------------------------------------------------
+ inline const std::string& strTC(TetraCorner tc) const;
+ inline const std::string& strTE(TetraEdge te) const;
+ inline const std::string& strTF(TetraFacet tf) const;
+ inline const std::string& strTriC(TriCorner tc) const;
+ inline const std::string& strTriS(TriSegment tc) const;
// ----------------------------------------------------------------------------------
// Member variables
/// following order in enumeration DoubleProduct
double _doubleProducts[24];
+ double _deltas[24];
+
/// Array containing the 4 triple products.
/// order : t_O, t_X, t_Y, t_Z
+ /// For example t_O represent the signed volume of the tetrahedron OPQR, and is positive if PQR is oriented clockwise
+ // when seen from the vertex O.
double _tripleProducts[4];
/// Vector holding the points of the intersection polygon A.
// by a given row
static const DoubleProduct DP_FOR_DETERMINANT_EXPANSION[12];
- // values used to decide how imprecise the double products
+ // values used to decide how/when imprecise the double products
// should be to set them to 0.0
- static const long double MACH_EPS; // machine epsilon
- static const long double MULT_PREC_F; // precision of multiplications (Grandy : f)
- static const long double THRESHOLD_F; // threshold for zeroing (Grandy : F/f)
+ static const double MACH_EPS; // machine epsilon
+ static const double MULT_PREC_F; // precision of multiplications (Grandy : f)
+ static const double THRESHOLD_F; // threshold for zeroing (Grandy : F/f)
static const double TRIPLE_PRODUCT_ANGLE_THRESHOLD;
// signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION
static const double SIGN_FOR_SEG_FACET_INTERSECTION[12];
-
+
// coordinates of corners of tetrahedron
static const double COORDS_TET_CORNER[12];
-
+
// indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION
// for the calculation of the coordinates (x,y,z) of the intersection points
// for Segment-Facet and Segment-Edge intersections
};
- // include definitions of inline methods
+ inline void TransformedTriangle::preCalculateTriangleSurroundsEdge()
+ {
+ for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
+ {
+ _triangleSurroundsEdgeCache[edge] = testTriangleSurroundsEdge(edge);
+ }
+ }
+
-#include "TransformedTriangleInline.hxx"
-}
+ // ----------------------------------------------------------------------------------
+ // TransformedTriangle_math.cxx
+ // ----------------------------------------------------------------------------------
+
+ inline void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner)
+ {
+ // set the three corresponding double products to 0.0
+ static const DoubleProduct DOUBLE_PRODUCTS[12] =
+ {
+ C_YZ, C_ZX, C_XY, // O
+ C_YZ, C_ZH, C_YH, // X
+ C_ZX, C_ZH, C_XH, // Y
+ C_XY, C_YH, C_XH // Z
+ };
+
+ for(int i = 0 ; i < 3 ; ++i) {
+ const DoubleProduct dp = DOUBLE_PRODUCTS[3*corner + i];
+
+ LOG(6, std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner);
+ _doubleProducts[8*seg + dp] = 0.0;
+ };
+ }
+
+ inline double TransformedTriangle::calcStableC(const TriSegment seg, const DoubleProduct dp) const
+ {
+ return _doubleProducts[8*seg + dp];
+ }
+
+ inline double TransformedTriangle::calcStableT(const TetraCorner corner) const
+ {
+ assert(_validTP[corner]);
+ return _tripleProducts[corner];
+ }
+
+ inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp, double& delta) 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];
+
+ const double prd1 = _coords[5*pt1 + off1] * _coords[5*pt2 + off2],
+ prd2 = _coords[5*pt1 + off2] * _coords[5*pt2 + off1];
+ delta = std::fabs(prd1) + std::fabs(prd2);
+ return prd1 - prd2;
+ }
+
+ // ----------------------------------------------------------------------------------
+ // TransformedTriangle_intersect.cxx
+ // ----------------------------------------------------------------------------------
+ inline bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const
+ {
+ return _triangleSurroundsEdgeCache[edge] && testEdgeIntersectsTriangle(edge);
+ }
+
+ inline bool TransformedTriangle::testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const
+ {
+ return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet);
+ }
+
+ inline bool TransformedTriangle::testSurfaceRayIntersection(const TetraCorner corner) const
+ {
+ return testTriangleSurroundsRay( corner ) && testSurfaceAboveCorner( corner );
+ }
+
+ inline bool TransformedTriangle::testCornerInTetrahedron(const TriCorner corner) const
+ {
+ const double pt[4] =
+ {
+ _coords[5*corner], // x
+ _coords[5*corner + 1], // y
+ _coords[5*corner + 2], // z
+ _coords[5*corner + 3] // z
+ };
+
+ for(int i = 0 ; i < 4 ; ++i)
+ {
+ if(pt[i] < 0.0 || pt[i] > 1.0)
+ {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ inline bool TransformedTriangle::testCornerOnXYZFacet(const TriCorner corner) const
+ {
+#if 0
+ const double pt[4] =
+ {
+ _coords[5*corner], // x
+ _coords[5*corner + 1], // y
+ _coords[5*corner + 2], // z
+ _coords[5*corner + 3] // h
+ };
+#endif
+ const double* pt = &_coords[5*corner];
+
+ if(pt[3] != 0.0)
+ {
+ return false;
+ }
+
+ for(int i = 0 ; i < 3 ; ++i)
+ {
+ if(pt[i] < 0.0 || pt[i] > 1.0)
+ {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ inline bool TransformedTriangle::testCornerAboveXYZFacet(const TriCorner corner) const
+ {
+ const double x = _coords[5*corner];
+ const double y = _coords[5*corner + 1];
+ const double h = _coords[5*corner + 3];
+ const double H = _coords[5*corner + 4];
+
+ return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0;
+
+ }
+
+ inline bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const
+ {
+ // correspondence edge - triple products for edges OX, ..., ZX (Grandy, table III)
+ static const TetraCorner TRIPLE_PRODUCTS[12] =
+ {
+ X, O, // OX
+ Y, O, // OY
+ Z, O, // OZ
+ X, Y, // XY
+ Y, Z, // YZ
+ Z, X, // ZX
+ };
+
+ // Grandy, [16]
+ const double t1 = calcStableT(TRIPLE_PRODUCTS[2*edge]);
+ const double t2 = calcStableT(TRIPLE_PRODUCTS[2*edge + 1]);
+
+ // [ABN] Okayyy: if either t1 or t2 exactly equal zero, then it can mean two things:
+ // - either PQR is very close to the corner -> this is OK, further computation of intersection point between
+ // surface and edge will produce a correct result
+ // - or, if the other triple prod is also very small, then this is a degenerate case: the edge is almost in PQR -> this is bad
+ // and leads to weird intersection point computation -> we avoid this.
+ // PS : here was written "// tuleap26461" -> whoo this helps :-)
+ if (t1 == 0.0 || t2 == 0.0)
+ if (std::fabs(t1+t2) < THRESHOLD_F*MULT_PREC_F)
+ return false;
+
+ return (t1*t2 <= 0.0) && !epsilonEqual(t1,t2, MULT_PREC_F);
+ }
+
+ inline bool TransformedTriangle::testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const
+ {
+#if 0
+ const double signs[3] =
+ {
+ SIGN_FOR_SEG_FACET_INTERSECTION[3*facet],
+ SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 1],
+ SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 2]
+ };
+#endif
+ const double* signs = &SIGN_FOR_SEG_FACET_INTERSECTION[3*facet];
+ const double c1 = signs[0]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet]);
+ const double c2 = signs[1]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]);
+ const double c3 = signs[2]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]);
+
+ return (c1*c3 > 0.0) && (c2*c3 > 0.0);
+ }
+
+ inline bool TransformedTriangle::testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const
+ {
+ // use correspondence facet a = 0 <=> offset for coordinate a in _coords
+ // and also correspondence segment AB => corner A
+ const double coord1 = _coords[5*seg + facet];
+ const double coord2 = _coords[5*( (seg + 1) % 3) + facet];
+
+ //? should we use epsilon-equality here in second test?
+ LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 );
+
+ return (coord1*coord2 <= 0.0) && (coord1 != coord2);
+ }
+
+ inline bool TransformedTriangle::testSegmentIntersectsHPlane(const TriSegment seg) const
+ {
+ // get the H - coordinates
+ const double coord1 = _coords[5*seg + 4];
+ const double coord2 = _coords[5*( (seg + 1) % 3) + 4];
+ //? should we use epsilon-equality here in second test?
+ LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 );
+
+ return (coord1*coord2 <= 0.0) && (coord1 != coord2);
+ }
+
+ inline bool TransformedTriangle::testSurfaceAboveCorner(const TetraCorner corner) const
+ {
+ // There is an error in Grandy -> it should be C_XY instead of C_YZ in [28].
+ //
+ // Idea: the nz value (Grandy [28] corrected!) can be interpreted as a special variant of the triple product t_O
+ // where the z coordinate has been set to 1. It represents the signed volume of the tet OP'Q'R' where P', Q' and R' are the
+ // projection of P, Q, R on the z=1 plane.
+ // Comparing the sign of this triple product with t_X (or t_Y, ... dep on the corner) indicates whether the corner is in the
+ // half-space above or below the PQR triangle, similarly to what is explained in [16].
+ // (this works even for the Z corner, since we could have chosen z=24 (instead of z=1) this would not have changed the final sign test).
+ const double nz = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY);
+
+ // Triple product might have not been computed, but here we need one:
+ const double tp = _validTP[corner] ? calcStableT(corner) : calcTByDevelopingRow(corner, 1, false);
+
+ return tp*nz >= 0.0;
+ }
+
+ inline bool TransformedTriangle::testTriangleSurroundsRay(const TetraCorner corner) const
+ {
+ // double products to use for the possible corners
+ static const DoubleProduct DP_FOR_RAY_INTERSECTION[4] =
+ {
+ DoubleProduct(0), // O - only here to fill out and make indices match
+ C_10, // X
+ C_01, // Y
+ C_XY // Z
+ };
+
+ const DoubleProduct dp = DP_FOR_RAY_INTERSECTION[corner];
+
+ const double cPQ = calcStableC(PQ, dp);
+ const double cQR = calcStableC(QR, dp);
+ const double cRP = calcStableC(RP, dp);
+
+ return ( cPQ*cQR > 0.0 ) && ( cPQ*cRP > 0.0 );
+ }
+
+ inline const std::string& TransformedTriangle::strTC(TetraCorner tc) const
+ {
+ static const std::map<TetraCorner, std::string> m = {{O, "O"}, {X, "X"}, {Y, "Y"}, {Z, "Z"}};
+ return m.at(tc);
+ }
+
+ inline const std::string& TransformedTriangle::strTE(TetraEdge te) const
+ {
+ static const std::map<TetraEdge, std::string> m = {{OX, "OX"}, {OY, "OY"}, {OZ, "OZ"}, {XY, "XY"}, {YZ, "YZ"},
+ {ZX, "ZX"}, {H01, "H01"}, {H10, "H10"}};
+ return m.at(te);
+ }
+
+ inline const std::string& TransformedTriangle::strTF(TetraFacet tf) const
+ {
+ static const std::map<TetraFacet, std::string> m = {{OYZ, "OYZ"}, {OZX, "OZX"}, {OXY, "OXY"}, {XYZ, "XYZ"}};
+ return m.at(tf);
+ }
+
+ inline const std::string& TransformedTriangle::strTriC(TriCorner tc) const
+ {
+ static const std::map<TriCorner, std::string> m = {{P, "P"}, {Q, "Q"}, {R, "R"}};
+ return m.at(tc);
+ }
+
+ inline const std::string& TransformedTriangle::strTriS(TriSegment ts) const
+ {
+ static const std::map<TriSegment, std::string> m = {{PQ, "PQ"}, {QR, "QR"}, {RP, "RP"}};
+ return m.at(ts);
+ }
+
+}
#endif
