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