src/INTERP_KERNEL/TransformedTriangleIntersect.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   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.
   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 #include "TransformedTriangle.hxx"
  20 #include <iostream>
  21 #include <fstream>
  22 #include <cassert>
  23 #include <cmath>
  24 #include "VectorUtils.hxx"
  25
  26 namespace INTERP_KERNEL
  27 {
  28
  29   // ----------------------------------------------------------------------------------
  30   //  Correspondance tables describing all the variations of formulas.
  31   // ----------------------------------------------------------------------------------
  32
  33   /// \brief Correspondance between facets and double products.
  34   ///
  35   /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
  36   const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_SEG_FACET_INTERSECTION[12] =
  37     {
  38       C_XH, C_XY, C_ZX, // OYZ
  39       C_YH, C_YZ, C_XY, // OZX
  40       C_ZH, C_ZX, C_YZ, // OXY
  41       C_XH, C_YH, C_ZH  // XYZ
  42     };
  43
  44   /// \brief Signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION.
  45   ///
  46   /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
  47   const double TransformedTriangle::SIGN_FOR_SEG_FACET_INTERSECTION[12] =
  48     {
  49       1.0, 1.0, -1.0,
  50       1.0, 1.0, -1.0,
  51       1.0, 1.0, -1.0,
  52       1.0, 1.0,  1.0
  53     };
  54
  55   /// \brief Coordinates of corners of tetrahedron.
  56   ///
  57   /// Use 3*Corner + coordinate as index
  58   const double TransformedTriangle::COORDS_TET_CORNER[12] =
  59     {
  60       0.0, 0.0, 0.0,
  61       1.0, 0.0, 0.0,
  62       0.0, 1.0, 0.0,
  63       0.0, 0.0, 1.0
  64     };
  65
  66   /// \brief Indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION
  67   /// for the calculation of the coordinates (x,y,z) of the intersection points
  68   /// for Segment-Facet and Segment-Edge intersections.
  69   ///
  70   /// Use 3*facet + coordinate as index. -1 indicates that the coordinate is 0.
  71   const int TransformedTriangle::DP_INDEX[12] =
  72     {
  73       // x, y, z
  74       -1, 1, 2,  // OYZ
  75       5, -1, 4,  // OZX
  76       7, 8, -1,  // OXY
  77       9, 10, 11  // XYZ
  78     };
  79
  80   /// \brief Correspondance edge - corners.
  81   ///
  82   /// Gives the two corners associated with each edge
  83   /// Use 2*edge + {0, 1} as index
  84   const TransformedTriangle::TetraCorner TransformedTriangle::CORNERS_FOR_EDGE[12] =
  85     {
  86       O, X, // OX
  87       O, Y, // OY
  88       O, Z, // OZ
  89       X, Y, // XY
  90       Y, Z, // YZ
  91       Z, X  // ZX
  92     };
  93
  94   /// \brief Correspondance edge - facets.
  95   ///
  96   /// Gives the two facets shared by and edge. Use 2*facet + {0, 1} as index
  97   const TransformedTriangle::TetraFacet TransformedTriangle::FACET_FOR_EDGE[12] =
  98     {
  99       OXY, OZX, // OX
 100       OXY, OYZ, // OY
 101       OZX, OYZ, // OZ
 102       OXY, XYZ, // XY
 103       OYZ, XYZ, // YZ
 104       OZX, XYZ  // ZX
 105     };
 106
 107   /// \brief Correspondance corners - edges.
 108   ///
 109   /// Gives edges meeting at a given corner. Use 3*corner + {0,1,2} as index
 110   const TransformedTriangle::TetraEdge TransformedTriangle::EDGES_FOR_CORNER[12] =
 111     {
 112       OX, OY, OZ, // O
 113       OX, XY, ZX, // X
 114       OY, XY, YZ, // Y
 115       OZ, ZX, YZ  // Z
 116     };
 117
 118   /// \brief Double products to use in halfstrip intersection tests.
 119   ///
 120   /// Use 4*(offset_edge) + {0,1,2,3} as index. offset_edge = edge - 3  (so that XY -> 0, YZ -> 1, ZX -> 2)
 121   /// Entries with offset 0 and 1 are for the first condition (positive product)
 122   /// and those with offset 2 and 3 are for the second condition (negative product).
 123   const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_HALFSTRIP_INTERSECTION[12] =
 124     {
 125       C_10, C_01, C_ZH, C_10, // XY
 126       C_01, C_XY, C_XH, C_01, // YZ
 127       C_XY, C_10, C_YH, C_XY  // ZX
 128     };
 129
 130   /// \brief Double products to use in segment-ray test.
 131   ///
 132   /// 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)
 133   /// Entries with offset 0 are for first condition (zero double product) and the rest are for condition 3 (in the same
 134   /// order as in the article)
 135   const TransformedTriangle::DoubleProduct TransformedTriangle::DP_SEGMENT_RAY_INTERSECTION[21] =
 136     {
 137       C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
 138       C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_10, // Y
 139       C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_01  // Z
 140     };
 141
 142   /**
 143    * Calculates the point of intersection between the given edge of the tetrahedron and the
 144    * triangle PQR. (Grandy, eq [22])
 145    *
 146    * @pre   testSurfaceEdgeIntersection(edge) returns true
 147    * @param edge   edge of tetrahedron
 148    * @param pt     array of three doubles in which to store the coordinates of the intersection point
 149    */
 150   void TransformedTriangle::calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const
 151   {
 152     assert(edge < H01);
 153
 154     // barycentric interpolation between points A and B
 155     // : (x,y,z)* = (1-alpha)*A + alpha*B where
 156     // alpha = t_A / (t_A - t_B)
 157
 158     const TetraCorner corners[2] =
 159       {
 160         CORNERS_FOR_EDGE[2*edge],
 161         CORNERS_FOR_EDGE[2*edge + 1]
 162       };
 163
 164     // calculate alpha
 165     const double tA = calcStableT(corners[0]);
 166     const double tB = calcStableT(corners[1]);
 167     const double alpha = tA / (tA - tB);
 168
 169     // calculate point
 170     LOG(4, "corner A = " << corners[0] << " corner B = " << corners[1] );
 171     LOG(4, "tA = " << tA << " tB = " << tB << " alpha= " << alpha );
 172     for(int i = 0; i < 3; ++i)
 173       {
 174
 175         pt[i] = (1 - alpha) * COORDS_TET_CORNER[3*corners[0] + i] +
 176           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    * If the OPTIMIZE is defined, it does not do the test the double product that should be zero.
 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             const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx];
 346             const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
 347             c[j] = dpIdx < 0.0 ? 0.0 : sign * calcStableC(seg, dp);
 348           }
 349
 350         // pt[i] = (c1*s1 + c2*s2) / (s1^2 + s2^2)
 351
 352         pt[i] = (c[0] * s[0] + c[1] * s[1]) / denominator;
 353
 354         // strange bug with -O2 enabled : assertion fails when we don't have the following
 355         // trace - line
 356         //std::cout << "pt[i] = " << pt[i] << std::endl;
 357         //assert(pt[i] >= 0.0); // check we are in tetraeder
 358         //assert(pt[i] <= 1.0);
 359
 360       }
 361   }
 362
 363
 364   /**
 365    * Tests if the given segment of the triangle intersects the given corner of the tetrahedron.
 366    * (Grandy, eq. [21]). If OPTIMIZE is defined, the double products that should be zero are not verified.
 367    *
 368    * @param seg    segment of the triangle
 369    * @param corner corner of the tetrahedron
 370    * @return      true if the segment intersects the corner
 371    */
 372   bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const
 373   {
 374
 375
 376     // facets meeting at a given corner
 377     static const TetraFacet FACETS_FOR_CORNER[12] =
 378       {
 379         OXY, OYZ, OZX, // O
 380         OZX, OXY, XYZ, // X
 381         OYZ, XYZ, OXY, // Y
 382         OZX, XYZ, OYZ  // Z
 383       };
 384
 385     // check segment intersect a facet
 386     for(int i = 0 ; i < 3 ; ++i)
 387       {
 388         const TetraFacet facet = FACETS_FOR_CORNER[3*corner + i];
 389         if(testSegmentIntersectsFacet(seg, facet))
 390           {
 391             return true;
 392           }
 393       }
 394
 395     return false;
 396   }
 397
 398   /**
 399    * Tests if the given segment of the triangle intersects the half-strip above the
 400    * given edge of the h = 0 plane. (Grandy, eq. [30])
 401    *
 402    * @param seg    segment of the triangle
 403    * @param edge   edge of the h = 0 plane of the tetrahedron (XY, YZ, ZX)
 404    * @return      true if the upwards ray from the corner intersects the triangle.
 405    */
 406   bool TransformedTriangle::testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edge)
 407   {
 408     // get right index here to avoid "filling out" array
 409     const int edgeIndex = static_cast<int>(edge) - 3;
 410
 411     // double products used in test
 412     // products 1 and 2 for each edge -> first condition in Grandy [30]
 413     // products 3 and 4 for each edge -> third condition
 414     // NB : some uncertainty whether these last are correct
 415     static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12] =
 416       {
 417         C_10, C_01, C_ZH, C_10, // XY
 418         C_01, C_XY, C_XH, C_01, // YZ
 419         C_XY, C_10, C_YH, C_XY  // ZX
 420       };
 421
 422     // facets to use in second condition (S_m)
 423     static const TetraFacet FACET_FOR_HALFSTRIP_INTERSECTION[3] =
 424       {
 425         NO_TET_FACET, // XY -> special case : test with plane H = 0
 426         OYZ, // YZ
 427         OZX  // ZX
 428       };
 429
 430     const double cVals[4] =
 431       {
 432         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex]),
 433         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 1]),
 434         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 2]),
 435         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 3])
 436       };
 437
 438     const TetraFacet facet =  FACET_FOR_HALFSTRIP_INTERSECTION[edgeIndex];
 439
 440
 441     // special case : facet H = 0
 442     const bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet);
 443     LOG(4, "Halfstrip tests (" << seg << ", " << edge << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) );
 444     LOG(4, "c2 = " << cVals[2] << ", c3 = " << cVals[3] );
 445
 446     return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0);
 447   }
 448
 449   /**
 450    * Calculates the point of intersection between the given segment of the triangle
 451    * and the halfstrip above the given edge of the tetrahedron. (Grandy, eq. [31])
 452    *
 453    * @pre   testSegmentHalfstripIntersection(seg, edge) returns true
 454    *
 455    * @param seg    segment of the triangle
 456    * @param edge   edge of the tetrahedron defining the halfstrip
 457    * @param pt     array of three doubles in which to store the coordinates of the intersection point
 458    */
 459   void TransformedTriangle::calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const
 460   {
 461     assert(edge > OZ);
 462     assert(edge < H01);
 463
 464     // get right index here to avoid "filling out" array
 465     const int edgeIndex = static_cast<int>(edge) - 3;
 466     assert(edgeIndex >= 0);
 467     assert(edgeIndex < 3);
 468
 469     // Barycentric interpolation on the edge
 470     // for edge AB : (x,y,z)* = (1-alpha) * A + alpha * B
 471     // where alpha = cB / (cB - cA)
 472
 473     const double cA = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex]);
 474     const double cB = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 1]);
 475     assert(cA != cB);
 476
 477     const double alpha = cA / (cA - cB);
 478
 479     for(int i = 0; i < 3; ++i)
 480       {
 481         const TetraCorner corners[2] =
 482           {
 483             CORNERS_FOR_EDGE[2*edge],
 484             CORNERS_FOR_EDGE[2*edge + 1]
 485           };
 486
 487         const double cornerCoords[2] =
 488           {
 489             COORDS_TET_CORNER[3*corners[0] + i],
 490             COORDS_TET_CORNER[3*corners[1] + i]
 491           };
 492
 493         pt[i] = (1 - alpha) * cornerCoords[0] + alpha * cornerCoords[1];
 494         LOG(6, pt[i] );
 495         assert(pt[i] >= 0.0);
 496         assert(pt[i] <= 1.0);
 497       }
 498     assert(epsilonEqualRelative(pt[0] + pt[1] + pt[2], 1.0));
 499   }
 500
 501   /**
 502    * Tests if the given segment of triangle PQR intersects the ray pointing
 503    * in the upwards z - direction from the given corner of the tetrahedron. (Grandy eq. [29])
 504    * If OPTIMIZE is defined, the double product that should be zero is not verified.
 505    *
 506    * @param seg    segment of the triangle PQR
 507    * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z)
 508    * @return      true if the upwards ray from the corner intersects the segment.
 509    */
 510   bool TransformedTriangle::testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const
 511   {
 512     assert(corner == X || corner == Y || corner == Z);
 513     LOG(4, "Testing seg - ray intersection for seg = " << seg << ", corner = " << corner );
 514
 515     // readjust index since O is not used
 516     const int cornerIdx = static_cast<int>(corner) - 1;
 517
 518     // facets to use
 519     //? not sure this is correct
 520     static const TetraFacet FIRST_FACET_SEGMENT_RAY_INTERSECTION[3] =
 521       {
 522         OZX, // X
 523         OYZ, // Y
 524         OZX, // Z
 525       };
 526
 527     // cond 2
 528     const bool cond21 = testSegmentIntersectsFacet(seg, FIRST_FACET_SEGMENT_RAY_INTERSECTION[cornerIdx]);
 529     const bool cond22  = (corner == Z) ? testSegmentIntersectsFacet(seg, OYZ) : testSegmentIntersectsHPlane(seg);
 530
 531     if(!(cond21 || cond22))
 532       {
 533         LOG(4, "SR fails at cond 2 : cond21 = " << cond21 << ", cond22 = " << cond22 );
 534         return false;
 535       }
 536
 537     // cond 3
 538     const double cVals[6] =
 539       {
 540         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 1]),
 541         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 2]),
 542         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 3]),
 543         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 4]),
 544         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 5]),
 545         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 6]),
 546       };
 547
 548     // cond. 3
 549     if(( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) >= 0.0)
 550       {
 551         LOG(4, "SR fails at cond 3 : " << (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5]  );
 552       }
 553     return ( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) < 0.0;
 554
 555   }
 556
 557   // /////////////////////////////////////////////////////////////////////////////////
 558   //  Utility methods used in intersection tests                       ///////////////
 559   // /////////////////////////////////////////////////////////////////////////////////
 560   /**
 561    * Tests if the triangle PQR surrounds the axis on which the
 562    * given edge of the tetrahedron lies.
 563    *
 564    * @param edge   edge of tetrahedron
 565    * @return      true if PQR surrounds edge, false if not (see Grandy, eq. [53])
 566    */
 567   bool TransformedTriangle::testTriangleSurroundsEdge(const TetraEdge edge) const
 568   {
 569     // NB DoubleProduct enum corresponds to TetraEdge enum according to Grandy, table III
 570     // so we can use the edge directly
 571
 572     const double cPQ = calcStableC(PQ, DoubleProduct(edge));
 573     const double cQR = calcStableC(QR, DoubleProduct(edge));
 574     const double cRP = calcStableC(RP, DoubleProduct(edge));
 575
 576     LOG(5, "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" );
 577
 578     // if two or more c-values are zero we disallow x-edge intersection
 579     // Grandy, p.446
 580     const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
 581
 582     if(numZeros >= 2 )
 583       {
 584         LOG(5, "TriangleSurroundsEdge test fails due to too many 0 dp" );
 585       }
 586
 587     return (cPQ*cQR >= 0.0) && (cQR*cRP >= 0.0) && (cRP*cPQ >= 0.0) && numZeros < 2;
 588   }
 589
 590 } // NAMESPACE