src/INTERP_KERNEL/TransformedTriangle.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 "VectorUtils.hxx"
  21
  22 #include <iostream>
  23 #include <fstream>
  24 #include <cassert>
  25 #include <algorithm>
  26 #include <functional>
  27 #include <iterator>
  28 #include <math.h>
  29 #include <vector>
  30
  31 namespace INTERP_KERNEL
  32 {
  33
  34   /**
  35    * \brief Class representing a circular order of a set of points around their barycenter.
  36    * It is used with the STL sort() algorithm to sort the point of the two polygons
  37    *
  38    */
  39   class ProjectedCentralCircularSortOrder
  40   {
  41   public:
  42
  43     /// Enumeration of different planes to project on when calculating order
  44     enum CoordType { XY, XZ, YZ };
  45
  46     /**
  47      * Constructor
  48      *
  49      * @param barycenter  double[3] containing the barycenter of the points to be compared
  50      * @param type        plane to project on when comparing. The comparison will not work if all the points are in a plane perpendicular
  51      *                    to the plane being projected on
  52      */
  53     ProjectedCentralCircularSortOrder(const double* barycenter, const CoordType type)
  54       : _aIdx((type == YZ) ? 1 : 0),
  55         _bIdx((type == XY) ? 1 : 2),
  56         _a(barycenter[_aIdx]),
  57         _b(barycenter[_bIdx])
  58     {
  59     }
  60
  61     /**
  62      * Comparison operator.
  63      * Compares the relative position between two points in their ordering around the barycenter.
  64      *
  65      * @param  pt1   a double[3] representing a point
  66      * @param  pt2   a double[3] representing a point
  67      * @return       true if the angle of the difference vector between pt1 and the barycenter is greater than that
  68      *               of the difference vector between pt2 and the barycenter.
  69      */
  70     bool operator()(const double* pt1, const double* pt2)
  71     {
  72       // calculate angles with the axis
  73       const double ang1 = atan2(pt1[_aIdx] - _a, pt1[_bIdx] - _b);
  74       const double ang2 = atan2(pt2[_aIdx] - _a, pt2[_bIdx] - _b);
  75
  76       return ang1 > ang2;
  77     }
  78
  79   private:
  80     /// index corresponding to first coordinate of plane on which points are projected
  81     const int _aIdx;
  82
  83     /// index corresponding to second coordinate of plane on which points are projected
  84     const int _bIdx;
  85
  86     /// value of first projected coordinate of the barycenter
  87     const double _a;
  88
  89     /// value of second projected coordinate of the barycenter
  90     const double _b;
  91   };
  92
  93   // ----------------------------------------------------------------------------------
  94   // TransformedTriangle PUBLIC
  95   // ----------------------------------------------------------------------------------
  96
  97   /**
  98    * Constructor
  99    *
 100    * The coordinates are copied to the internal member variables
 101    *
 102    * @param p   array of three doubles containing coordinates of P
 103    * @param q   array of three doubles containing coordinates of Q
 104    * @param r   array of three doubles containing coordinates of R
 105    */
 106   TransformedTriangle::TransformedTriangle(double* p, double* q, double* r)
 107     : _is_double_products_calculated(false),  _is_triple_products_calculated(false), _volume(0)
 108   {
 109
 110     for(int i = 0 ; i < 3 ; ++i)
 111       {
 112         // xyz coordinates
 113         _coords[5*P + i] = p[i];
 114         _coords[5*Q + i] = q[i];
 115         _coords[5*R + i] = r[i];
 116       }
 117
 118     // h coordinate
 119
 120     _coords[5*P + 3] = 1 - p[0] - p[1] - p[2];
 121     _coords[5*Q + 3] = 1 - q[0] - q[1] - q[2];
 122     _coords[5*R + 3] = 1 - r[0] - r[1] - r[2];
 123
 124     // H coordinate
 125     _coords[5*P + 4] = 1 - p[0] - p[1];
 126     _coords[5*Q + 4] = 1 - q[0] - q[1];
 127     _coords[5*R + 4] = 1 - r[0] - r[1];
 128
 129     // initialise rest of data
 130     preCalculateDoubleProducts();
 131
 132     preCalculateTriangleSurroundsEdge();
 133
 134     preCalculateTripleProducts();
 135
 136   }
 137
 138   /**
 139    * Destructor
 140    *
 141    * Deallocates the memory used to store the points of the polygons.
 142    * This memory is allocated in calculateIntersectionPolygons().
 143    */
 144   TransformedTriangle::~TransformedTriangle()
 145   {
 146     // delete elements of polygons
 147     for(std::vector<double*>::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it)
 148       {
 149         delete[] *it;
 150       }
 151     for(std::vector<double*>::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it)
 152       {
 153         delete[] *it;
 154       }
 155   }
 156
 157   /**
 158    * Calculates the volume of intersection between the triangle and the
 159    * unit tetrahedron.
 160    *
 161    * @return   volume of intersection of this triangle with unit tetrahedron,
 162    *            as described in Grandy
 163    *
 164    */
 165   double TransformedTriangle::calculateIntersectionVolume()
 166   {
 167     // check first that we are not below z - plane
 168     if(isTriangleBelowTetraeder())
 169       {
 170         LOG(2, " --- Triangle is below tetraeder - V = 0.0");
 171         return 0.0;
 172       }
 173
 174     // get the sign of the volume -  equal to the sign of the z-component of the normal
 175     // of the triangle, u_x * v_y - u_y * v_x, where u = q - p and v = r - p
 176     // if it is zero, the triangle is perpendicular to the z - plane and so the volume is zero
 177 //     const double uv_xy[4] =
 178 //       {
 179 //         _coords[5*Q] - _coords[5*P], _coords[5*Q + 1] - _coords[5*P + 1], // u_x, u_y
 180 //         _coords[5*R] - _coords[5*P], _coords[5*R + 1] - _coords[5*P + 1]  // v_x, v_y
 181 //       };
 182
 183 //     double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
 184     int sign = isTriangleInclinedToFacet( OXY );
 185
 186     if(sign == 0 )
 187       {
 188         LOG(2, " --- Triangle is perpendicular to z-plane - V = 0.0");
 189         return _volume = 0.0;
 190       }
 191
 192
 193     // normalize sign
 194     //sign = sign > 0.0 ? 1.0 : -1.0;
 195
 196     LOG(2, "-- Calculating intersection polygons ... ");
 197     calculateIntersectionPolygons();
 198
 199     double barycenter[3];
 200
 201     // calculate volume under A
 202     double volA = 0.0;
 203     if(_polygonA.size() > 2)
 204       {
 205         LOG(2, "---- Treating polygon A ... ");
 206         calculatePolygonBarycenter(A, barycenter);
 207         sortIntersectionPolygon(A, barycenter);
 208         volA = calculateVolumeUnderPolygon(A, barycenter);
 209         LOG(2, "Volume is " << sign * volA);
 210       }
 211
 212     double volB = 0.0;
 213     // if triangle is not in h = 0 plane, calculate volume under B
 214     if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ))
 215       {
 216         LOG(2, "---- Treating polygon B ... ");
 217
 218         calculatePolygonBarycenter(B, barycenter);
 219         sortIntersectionPolygon(B, barycenter);
 220         volB = calculateVolumeUnderPolygon(B, barycenter);
 221         LOG(2, "Volume is " << sign * volB);
 222       }
 223
 224     LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********");
 225
 226     return _volume = sign * (volA + volB);
 227
 228   }
 229
 230   // ----------------------------------------------------------------------------------
 231   // TransformedTriangle PRIVATE
 232   // ----------------------------------------------------------------------------------
 233
 234   /**
 235    * Calculates the intersection polygons A and B, performing the intersection tests
 236    * and storing the corresponding points in the vectors _polygonA and _polygonB.
 237    *
 238    * @post _polygonA contains the intersection polygon A and _polygonB contains the
 239    *       intersection polygon B.
 240    *
 241    */
 242   void TransformedTriangle::calculateIntersectionPolygons()
 243   {
 244     assert(_polygonA.size() == 0);
 245     assert(_polygonB.size() == 0);
 246     // avoid reallocations in push_back() by pre-allocating enough memory
 247     // we should never have more than 20 points
 248     _polygonA.reserve(20);
 249     _polygonB.reserve(20);
 250     // -- surface intersections
 251     // surface - edge
 252     for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
 253       {
 254         if(testSurfaceEdgeIntersection(edge))
 255           {
 256             double* ptA = new double[3];
 257             calcIntersectionPtSurfaceEdge(edge, ptA);
 258             _polygonA.push_back(ptA);
 259             LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A ");
 260             if(edge >= XY)
 261               {
 262                 double* ptB = new double[3];
 263                 copyVector3(ptA, ptB);
 264                 _polygonB.push_back(ptB);
 265                 LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B ");
 266               }
 267
 268           }
 269       }
 270
 271     // surface - ray
 272     for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
 273       {
 274         if(testSurfaceRayIntersection(corner))
 275           {
 276             double* ptB = new double[3];
 277             copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
 278             _polygonB.push_back(ptB);
 279             LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B");
 280           }
 281       }
 282
 283     // -- segment intersections
 284     for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
 285       {
 286
 287         bool isZero[NO_DP];
 288
 289         // check beforehand which double-products are zero
 290         for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
 291           {
 292             isZero[dp] = (calcStableC(seg, dp) == 0.0);
 293           }
 294
 295         // segment - facet
 296         for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
 297           {
 298             // is this test worth it?
 299             const bool doTest =
 300               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
 301               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
 302               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
 303
 304             if(doTest && testSegmentFacetIntersection(seg, facet))
 305               {
 306                 double* ptA = new double[3];
 307                 calcIntersectionPtSegmentFacet(seg, facet, ptA);
 308                 _polygonA.push_back(ptA);
 309                 LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
 310                 if(facet == XYZ)
 311                   {
 312                     double* ptB = new double[3];
 313                     copyVector3(ptA, ptB);
 314                     _polygonB.push_back(ptB);
 315                     LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B");
 316                   }
 317
 318               }
 319           }
 320
 321         // segment - edge
 322         for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
 323           {
 324             const DoubleProduct edge_dp = DoubleProduct(edge);
 325
 326             if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
 327               {
 328                 double* ptA = new double[3];
 329                 calcIntersectionPtSegmentEdge(seg, edge, ptA);
 330                 _polygonA.push_back(ptA);
 331                 LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
 332                 if(edge >= XY)
 333                   {
 334                     double* ptB = new double[3];
 335                     copyVector3(ptA, ptB);
 336                     _polygonB.push_back(ptB);
 337                   }
 338               }
 339           }
 340
 341         // segment - corner
 342         for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
 343           {
 344             const bool doTest =
 345               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
 346               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
 347               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
 348
 349             if(doTest && testSegmentCornerIntersection(seg, corner))
 350               {
 351                 double* ptA = new double[3];
 352                 copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
 353                 _polygonA.push_back(ptA);
 354                 LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
 355                 if(corner != O)
 356                   {
 357                     double* ptB = new double[3];
 358                     _polygonB.push_back(ptB);
 359                     copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
 360                     LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B");
 361                   }
 362               }
 363           }
 364
 365             // segment - ray
 366             for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
 367               {
 368                 if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner))
 369                   {
 370                     double* ptB = new double[3];
 371                     copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
 372                     _polygonB.push_back(ptB);
 373                     LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
 374                   }
 375               }
 376
 377             // segment - halfstrip
 378             for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
 379               {
 380
 381 #if 0
 382                 const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX
 383                 const bool doTest =
 384                   !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] &&
 385                   !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]];
 386
 387
 388                 if(doTest && testSegmentHalfstripIntersection(seg, edge))
 389 #endif
 390                   if(testSegmentHalfstripIntersection(seg, edge))
 391                     {
 392                       double* ptB = new double[3];
 393                       calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
 394                       _polygonB.push_back(ptB);
 395                       LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
 396                     }
 397               }
 398       }
 399
 400         // inclusion tests
 401         for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
 402           {
 403             // { XYZ - inclusion only possible if in Tetrahedron?
 404             // tetrahedron
 405             if(testCornerInTetrahedron(corner))
 406               {
 407                 double* ptA = new double[3];
 408                 copyVector3(&_coords[5*corner], ptA);
 409                 _polygonA.push_back(ptA);
 410                 LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
 411               }
 412
 413             // XYZ - plane
 414             if(testCornerOnXYZFacet(corner))
 415               {
 416                 double* ptB = new double[3];
 417                 copyVector3(&_coords[5*corner], ptB);
 418                 _polygonB.push_back(ptB);
 419                 LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B");
 420               }
 421
 422             // projection on XYZ - facet
 423             if(testCornerAboveXYZFacet(corner))
 424               {
 425                 double* ptB = new double[3];
 426                 copyVector3(&_coords[5*corner], ptB);
 427                 ptB[2] = 1 - ptB[0] - ptB[1];
 428                 assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
 429                 _polygonB.push_back(ptB);
 430                 LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B");
 431               }
 432
 433           }
 434
 435       }
 436
 437     /**
 438      * Calculates the barycenters of the given intersection polygon.
 439      *
 440      * @pre  the intersection polygons have been calculated with calculateIntersectionPolygons()
 441      *
 442      * @param poly        one of the two intersection polygons
 443      * @param barycenter  array of three doubles where barycenter is stored
 444      *
 445      */
 446     void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter)
 447     {
 448       LOG(3,"--- Calculating polygon barycenter");
 449
 450       // get the polygon points
 451       std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
 452
 453       // calculate barycenter
 454       const int m = polygon.size();
 455
 456       for(int j = 0 ; j < 3 ; ++j)
 457         {
 458           barycenter[j] = 0.0;
 459         }
 460
 461       if(m != 0)
 462         {
 463           for(int i = 0 ; i < m ; ++i)
 464             {
 465               const double* pt = polygon[i];
 466               for(int j = 0 ; j < 3 ; ++j)
 467                 {
 468                   barycenter[j] += pt[j] / double(m);
 469                 }
 470             }
 471         }
 472       LOG(3,"Barycenter is " << vToStr(barycenter));
 473     }
 474
 475     /**
 476      * Sorts the given intersection polygon in circular order around its barycenter.
 477      * @pre  the intersection polygons have been calculated with calculateIntersectionPolygons()
 478      * @post the vertices in _polygonA and _polygonB are sorted in circular order around their
 479      *       respective barycenters
 480      *
 481      * @param poly        one of the two intersection polygons
 482      * @param barycenter  array of three doubles with the coordinates of the barycenter
 483      *
 484      */
 485     void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter)
 486     {
 487       LOG(3,"--- Sorting polygon ...");
 488
 489       using INTERP_KERNEL::ProjectedCentralCircularSortOrder;
 490       typedef ProjectedCentralCircularSortOrder SortOrder; // change is only necessary here and in constructor
 491       typedef SortOrder::CoordType CoordType;
 492
 493       // get the polygon points
 494       std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
 495
 496       if(polygon.size() == 0)
 497         return;
 498
 499       // determine type of sorting
 500       CoordType type = SortOrder::XY;
 501       if(poly == A && !isTriangleInclinedToFacet( OXY )) // B is on h = 0 plane -> ok
 502         {
 503           // NB : the following test is never true if we have eliminated the
 504           // triangles parallel to x == 0 and y == 0 in calculateIntersectionVolume().
 505           // We keep the test here anyway, to avoid interdependency.
 506
 507           // is triangle inclined to x == 0 ?
 508           if(isTriangleInclinedToFacet(OZX))
 509             {
 510               type = SortOrder::XZ;
 511             }
 512           else //if(isTriangleParallelToFacet(OYZ))
 513             {
 514               type = SortOrder::YZ;
 515             }
 516         }
 517
 518       // create order object
 519       SortOrder order(barycenter, type);
 520
 521       // sort vector with this object
 522       // NB : do not change place of first object, with respect to which the order
 523       // is defined
 524       sort((polygon.begin()), polygon.end(), order);
 525
 526       LOG(3,"Sorted polygon is ");
 527       for(size_t i = 0 ; i < polygon.size() ; ++i)
 528         {
 529           LOG(3,vToStr(polygon[i]));
 530         }
 531
 532     }
 533
 534     /**
 535      * Calculates the volume between the given polygon and the z = 0 plane.
 536      *
 537      * @pre  the intersection polygones have been calculated with calculateIntersectionPolygons(),
 538      *       and they have been sorted in circular order with sortIntersectionPolygons(void)
 539      *
 540      * @param poly        one of the two intersection polygons
 541      * @param barycenter  array of three doubles with the coordinates of the barycenter
 542      * @return           the volume between the polygon and the z = 0 plane
 543      *
 544      */
 545     double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter)
 546     {
 547       LOG(2,"--- Calculating volume under polygon");
 548
 549       // get the polygon points
 550       std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
 551
 552       double vol = 0.0;
 553       const int m = polygon.size();
 554
 555       for(int i = 0 ; i < m ; ++i)
 556         {
 557           const double* ptCurr = polygon[i];  // pt "i"
 558           const double* ptNext = polygon[(i + 1) % m]; // pt "i+1" (pt m == pt 0)
 559
 560           const double factor1 = ptCurr[2] + ptNext[2] + barycenter[2];
 561           const double factor2 =
 562             ptCurr[0]*(ptNext[1] - barycenter[1])
 563             + ptNext[0]*(barycenter[1] - ptCurr[1])
 564             + barycenter[0]*(ptCurr[1] - ptNext[1]);
 565           vol += (factor1 * factor2) / 6.0;
 566         }
 567
 568       LOG(2,"Abs. Volume is " << vol);
 569       return vol;
 570     }
 571
 572
 573     ////////////////////////////////////////////////////////////////////////////////////
 574     // Detection of (very) degenerate cases                                /////////////
 575     ////////////////////////////////////////////////////////////////////////////////////
 576
 577     /**
 578      * Checks if the triangle lies in the plane of a given facet
 579      *
 580      * @param facet     one of the facets of the tetrahedron
 581      * @return         true if PQR lies in the plane of the facet, false if not
 582      */
 583     bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) const
 584     {
 585
 586       // coordinate to check
 587       const int coord = static_cast<int>(facet);
 588
 589       for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
 590         {
 591           if(_coords[5*c + coord] != 0.0)
 592             {
 593               return false;
 594             }
 595         }
 596
 597       return true;
 598     }
 599
 600     /**
 601      * Checks if the triangle is parallel to the given facet
 602      *
 603      * @param facet  one of the facets of the unit tetrahedron
 604      * @return       true if triangle is parallel to facet, false if not
 605      */
 606     bool TransformedTriangle::isTriangleParallelToFacet(const TetraFacet facet) const
 607     {
 608       // coordinate to check
 609       const int coord = static_cast<int>(facet);
 610       return (_coords[5*P + coord] == _coords[5*Q + coord]) && (_coords[5*P + coord] == _coords[5*R + coord]);
 611     }
 612
 613     /**
 614      * Checks if the triangle is not perpedicular to the given facet
 615      *
 616      * @param facet  one of the facets of the unit tetrahedron
 617      * @return       zero if the triangle is perpendicular to the facet,
 618      *               else 1 or -1 depending on the sign of cross product of facet edges
 619      */
 620     int TransformedTriangle::isTriangleInclinedToFacet(const TetraFacet facet) const
 621     {
 622       // coordinate to check
 623       const int coord = static_cast<int>(facet);
 624       const int ind1 = ( coord+1 ) % 3, ind2 = ( coord+2 ) % 3;
 625       const double uv_xy[4] =
 626         {
 627           // u_x, u_y
 628           _coords[5*Q+ind1] - _coords[5*P+ind1], _coords[5*Q+ind2] - _coords[5*P+ind2],
 629           // v_x, v_y
 630           _coords[5*R+ind1] - _coords[5*P+ind1], _coords[5*R+ind2] - _coords[5*P+ind2]
 631         };
 632
 633       double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
 634       return (sign < 0.) ? -1 : (sign > 0.) ? 1 : 0;
 635     }
 636
 637     /**
 638      * Determines whether the triangle is below the z-plane.
 639      *
 640      * @return true if the z-coordinate of the three corners of the triangle are all less than 0, false otherwise.
 641      */
 642     bool TransformedTriangle::isTriangleBelowTetraeder() const
 643     {
 644       for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
 645         {
 646           // check z-coords for all points
 647           if(_coords[5*c + 2] >= 0.0)
 648             {
 649               return false;
 650             }
 651         }
 652       return true;
 653     }
 654
 655     /**
 656      * Prints the coordinates of the triangle to std::cout
 657      *
 658      */
 659     void TransformedTriangle::dumpCoords() const
 660     {
 661       std::cout << "Coords : ";
 662       for(int i = 0 ; i < 3; ++i)
 663         {
 664           std::cout << vToStr(&_coords[5*i]) << ",";
 665         }
 666       std::cout << std::endl;
 667     }
 668
 669   } // NAMESPACE