src/INTERP_KERNEL/TransformedTriangleIntersect.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "TransformedTriangle.hxx"
  21 #include <iostream>
  22 #include <fstream>
  23 #include <cassert>
  24 #include <cmath>
  25 #include "VectorUtils.hxx"
  26
  27 namespace INTERP_KERNEL
  28 {
  29
  30   // ----------------------------------------------------------------------------------
  31   //  Correspondence tables describing all the variations of formulas.
  32   // ----------------------------------------------------------------------------------
  33
  34   /// \brief Correspondence between facets and double products.
  35   ///
  36   /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
  37   const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_SEG_FACET_INTERSECTION[12] =
  38     {
  39       C_XH, C_XY, C_ZX, // OYZ
  40       C_YH, C_YZ, C_XY, // OZX
  41       C_ZH, C_ZX, C_YZ, // OXY
  42       C_XH, C_YH, C_ZH  // XYZ
  43     };
  44
  45   /// \brief Signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION.
  46   ///
  47   /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
  48   const double TransformedTriangle::SIGN_FOR_SEG_FACET_INTERSECTION[12] =
  49     {
  50       1.0, 1.0, -1.0,
  51       1.0, 1.0, -1.0,
  52       1.0, 1.0, -1.0,
  53       1.0, 1.0,  1.0
  54     };
  55
  56   /// \brief Coordinates of corners of tetrahedron.
  57   ///
  58   /// Use 3*Corner + coordinate as index
  59   const double TransformedTriangle::COORDS_TET_CORNER[12] =
  60     {
  61       0.0, 0.0, 0.0,
  62       1.0, 0.0, 0.0,
  63       0.0, 1.0, 0.0,
  64       0.0, 0.0, 1.0
  65     };
  66
  67   /// \brief Indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION
  68   /// for the calculation of the coordinates (x,y,z) of the intersection points
  69   /// for Segment-Facet and Segment-Edge intersections.
  70   ///
  71   /// Use 3*facet + coordinate as index. -1 indicates that the coordinate is 0.
  72   const int TransformedTriangle::DP_INDEX[12] =
  73     {
  74       // x, y, z
  75       -1, 1, 2,  // OYZ
  76       5, -1, 4,  // OZX
  77       7, 8, -1,  // OXY
  78       9, 10, 11  // XYZ
  79     };
  80
  81   /// \brief Correspondence edge - corners.
  82   ///
  83   /// Gives the two corners associated with each edge
  84   /// Use 2*edge + {0, 1} as index
  85   const TransformedTriangle::TetraCorner TransformedTriangle::CORNERS_FOR_EDGE[12] =
  86     {
  87       O, X, // OX
  88       O, Y, // OY
  89       O, Z, // OZ
  90       X, Y, // XY
  91       Y, Z, // YZ
  92       Z, X  // ZX
  93     };
  94
  95   /// \brief Correspondence edge - facets.
  96   ///
  97   /// Gives the two facets shared by and edge. Use 2*facet + {0, 1} as index
  98   const TransformedTriangle::TetraFacet TransformedTriangle::FACET_FOR_EDGE[12] =
  99     {
 100       OXY, OZX, // OX
 101       OXY, OYZ, // OY
 102       OZX, OYZ, // OZ
 103       OXY, XYZ, // XY
 104       OYZ, XYZ, // YZ
 105       OZX, XYZ  // ZX
 106     };
 107
 108   /// \brief Correspondence corners - edges.
 109   ///
 110   /// Gives edges meeting at a given corner. Use 3*corner + {0,1,2} as index
 111   const TransformedTriangle::TetraEdge TransformedTriangle::EDGES_FOR_CORNER[12] =
 112     {
 113       OX, OY, OZ, // O
 114       OX, XY, ZX, // X
 115       OY, XY, YZ, // Y
 116       OZ, ZX, YZ  // Z
 117     };
 118
 119   /// \brief Double products to use in halfstrip intersection tests.
 120   ///
 121   /// Use 4*(offset_edge) + {0,1,2,3} as index. offset_edge = edge - 3  (so that XY -> 0, YZ -> 1, ZX -> 2)
 122   /// Entries with offset 0 and 1 are for the first condition (positive product)
 123   /// and those with offset 2 and 3 are for the second condition (negative product).
 124   const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_HALFSTRIP_INTERSECTION[12] =
 125     {
 126       C_10, C_01, C_ZH, C_10, // XY
 127       C_01, C_XY, C_XH, C_01, // YZ
 128       C_XY, C_10, C_YH, C_XY  // ZX
 129     };
 130
 131   /// \brief Double products to use in segment-ray test.
 132   ///
 133   /// Use 7*corner_offset + {0,1,2,3,4,5,6} as index. corner_offset = corner - 1 (so that X -> 0, Y-> 1, Z->2)
 134   /// Entries with offset 0 are for first condition (zero double product) and the rest are for condition 3 (in the same
 135   /// order as in the article)
 136   const TransformedTriangle::DoubleProduct TransformedTriangle::DP_SEGMENT_RAY_INTERSECTION[21] =
 137     {
 138       C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
 139       C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_10, // Y
 140       C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_01  // Z
 141     };
 142
 143   /**
 144    * Calculates the point of intersection between the given edge of the tetrahedron and the
 145    * triangle PQR. (Grandy, eq [22])
 146    *
 147    * @pre   testSurfaceEdgeIntersection(edge) returns true
 148    * @param edge   edge of tetrahedron
 149    * @param pt     array of three doubles in which to store the coordinates of the intersection point
 150    */
 151   void TransformedTriangle::calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const
 152   {
 153     assert(edge < H01);
 154
 155     // barycentric interpolation between points A and B
 156     // : (x,y,z)* = (1-alpha)*A + alpha*B where
 157     // alpha = t_A / (t_A - t_B)
 158
 159     const TetraCorner corners[2] =
 160       {
 161         CORNERS_FOR_EDGE[2*edge],
 162         CORNERS_FOR_EDGE[2*edge + 1]
 163       };
 164
 165     // calculate alpha
 166     const double tA = calcStableT(corners[0]);
 167     const double tB = calcStableT(corners[1]);
 168     const double alpha = tA / (tA - tB);
 169
 170     // calculate point
 171     LOG(4, "corner A = " << strTC(corners[0]) << " corner B = " << strTC(corners[1]) );
 172     LOG(4, "tA = " << tA << " tB = " << tB << " alpha= " << alpha );
 173     for(int i = 0; i < 3; ++i)
 174       {
 175
 176         pt[i] = (1 - alpha)*COORDS_TET_CORNER[3*corners[0] + i] + alpha*COORDS_TET_CORNER[3*corners[1] + i];
 177 #if 0
 178         pt[i] = (1 - alpha) * getCoordinateForTetCorner<corners[0], i>() +
 179           alpha * getCoordinateForTetCorner<corners[0], i>();
 180 #endif
 181         LOG(6, pt[i] );
 182         assert(pt[i] >= 0.0);
 183         assert(pt[i] <= 1.0);
 184       }
 185   }
 186
 187   /**
 188    * Calculates the point of intersection between the given segment of the triangle
 189    * and the given facet of the tetrahedron. (Grandy, eq. [23])
 190    *
 191    * @pre   testSurfaceEdgeIntersection(seg, facet) returns true
 192    *
 193    * @param seg    segment of the triangle
 194    * @param facet  facet of the tetrahedron
 195    * @param pt     array of three doubles in which to store the coordinates of the intersection point
 196    */
 197   void TransformedTriangle::calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const
 198   {
 199     // calculate s
 200     double s = 0.0;
 201     for(int i = 0; i < 3; ++i)
 202       {
 203         const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[3*facet + i];
 204         const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + i];
 205         s -= sign * calcStableC(seg, dp);
 206       }
 207
 208     assert(s != 0.0);
 209
 210     // calculate coordinates of intersection point
 211     for(int i = 0 ; i < 3; ++i)
 212       {
 213         const int dpIdx = DP_INDEX[3*facet + i];
 214
 215         if(dpIdx < 0)
 216           {
 217             pt[i] = 0.0;
 218           }
 219         else
 220           {
 221             const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx];
 222             const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
 223             pt[i] = -( sign * calcStableC(seg, dp) ) / s;
 224
 225             LOG(4, "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i]  );
 226             LOG(4, "c(" << seg << ", " << dp << ") = " <<  sign * calcStableC(seg, dp) );
 227             assert(pt[i] >= 0.0);
 228             assert(pt[i] <= 1.0);
 229           }
 230       }
 231
 232   }
 233
 234   /**
 235    * Tests if the given segment of the triangle intersects the given edge of the tetrahedron (Grandy, eq. [20]
 236    *
 237    * @param seg    segment of the triangle
 238    * @param edge   edge of tetrahedron
 239    * @return      true if the segment intersects the edge
 240    */
 241   bool TransformedTriangle::testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const
 242   {
 243       {
 244         // check condition that the double products for one of the two
 245         // facets adjacent to the edge has a positive product
 246         bool isFacetCondVerified = false;
 247         TetraFacet facet[2];
 248         for(int i = 0 ; i < 2 ; ++i)
 249           {
 250             facet[i] = FACET_FOR_EDGE[2*edge + i];
 251
 252             // find the two c-values -> the two for the other edges of the facet
 253             int idx1 = 0 ;
 254             int idx2 = 1;
 255             DoubleProduct dp1 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1];
 256             DoubleProduct dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2];
 257
 258             if(dp1 == DoubleProduct( edge ))
 259               {
 260                 idx1 = 2;
 261                 dp1 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1];
 262               }
 263             else if(dp2 == DoubleProduct( edge ))
 264               {
 265                 idx2 = 2;
 266                 dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2];
 267               }
 268
 269             const double c1 = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1]*calcStableC(seg, dp1);
 270             const double c2 = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2]*calcStableC(seg, dp2);
 271
 272             //isFacetCondVerified = isFacetCondVerified || c1*c2 > 0.0;
 273             if(c1*c2 > 0.0)
 274               {
 275                 isFacetCondVerified = true;
 276               }
 277           }
 278
 279         if(!isFacetCondVerified)
 280           {
 281             return false;
 282           }
 283         else
 284           {
 285             return testSegmentIntersectsFacet(seg, facet[0]) || testSegmentIntersectsFacet(seg, facet[1]);
 286           }
 287       }
 288   }
 289
 290   /**
 291    * Calculates the point of intersection between the given segment of the triangle
 292    * and the given edge of the tetrahedron. (Grandy, eq. [25])
 293    *
 294    * @pre   testSegmentEdgeIntersection(seg, edge) returns true
 295    *
 296    * @param seg    segment of the triangle
 297    * @param edge   edge of the tetrahedron
 298    * @param pt     array of three doubles in which to store the coordinates of the intersection point
 299    */
 300   void TransformedTriangle::calcIntersectionPtSegmentEdge(const TriSegment seg, const TetraEdge edge, double* pt) const
 301   {
 302     assert(edge < H01);
 303
 304     // get the two facets associated with the edge
 305     static const TetraFacet FACETS_FOR_EDGE[12] =
 306       {
 307         OXY, OZX, // OX
 308         OXY, OYZ, // OY
 309         OZX, OYZ, // OZ
 310         OXY, XYZ, // XY
 311         OYZ, XYZ, // YZ
 312         OZX, XYZ  // ZX
 313       };
 314
 315     const TetraFacet facets[2] =
 316       {
 317         FACETS_FOR_EDGE[2*edge],
 318         FACETS_FOR_EDGE[2*edge + 1]
 319       };
 320
 321     // calculate s for the two edges
 322     double s[2];
 323     for(int i = 0; i < 2; ++i)
 324       {
 325         s[i] = 0.0;
 326         for(int j = 0; j < 3; ++j)
 327           {
 328             const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[3*facets[i] + j];
 329             const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[3*facets[i] + j];
 330             s[i] += sign * calcStableC(seg, dp);
 331           }
 332       }
 333
 334     // calculate denominator
 335     const double denominator = s[0]*s[0] + s[1]*s[1];
 336
 337     // calculate intersection point
 338     for(int i = 0; i < 3; ++i)
 339       {
 340         // calculate double product values for the two faces
 341         double c[2];
 342         for(int j = 0 ; j < 2; ++j)
 343           {
 344             const int dpIdx = DP_INDEX[3*facets[j] + i];
 345             if (dpIdx != -1)
 346               {
 347                 const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx];
 348                 const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
 349                 c[j] = sign * calcStableC(seg, dp);
 350               }
 351             else
 352               c[j] = 0.0;
 353           }
 354
 355         pt[i] = (c[0] * s[0] + c[1] * s[1]) / denominator;
 356       }
 357   }
 358
 359
 360   /**
 361    * Tests if the given segment of the triangle intersects the given corner of the tetrahedron.
 362    * (Grandy, eq. [21]).
 363    *
 364    * @param seg    segment of the triangle
 365    * @param corner corner of the tetrahedron
 366    * @return      true if the segment intersects the corner
 367    */
 368   bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const
 369   {
 370     // facets meeting at a given corner
 371     static const TetraFacet FACETS_FOR_CORNER[12] =
 372       {
 373         OXY, OYZ, OZX, // O
 374         OZX, OXY, XYZ, // X
 375         OYZ, XYZ, OXY, // Y
 376         OZX, XYZ, OYZ  // Z
 377       };
 378
 379     // check segment intersect a facet
 380     for(int i = 0 ; i < 3 ; ++i)
 381       {
 382         const TetraFacet facet = FACETS_FOR_CORNER[3*corner + i];
 383         if(testSegmentIntersectsFacet(seg, facet))
 384           return true;
 385       }
 386
 387     return false;
 388   }
 389
 390   /**
 391    * Tests if the given segment of the triangle intersects the half-strip above the
 392    * given edge of the h = 0 plane. (Grandy, eq. [30])
 393    *
 394    * @param seg    segment of the triangle
 395    * @param edge   edge of the h = 0 plane of the tetrahedron (XY, YZ, ZX)
 396    * @return      true if the upwards ray from the corner intersects the triangle.
 397    */
 398   bool TransformedTriangle::testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edge)
 399   {
 400     // get right index here to avoid "filling out" array
 401     const int edgeIndex = static_cast<int>(edge) - 3;
 402
 403     // double products used in test
 404     // products 1 and 2 for each edge -> first condition in Grandy [30]
 405     // products 3 and 4 for each edge -> third condition
 406     // NB : some uncertainty whether these last are correct
 407     // DP_FOR_HALFSTRIP_INTERSECTION
 408
 409     // facets to use in second condition (S_m)
 410     static const TetraFacet FACET_FOR_HALFSTRIP_INTERSECTION[3] =
 411       {
 412         NO_TET_FACET, // XY -> special case : test with plane H = 0
 413         OYZ, // YZ
 414         OZX  // ZX
 415       };
 416
 417     const double cVals[4] =
 418       {
 419         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex]),
 420         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 1]),
 421         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 2]),
 422         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 3])
 423       };
 424
 425     const TetraFacet facet =  FACET_FOR_HALFSTRIP_INTERSECTION[edgeIndex];
 426
 427     // special case : facet H = 0
 428     const bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet);
 429     LOG(4, "Halfstrip tests (" << strTriS(seg) << ", " << strTE(edge) << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) );
 430     LOG(4, "c2 = " << cVals[2] << ", c3 = " << cVals[3] );
 431
 432     return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0);
 433   }
 434
 435   /**
 436    * Calculates the point of intersection between the given segment of the triangle
 437    * and the halfstrip above the given edge of the tetrahedron. (Grandy, eq. [31])
 438    *
 439    * @pre   testSegmentHalfstripIntersection(seg, edge) returns true
 440    *
 441    * @param seg    segment of the triangle
 442    * @param edge   edge of the tetrahedron defining the halfstrip
 443    * @param pt     array of three doubles in which to store the coordinates of the intersection point
 444    */
 445   void TransformedTriangle::calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const
 446   {
 447     assert(edge > OZ);
 448     assert(edge < H01);
 449
 450     // get right index here to avoid "filling out" array
 451     const int edgeIndex = static_cast<int>(edge) - 3;
 452     assert(edgeIndex >= 0);
 453     assert(edgeIndex < 3);
 454
 455     // Barycentric interpolation on the edge
 456     // for edge AB : (x,y,z)* = (1-alpha) * A + alpha * B
 457     // where alpha = cB / (cB - cA)
 458
 459     const double cA = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex]);
 460     const double cB = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 1]);
 461     assert(cA != cB);
 462
 463     const double alpha = cA / (cA - cB);
 464
 465     for(int i = 0; i < 3; ++i)
 466       {
 467         const TetraCorner corners[2] =
 468           {
 469             CORNERS_FOR_EDGE[2*edge],
 470             CORNERS_FOR_EDGE[2*edge + 1]
 471           };
 472
 473         const double cornerCoords[2] =
 474           {
 475             COORDS_TET_CORNER[3*corners[0] + i],
 476             COORDS_TET_CORNER[3*corners[1] + i]
 477           };
 478
 479         pt[i] = (1 - alpha) * cornerCoords[0] + alpha * cornerCoords[1];
 480         LOG(6, pt[i] );
 481         assert(pt[i] >= 0.0);
 482         assert(pt[i] <= 1.0);
 483       }
 484     assert(epsilonEqualRelative(pt[0] + pt[1] + pt[2], 1.0));
 485   }
 486
 487   /**
 488    * Tests if the given segment of triangle PQR intersects the ray pointing
 489    * in the upwards z - direction from the given corner of the tetrahedron. (Grandy eq. [29])
 490    *
 491    * @param seg    segment of the triangle PQR
 492    * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z)
 493    * @return      true if the upwards ray from the corner intersects the segment.
 494    */
 495   bool TransformedTriangle::testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const
 496   {
 497     assert(corner == X || corner == Y || corner == Z);
 498     LOG(4, "Testing seg - ray intersection for seg = " << seg << ", corner = " << corner );
 499
 500     // readjust index since O is not used
 501     const int cornerIdx = static_cast<int>(corner) - 1;
 502
 503     // facets to use
 504     //? not sure this is correct
 505     static const TetraFacet FIRST_FACET_SEGMENT_RAY_INTERSECTION[3] =
 506       {
 507         OZX, // X
 508         OYZ, // Y
 509         OZX, // Z
 510       };
 511
 512     // cond 2
 513     const bool cond21 = testSegmentIntersectsFacet(seg, FIRST_FACET_SEGMENT_RAY_INTERSECTION[cornerIdx]);
 514     const bool cond22  = (corner == Z) ? testSegmentIntersectsFacet(seg, OYZ) : testSegmentIntersectsHPlane(seg);
 515
 516     if(!(cond21 || cond22))
 517       {
 518         LOG(4, "SR fails at cond 2 : cond21 = " << cond21 << ", cond22 = " << cond22 );
 519         return false;
 520       }
 521
 522     // cond 3
 523     const double cVals[6] =
 524       {
 525         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 1]),
 526         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 2]),
 527         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 3]),
 528         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 4]),
 529         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 5]),
 530         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 6]),
 531       };
 532
 533     // cond. 3
 534     if(( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) >= 0.0)
 535       {
 536         LOG(4, "SR fails at cond 3 : " << (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5]  );
 537       }
 538     return ( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) < 0.0;
 539
 540   }
 541
 542   // /////////////////////////////////////////////////////////////////////////////////
 543   //  Utility methods used in intersection tests                       ///////////////
 544   // /////////////////////////////////////////////////////////////////////////////////
 545   /**
 546    * Tests if the triangle PQR surrounds the axis on which the
 547    * given edge of the tetrahedron lies.
 548    *
 549    * @param edge   edge of tetrahedron
 550    * @return      true if PQR surrounds edge, false if not (see Grandy, eq. [53])
 551    */
 552   bool TransformedTriangle::testTriangleSurroundsEdge(const TetraEdge edge) const
 553   {
 554     // NB DoubleProduct enum corresponds to TetraEdge enum according to Grandy, table III
 555     // so we can use the edge directly
 556
 557     const double cPQ = calcStableC(PQ, DoubleProduct(edge));
 558     const double cQR = calcStableC(QR, DoubleProduct(edge));
 559     const double cRP = calcStableC(RP, DoubleProduct(edge));
 560
 561     LOG(5, "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" );
 562
 563     // if two or more c-values are zero we disallow x-edge intersection
 564     // Grandy, p.446
 565     // test for == 0.0 is OK here, if you look at the way double products were pre-comuted:
 566     const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
 567
 568     if(numZeros >= 2 )
 569       {
 570         LOG(5, "TriangleSurroundsEdge test fails due to too many 0 dp" );
 571       }
 572
 573     return (cPQ*cQR >= 0.0) && (cQR*cRP >= 0.0) && (cRP*cPQ >= 0.0) && numZeros < 2;
 574   }
 575
 576 } // NAMESPACE