src/INTERP_KERNEL/TransformedTriangle.cxx

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