+ triple product inconsistency was not properly detected (now using original deltas from douple prod computation)
+ f/F factor set to 20 (not 500) as in Grandy's article -> (will be set to 100, see next commit)
+ epsilonEqual used when necessary: when it provides more points in polygon A or B (risk is completely missing a point in polygon!)
+ better handling of degenerated case where PQR triangle is in h=0
plane, or when P,Q or R close to a tet corner
+ fixed surface-edge intersection test: triple product equality (and zeroing) must be checked with care
+ using 'long double' is not necessary -> double are enough
_coords[5*Q + 3] = 1 - q[0] - q[1] - q[2];
_coords[5*R + 3] = 1 - r[0] - r[1] - r[2];
+ // Handle degenerated case where one of the seg of PQR is (almost) inside XYZ plane,
+ // and hence by extension when the whole PQR triangle is in the XYZ plane
+ handleDegenerateCases();
+
// 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];
- resetNearZeroCoordinates();
-
// initialise rest of data
preCalculateDoubleProducts();
{
// coordinate to check
const int coord = static_cast<int>(facet);
- return (_coords[5*P + coord] == _coords[5*Q + coord]) && (_coords[5*P + coord] == _coords[5*R + coord]);
+ return (epsilonEqual(_coords[5*P + coord], _coords[5*Q + coord])) && (epsilonEqual(_coords[5*P + coord], _coords[5*R + coord]));
}
/**
* 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
return _tripleProducts[corner];
}
- inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp) const
+ 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 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];
+ 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;
}
// ----------------------------------------------------------------------------------
inline bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const
{
- // correspondence edge - triple products
- // for edges OX, ..., ZX (Grandy, table III)
+ // correspondence edge - triple products for edges OX, ..., ZX (Grandy, table III)
static const TetraCorner TRIPLE_PRODUCTS[12] =
{
X, O, // OX
const double t1 = calcStableT(TRIPLE_PRODUCTS[2*edge]);
const double t2 = calcStableT(TRIPLE_PRODUCTS[2*edge + 1]);
- //? should equality with zero use epsilon?
- LOG(5, "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 );
- return (t1*t2 <= 0.0) && !epsilonEqual(t1 - t2, 0.0); // tuleap26461
+ // [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, (double)MULT_PREC_F);
}
inline bool TransformedTriangle::testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const
inline bool TransformedTriangle::testSurfaceAboveCorner(const TetraCorner corner) const
{
- // ? There seems to be an error in Grandy -> it should be C_XY instead of C_YZ in [28].
- // ? I haven't really figured out why, but it seems to work.
- const double normal = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY);
-
- LOG(6, "surface above corner " << corner << " : " << "n = " << normal << ", t = [" << calcTByDevelopingRow(corner, 1, false) << ", " << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) );
- LOG(6, "] - stable : " << calcStableT(corner) );
-
- //? 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])
- if(!_validTP[corner])
- return ( calcTByDevelopingRow(corner, 1, false) * normal ) >= 0.0;
- else
- return ( calcStableT(corner) * normal ) >= 0.0;
+ // 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
}
-
#endif
};
/// The machine epsilon, used in precision corrections
- const long double TransformedTriangle::MACH_EPS = std::numeric_limits<double>::epsilon();
+ const double TransformedTriangle::MACH_EPS = std::numeric_limits<double>::epsilon();
/// 4.0 * the machine epsilon, represents the precision of multiplication when performing corrections corrections ( f in Grandy )
- const long double TransformedTriangle::MULT_PREC_F = 4.0 * TransformedTriangle::MACH_EPS;
+ const double TransformedTriangle::MULT_PREC_F = 4.0 * TransformedTriangle::MACH_EPS;
/// Threshold for resetting double and triple products to zero; ( F / f in Grandy )
- const long double TransformedTriangle::THRESHOLD_F = 500.0;
+ const double TransformedTriangle::THRESHOLD_F = 20.0;
/// Threshold for what is considered a small enough angle to warrant correction of triple products by Grandy, [57]
const double TransformedTriangle::TRIPLE_PRODUCT_ANGLE_THRESHOLD = 0.1;
- // after transformation to the U-space, coordinates are inaccurate
- // small variations around zero should not be taken into account
- void TransformedTriangle::resetNearZeroCoordinates()
+ // Handle cases where one of the segment (or all) is (almost) in XYZ plane.
+ // We follow Grandy's suggestion and perturb slightly to have exactly h=0 for the segment (Grandy p.447)
+ // Note that if PQR is == to the upper facet of the unit tetra (XYZ), the tetra-corner-inclusion test should take it in,
+ // thanks to Grandy [21] and the fact that S_x test is "<=0" (not <0)
+ // After that, we also snap P,Q,R to the corners if they're very close.
+ void TransformedTriangle::handleDegenerateCases()
{
- for (int i=0; i<15; i++)
- if (fabs(_coords[i])<TransformedTriangle::MACH_EPS*40.0) _coords[i]=0.0;
+ static const TriCorner PT_SEG_MAP[] = {
+ P, Q,
+ Q, R,
+ R, P
+ };
+
+ const double eps = THRESHOLD_F*TransformedTriangle::MULT_PREC_F;
+ for (TriSegment seg = PQ; seg <= RP; seg = TriSegment(seg+1))
+ {
+ // Is h coordinate for both end of segment small enough?
+ int pt1 = PT_SEG_MAP[2*seg], pt2 = PT_SEG_MAP[2*seg+1];
+ if (fabs(_coords[5*pt1+3]) < eps && fabs(_coords[5*pt2+3]) < eps)
+ {
+ // If so, perturb x,y and z to reset h to exactly zero.
+ for (auto pt: {pt1, pt2}) // thx C++17
+ {
+ const double correc = _coords[pt*5+3]/3.; // this should be really small!
+ _coords[pt*5+0] += correc;
+ _coords[pt*5+1] += correc;
+ _coords[pt*5+2] += correc;
+ // And then, if x,y or z very close to 0 or 1, snap exactly to tetra corner:
+ for(int d=0; d < 3; d++)
+ {
+ if (fabs(_coords[5*pt+d]) < eps) _coords[5*pt+d] = 0.0;
+ if (fabs(_coords[5*pt+d]-1) < eps) _coords[5*pt+d] = 1.0;
+ }
+ _coords[pt*5+3] = 0.0;
+ }
+ }
+ }
}
// ----------------------------------------------------------------------------------
* and it is thus the "stable" double products that are stored.
*
*/
- void TransformedTriangle::preCalculateDoubleProducts(void)
+ void TransformedTriangle::preCalculateDoubleProducts()
{
if(_is_double_products_calculated)
return;
for(TriSegment seg = PQ ; seg <= RP ; seg = TriSegment(seg + 1))
{
for(DoubleProduct dp = C_YZ ; dp <= C_10 ; dp = DoubleProduct(dp + 1))
- _doubleProducts[8*seg + dp] = calcUnstableC(seg, dp);
+ {
+ const int idx = 8*seg + dp;
+ _doubleProducts[idx] = calcUnstableC(seg, dp, _deltas[idx]);
+ }
}
std::map<double, TetraCorner> distances;
distances.insert( std::make_pair( dist, corner ) );
}
- // first element -> minimum distance
+ // first element -> minimum distance (because map is sorted)
const TetraCorner minCorner = distances.begin()->second;
resetDoubleProducts(seg, minCorner);
distances.clear();
}
}
-
+
// -- (2) check that each double product satisfies Grandy, [47], else set to 0
for(TriSegment seg = PQ ; seg <= RP ; seg = TriSegment(seg + 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
- 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 int idx = 8*seg+dp;
- const double term1 = _coords[5*pt1 + off1] * _coords[5*pt2 + off2];
- const double term2 = _coords[5*pt1 + off2] * _coords[5*pt2 + off1];
-
- const long double delta = MULT_PREC_F * ( std::fabs(term1) + std::fabs(term2) );
-
- if( epsilonEqual(_doubleProducts[8*seg + dp], 0.0, (double)(THRESHOLD_F * delta)))
+ if( epsilonEqual(_doubleProducts[idx], 0.0, THRESHOLD_F * MULT_PREC_F * _deltas[idx]))
{
// debug output
#if LOG_LEVEL >= 5
}
#endif
-
- _doubleProducts[8*seg + dp] = 0.0;
-
+ _doubleProducts[idx] = 0.0;
}
}
}
-
+
_is_double_products_calculated = true;
}
*/
bool TransformedTriangle::areDoubleProductsConsistent(const TriSegment seg) const
{
- const double term1 = _doubleProducts[8*seg + C_YZ] * _doubleProducts[8*seg + C_XH];
- const double term2 = _doubleProducts[8*seg + C_ZX] * _doubleProducts[8*seg + C_YH];
- const double term3 = _doubleProducts[8*seg + C_XY] * _doubleProducts[8*seg + C_ZH];
+ // Careful! Here doubleProducts have not yet been corrected for roundoff errors!
+ // So we need to epsilon-adjust to correctly identify zeros:
+ static const DoubleProduct DP_LST[6] = {C_YZ, C_XH,
+ C_ZX, C_YH,
+ C_XY, C_ZH};
+ double dps[6];
+ for (int i = 0; i < 6; i++)
+ {
+ const double dp = _doubleProducts[8*seg + DP_LST[i]];
+ dps[i] = dp;
+ }
+
+ const double term1 = dps[0] * dps[1];
+ const double term2 = dps[2] * dps[3];
+ const double term3 = dps[4] * dps[5];
LOG(2, "for seg " << seg << " consistency " << term1 + term2 + term3 );
LOG(2, "term1 :" << term1 << " term2 :" << term2 << " term3: " << term3 );
+ // Test for "== 0.0" here is OK since doubleProduct has been fixed for rounding to zero already.
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);
const int num_pos = (term1 > 0.0 ? 1 : 0) + (term2 > 0.0 ? 1 : 0) + (term3 > 0.0 ? 1 : 0);
assert( num_zero + num_neg + num_pos == 3 );
- // calculated geometry is inconsistent if we have one of the following cases
+ // Calculated geometry is inconsistent if we have one of the following cases
// * one term zero and the other two of the same sign
// * two terms zero
// * all terms positive
// * all terms negative
- if(((num_zero == 1 && num_neg != 1) || num_zero == 2 || (num_neg == 0 && num_zero != 3) || num_neg == 3 ))
- {
+ const bool inconsist = (num_zero == 1 && num_neg != 1) ||
+ num_zero == 2 ||
+ (num_neg == 0 && num_zero != 3) ||
+ num_neg == 3;
+ if(inconsist) {
LOG(4, "inconsistent dp found" );
}
- return !((num_zero == 1 && num_neg != 1) || num_zero == 2 || (num_neg == 0 && num_zero != 3) || num_neg == 3 );
-
+ return !inconsist;
}
/**
* the problem of errors due to cancellation.
*
*/
- void TransformedTriangle::preCalculateTripleProducts(void)
+ void TransformedTriangle::preCalculateTripleProducts()
{
if(_is_triple_products_calculated)
- {
- return;
- }
+ return;
// find edge / row to use -> that whose edge makes the smallest angle to the triangle
// use a map to find the minimum
// get edge by using correspondence between Double Product and Edge
TetraEdge edge = TetraEdge(dp);
-
+
// use edge only if it is surrounded by the surface
if( _triangleSurroundsEdgeCache[edge] )
{
const int minRow = anglesForRows.begin()->second;
if(minAngle < TRIPLE_PRODUCT_ANGLE_THRESHOLD)
- {
- _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, true);
- }
+ _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, true);
else
- {
- _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, false);
- }
+ _tripleProducts[corner] = calcTByDevelopingRow(corner, minRow, false);
+
_validTP[corner] = true;
}
else
LOG(6, "Triple product not calculated for corner " << corner );
_tripleProducts[corner] = -3.14159265;
_validTP[corner] = false;
-
}
anglesForRows.clear();
-
}
-
_is_triple_products_calculated = true;
}
}
/**
- * Calculates triple product associated with the given corner of tetrahedron, developing
+ * 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.
*
+ * Consistency with the double product computation and potential cancellation is also done here.
+ *
+ *
* @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
const double cPQ = calcStableC(PQ, dp);
double alpha = 0.0;
-
+
// coordinate to use for projection (Grandy, [57]) with edges
// OX, OY, OZ, XY, YZ, ZX in order :
// (y, z, x, h, h, h)
// for the first three we could also use {2, 0, 1}
static const int PROJECTION_COORDS[6] = { 1, 2, 0, 3, 3, 3 } ;
-
+
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[2] * cPQ };
-
+
const double sumDPProd = coordDPProd[0] + coordDPProd[1] + coordDPProd[2];
const double sumDPProdSq = dot(coordDPProd, coordDPProd);
const double cRPbar = cRP * (1.0 - alpha * coordValues[1] * cRP);
const double cPQbar = cPQ * (1.0 - alpha * coordValues[2] * cPQ);
+ // [ABN] Triple product cancellation logic:
+ // This part is not well described in Grandy (end of p.446) :
+ // "We use a method analogous to (47) to remove imprecise triple products,..."
+ //
+ // Our algo for cancelling a triple product:
+ // - retrieve the deltas associated with each DP involved (because a DP itself is a sum of two terms - see [42]
+ // - multiply them by the coordinate coming from the determinant expansion
+ // - and finally sum the 3 corresponding terms of the developement
+ //
+ // Using directly the DP (as was done before here) leads to issues, since this DP might have been cancelled
+ // already earlier on, and we lost the delta information -> doing this, we risk not cancelling the triple prod
+ // when we should have.
+ const double cQRbar_delta = _deltas[8*QR + dp],
+ cRPbar_delta = _deltas[8*RP + dp],
+ cPQbar_delta = _deltas[8*PQ + dp];
+
// check sign of barred products - should not change
// assert(cQRbar * cQR >= 0.0);
//assert(cRPbar * cRP >= 0.0);
const double q_term = _coords[5*Q + offset] * cRPbar;
const double r_term = _coords[5*R + offset] * cPQbar;
- // 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::fabs(p_term) + std::fabs(q_term) + std::fabs(r_term));
-#endif
+ const double p_delta = std::fabs(_coords[5*P + offset] * cQRbar_delta),
+ q_delta = std::fabs(_coords[5*Q + offset] * cRPbar_delta),
+ r_delta = std::fabs(_coords[5*R + offset] * cPQbar_delta);
+
+ const double delta = p_delta + q_delta + r_delta;
- if( epsilonEqual( p_term + q_term + r_term, 0.0, (double)(THRESHOLD_F * delta)) )
+ if( epsilonEqual( p_term + q_term + r_term, 0.0, THRESHOLD_F * MULT_PREC_F * delta) )
{
LOG(4, "Reset imprecise triple product for corner " << corner << " to zero" );
return 0.0;
#include "TransformedTriangleIntersectTest.hxx"
#include <iostream>
+#include <iomanip>
#include "Log.hxx"
CPPUNIT_ASSERT_EQUAL(false, tri->testSurfaceEdgeIntersection(TT::XY));
// surface-ray (3 possibilities)
- CPPUNIT_ASSERT_EQUAL(true , tri->testSurfaceRayIntersection(TT::X));
+ CPPUNIT_ASSERT_EQUAL(true, tri->testSurfaceRayIntersection(TT::X));
CPPUNIT_ASSERT_EQUAL(false, tri->testSurfaceRayIntersection(TT::Y));
CPPUNIT_ASSERT_EQUAL(false, tri->testSurfaceRayIntersection(TT::Z));
delete tri;
}
+ // Extract from BoxHexa1.med (or was it BoxHexa2.med ?)
+ void TransformedTriangleIntersectTest::testTriangle_vol1()
+ {
+ LOG(1, "+++++++ Testing testTriangle_vol1" );
-} // NAMESPACE
+ double coords[9] = {
+ -0.8088235294117645, 0, 0.55882352941176472,
+ 0.44117647058823506, 0, 0.55882352941176483,
+ -0.89215686274509864, 1.3333333333333339, 0.55882352941176483};
+
+ double refVol = 0.054383777732546296;
+
+ TransformedTriangle tri(&coords[0], &coords[3], &coords[6]);
+ const double vol = tri.calculateIntersectionVolume();
+
+ LOG(3, " --- Final points in polygon A");
+ for(const auto& pt: tri._polygonA)
+ LOG(3,vToStr(pt));
+ LOG(3, " --- Final points in polygon B");
+ for(const auto& pt: tri._polygonB)
+ LOG(3,vToStr(pt));
+
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10);
+ }
+
+ // Extract from BoxHexa1.med (or was it BoxHexa2.med ?)
+ void TransformedTriangleIntersectTest::testTriangle_vol2()
+ {
+ LOG(1, "+++++++ Testing testTriangle_vol2" );
+
+ double coords[9] = {
+ 0.44117647058823506, 0, 0.55882352941176483,
+ -0.55882352941176472, 0, 1.5588235294117649,
+ -0.89215686274509864, 1.3333333333333339, 0.55882352941176483 };
+
+ double refVol = -0.06869529818848;
+
+ TransformedTriangle tri(&coords[0], &coords[3], &coords[6]);
+ const double vol = tri.calculateIntersectionVolume();
+
+ LOG(3, " --- Final points in polygon A");
+ for(const auto& pt: tri._polygonA)
+ LOG(3,vToStr(pt));
+ LOG(3, " --- Final points in polygon B");
+ for(const auto& pt: tri._polygonB)
+ LOG(3,vToStr(pt));
+
+ CPPUNIT_ASSERT(tri.isTriangleInPlaneOfFacet(TransformedTriangle::XYZ));
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10);
+ }
+
+ // Extract from BoxModSmall1.med
+ void TransformedTriangleIntersectTest::testTriangle_vol3()
+ {
+ LOG(1, "+++++++ Testing testTriangle_vol3" );
+
+ double coords[9] = {
+ 0, -4.4408920985006262e-16, 1,
+ -1.2062474433365091, -0.037350951323461778, 2.1879983126221099,
+ 0.49877186496532655, 0.59827776894780405, 0.79353793765518521
+ };
+ double refVol = -0.051135429735185;
+
+ TransformedTriangle tri(&coords[0], &coords[3], &coords[6]);
+ const double vol = tri.calculateIntersectionVolume();
+
+ LOG(3, " --- Final points in polygon A");
+ for(const auto& pt: tri._polygonA)
+ LOG(3,vToStr(pt));
+ LOG(3, " --- Final points in polygon B");
+ for(const auto& pt: tri._polygonB)
+ LOG(3,vToStr(pt));
+
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10);
+ }
+
+ // Taken from Box1Moderate.med
+ void TransformedTriangleIntersectTest::testTriangle_vol4()
+ {
+ LOG(1, "+++++++ Testing testTriangle_vol4" );
+
+ double coords[9] = {
+ 1, 3.552713678800501e-15, 0,
+ 2.022774182629973, -1.020222639063029, -0.01375178680446254,
+ 0.7495960843059706, 0.1125313911637846, 0.7430770879625861
+ };
+ double refVol = -0.00060846166394417;
+
+ TransformedTriangle tri(&coords[0], &coords[3], &coords[6]);
+ const double vol = tri.calculateIntersectionVolume();
+
+ LOG(3, " --- Final points in polygon A");
+ for(const auto& pt: tri._polygonA)
+ LOG(3,vToStr(pt));
+ LOG(3, " --- Final points in polygon B");
+ for(const auto& pt: tri._polygonB)
+ LOG(3,vToStr(pt));
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10);
+ }
+
+ // Taken from Box1Moderate.med too
+ void TransformedTriangleIntersectTest::testTriangle_vol5()
+ {
+ LOG(1, "+++++++ Testing testTriangle_vol5" );
+
+ double coords[9] = {
+ 1, 3.552713678800501e-15, 0,
+ -3.552713678800501e-15, 0, 0.9999999999999982,
+ 0, 1.000000000000004, -8.881784197001252e-16
+ };
+ double refVol = -1/6.;
+
+ TransformedTriangle tri(&coords[0], &coords[3], &coords[6]);
+ const double vol = tri.calculateIntersectionVolume();
+ LOG(3, " --- Final points in polygon A");
+ for(const auto& pt: tri._polygonA)
+ LOG(3,vToStr(pt));
+ LOG(3, " --- Final points in polygon B");
+ for(const auto& pt: tri._polygonB)
+ LOG(3,vToStr(pt));
+
+ CPPUNIT_ASSERT(tri.isTriangleInPlaneOfFacet(TransformedTriangle::XYZ));
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10);
+ }
+
+ // Taken from Box1Moderate.med again
+ void TransformedTriangleIntersectTest::testTriangle_vol6()
+ {
+ LOG(1, "+++++++ Testing testTriangle_vol6" );
+
+ double coords[9] = {
+ 1.000000000000004, 0, 0,
+ 0, 0, 0.9999999999999929,
+ 3.552713678800501e-15, 1, 0};
+ double refVol = -1/6.;
+
+ TransformedTriangle tri(&coords[0], &coords[3], &coords[6]);
+ const double vol = tri.calculateIntersectionVolume();
+
+ LOG(3, " --- Final points in polygon A");
+ for(const auto& pt: tri._polygonA)
+ LOG(3,vToStr(pt));
+ LOG(3, " --- Final points in polygon B");
+ for(const auto& pt: tri._polygonB)
+ LOG(3,vToStr(pt));
+
+ CPPUNIT_ASSERT(tri.isTriangleInPlaneOfFacet(TransformedTriangle::XYZ));
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10);
+ }
+
+ // Taken from MEDCouplingRemapperTest.py::testCellToNodeReverse3D()
+ void TransformedTriangleIntersectTest::testTriangle_vol7()
+ {
+ LOG(1, "+++++++ Testing testTriangle_vol7" );
+
+ double coords[9] = {
+ 0.9999999999999999, 3.700743415417188e-17, 3.700743415417188e-17,
+ -5.551115123125783e-17, 0, 1,
+ 3.700743415417188e-17, 0.9999999999999999, 3.700743415417188e-17
+ };
+
+ double refVol = -1/6.;
+
+ TransformedTriangle tri(&coords[0], &coords[3], &coords[6]);
+ const double vol = tri.calculateIntersectionVolume();
+
+ LOG(3, " --- Final points in polygon A");
+ for(const auto& pt: tri._polygonA)
+ LOG(3,vToStr(pt));
+ LOG(3, " --- Final points in polygon B");
+ for(const auto& pt: tri._polygonB)
+ LOG(3,vToStr(pt));
+
+ CPPUNIT_ASSERT(tri.isTriangleInPlaneOfFacet(TransformedTriangle::XYZ));
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10);
+
+ // Polygon A should not contain anything else than tetra corners:
+ double barycenter[3];
+ tri.calculatePolygonBarycenter(TransformedTriangle::A, barycenter);
+ tri.sortIntersectionPolygon(TransformedTriangle::A, barycenter);
+ const double eps = TransformedTriangle::THRESHOLD_F * TransformedTriangle::MULT_PREC_F;
+ for (auto& p: tri._polygonA)
+ for (int d=0; d<3; d++)
+ {
+ CPPUNIT_ASSERT(fabs(p[d]) < eps || fabs(p[d]-1) < eps);
+ }
+ }
+
+ // This test is useful if TransformedTriangle::THRESHOLD_F = 20 (like in Grandy's article)
+ // With this value, it is an almost degenerate case : the triangle is very close to XYZ plane, but this still
+ // falls in the 'non-degenerate' category. Total volume is 1/6 but coming from
+ // contributions from A and B.
+ // With current value of threshold (=100, chosen because of P1P1 intersector), triangle becomes completely degenerated.
+ void TransformedTriangleIntersectTest::testTriangle_vol8()
+ {
+ LOG(1, "+++++++ Testing testTriangle_vol8" );
+
+ double coords[9] = {
+ 0.9999999999999787, 7.105427357601002e-15, -3.552713678800501e-15,
+ -1.4210854715202e-14, 0, 0.9999999999999964,
+ -7.105427357601002e-15, 1.000000000000014, -3.552713678800501e-15
+ };
+
+ double refVol = -1/6.;
+
+ TransformedTriangle tri(&coords[0], &coords[3], &coords[6]);
+ const double vol = tri.calculateIntersectionVolume();
+
+ LOG(3, " --- Final points in polygon A");
+ for(const auto& pt: tri._polygonA)
+ LOG(3,vToStr(pt));
+ LOG(3, " --- Final points in polygon B");
+ for(const auto& pt: tri._polygonB)
+ LOG(3,vToStr(pt));
+
+// CPPUNIT_ASSERT (! tri.isTriangleInPlaneOfFacet(TransformedTriangle::XYZ) ); // not degenerate if threshold == 20.0
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vol, refVol, 1.0e-10);
+ }
+
+
+} // NAMESPACE
CPPUNIT_TEST( testTriangle12 );
CPPUNIT_TEST( testTriangle13 );
+
+ // Tests for degenerated cases where PQR is almost in XYZ plane:
+ CPPUNIT_TEST( testTriangle_vol1 );
+ CPPUNIT_TEST( testTriangle_vol2 );
+ CPPUNIT_TEST( testTriangle_vol3 );
+ CPPUNIT_TEST( testTriangle_vol4 );
+ CPPUNIT_TEST( testTriangle_vol5 );
+ CPPUNIT_TEST( testTriangle_vol6 );
+ CPPUNIT_TEST( testTriangle_vol7 );
+ CPPUNIT_TEST( testTriangle_vol8 );
+
CPPUNIT_TEST_SUITE_END();
typedef INTERP_KERNEL::TransformedTriangle::TriSegment TriSegment;
public:
void testTriangle1();
-
void testTriangle2();
-
void testTriangle3();
-
void testTriangle4();
-
void testTriangle5();
-
void testTriangle6();
-
void testTriangle7();
-
void testTriangle8();
-
void testTriangle9();
-
void testTriangle10();
-
void testTriangle11();
-
void testTriangle12();
-
void testTriangle13();
- private:
-
+ void testTriangle_vol1();
+ void testTriangle_vol2();
+ void testTriangle_vol3();
+ void testTriangle_vol4();
+ void testTriangle_vol5();
+ void testTriangle_vol6();
+ void testTriangle_vol7();
+ void testTriangle_vol8();
+
};
}
-
-
-
-
-
#endif
double c_vals[3 * 8];
for(TriSegment seg = TransformedTriangle::PQ ; seg <= TransformedTriangle::RP ; seg = TriSegment(seg + 1)) {
-
- c_vals[seg*8 + 0] = tri1->calcUnstableC(seg, TransformedTriangle::C_XY);
- c_vals[seg*8 + 1] = tri1->calcUnstableC(seg, TransformedTriangle::C_YZ);
- c_vals[seg*8 + 2] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZX);
- c_vals[seg*8 + 3] = tri1->calcUnstableC(seg, TransformedTriangle::C_XH);
- c_vals[seg*8 + 4] = tri1->calcUnstableC(seg, TransformedTriangle::C_YH);
- c_vals[seg*8 + 5] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZH);
- c_vals[seg*8 + 6] = tri1->calcUnstableC(seg, TransformedTriangle::C_01);
- c_vals[seg*8 + 7] = tri1->calcUnstableC(seg, TransformedTriangle::C_10);
+ double dnu;
+ c_vals[seg*8 + 0] = tri1->calcUnstableC(seg, TransformedTriangle::C_XY, dnu);
+ c_vals[seg*8 + 1] = tri1->calcUnstableC(seg, TransformedTriangle::C_YZ, dnu);
+ c_vals[seg*8 + 2] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZX, dnu);
+ c_vals[seg*8 + 3] = tri1->calcUnstableC(seg, TransformedTriangle::C_XH, dnu);
+ c_vals[seg*8 + 4] = tri1->calcUnstableC(seg, TransformedTriangle::C_YH, dnu);
+ c_vals[seg*8 + 5] = tri1->calcUnstableC(seg, TransformedTriangle::C_ZH, dnu);
+ c_vals[seg*8 + 6] = tri1->calcUnstableC(seg, TransformedTriangle::C_01, dnu);
+ c_vals[seg*8 + 7] = tri1->calcUnstableC(seg, TransformedTriangle::C_10, dnu);
}
// find unstable products to check for consistency (Grandy [46])
for(TriSegment seg = TransformedTriangle::PQ ; seg <= TransformedTriangle::RP ; seg = TriSegment(seg + 1))
{
- const double c_xy = tri2->calcUnstableC(seg, TransformedTriangle::C_XY);
- const double c_yz = tri2->calcUnstableC(seg, TransformedTriangle::C_YZ);
- const double c_zx = tri2->calcUnstableC(seg, TransformedTriangle::C_ZX);
- const double c_xh = tri2->calcUnstableC(seg, TransformedTriangle::C_XH);
- const double c_yh = tri2->calcUnstableC(seg, TransformedTriangle::C_YH);
- const double c_zh = tri2->calcUnstableC(seg, TransformedTriangle::C_ZH);
+ double dnu;
+ const double c_xy = tri2->calcUnstableC(seg, TransformedTriangle::C_XY, dnu);
+ const double c_yz = tri2->calcUnstableC(seg, TransformedTriangle::C_YZ, dnu);
+ const double c_zx = tri2->calcUnstableC(seg, TransformedTriangle::C_ZX, dnu);
+ const double c_xh = tri2->calcUnstableC(seg, TransformedTriangle::C_XH, dnu);
+ const double c_yh = tri2->calcUnstableC(seg, TransformedTriangle::C_YH, dnu);
+ const double c_zh = tri2->calcUnstableC(seg, TransformedTriangle::C_ZH, dnu);
const int num_zero = (c_yz*c_xh == 0.0 ? 1 : 0) + (c_zx*c_yh == 0.0 ? 1 : 0) + (c_xy*c_zh == 0.0 ? 1 : 0);
const int num_neg = (c_yz*c_xh < 0.0 ? 1 : 0) + (c_zx*c_yh < 0.0 ? 1 : 0) + (c_xy*c_zh < 0.0 ? 1 : 0);
CPPUNIT_ASSERT_DOUBLES_EQUAL(6944.4444444444443,volume,1e-9);
}
+ /**
+ * Extract from a single tetra from BoxHexa1.med and BoxHexa2.med.
+ * Symetry test was failing.
+ */
+ void UnitTetraIntersectionBaryTest::test_UnitTetraIntersectionBary_14()
+ {
+ double S[] = {
+ 25.0, 96.0, 0.0,
+ 25.0, 120.0, 0.0,
+ 37.5, 120.0, 0.0,
+ 25.0, 120.0, 11.333333333333339255};
+
+ double T[] = {
+ 25.0, 90.0, 6.333333333333339255,
+ 25.0, 120.0, 6.3333333333333357018,
+ 41.6666666666666714036, 120.0, 6.3333333333333348136,
+ 25.0, 120.0, 17.6666666666666714036};
+
+ mcIdType conn[4] = { 0,1,2,3 };
+ const double* tnodes[4]={ T, T+3, T+6, T+9 };
+ const double* snodes[4]={ S, S+3, S+6, S+9 };
+ const double refVol = 48.6591695501729;
+
+ __MESH_DUMMY dummyMesh;
+ SplitterTetra<__MESH_DUMMY> src( dummyMesh, snodes, conn );
+ double volume = src.intersectTetra( tnodes );
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(refVol,volume,1e-9);
+
+ // Now the other way round:
+ SplitterTetra<__MESH_DUMMY> tgt( dummyMesh, tnodes, conn );
+ double volume2 = tgt.intersectTetra( snodes );
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(refVol,volume2,1e-9);
+ }
+
+
void UnitTetraIntersectionBaryTest::test_TetraAffineTransform_reverseApply()
{
double nodes[12] = { -4.0, 9.0, 3.0,
class INTERPKERNELTEST_EXPORT UnitTetraIntersectionBaryTest : public CppUnit::TestFixture
{
CPPUNIT_TEST_SUITE( UnitTetraIntersectionBaryTest );
+ CPPUNIT_TEST( test_UnitTetraIntersectionBary_14 );
CPPUNIT_TEST( test_UnitTetraIntersectionBary_13 );
CPPUNIT_TEST( test_UnitTetraIntersectionBary_12 );
CPPUNIT_TEST( test_UnitTetraIntersectionBary_1 );
void test_UnitTetraIntersectionBary_11();
void test_UnitTetraIntersectionBary_12();
void test_UnitTetraIntersectionBary_13();
+ void test_UnitTetraIntersectionBary_14();
void test_TetraAffineTransform_reverseApply();
void test_barycentric_coords();
void test_cuboid_mapped_coords_3D();
