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       // 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(auto& it: _polygonA)
 151       delete[] it;
 152     for(auto& it: _polygonB)
 153       delete[] it;
 154   }
 155
 156   /**
 157    * Calculates the volume of intersection between the triangle and the
 158    * unit tetrahedron.
 159    *
 160    * @return   volume of intersection of this triangle with unit tetrahedron,
 161    *            as described in Grandy
 162    *
 163    */
 164   double TransformedTriangle::calculateIntersectionVolume()
 165   {
 166     // check first that we are not below z - plane
 167     if(isTriangleBelowTetraeder())
 168       {
 169         LOG(2, " --- Triangle is below tetraeder - V = 0.0");
 170         return 0.0;
 171       }
 172
 173     // get the sign of the volume -  equal to the sign of the z-component of the normal
 174     // of the triangle, u_x * v_y - u_y * v_x, where u = q - p and v = r - p
 175     // if it is zero, the triangle is perpendicular to the z - plane and so the volume is zero
 176 //     const double uv_xy[4] =
 177 //       {
 178 //         _coords[5*Q] - _coords[5*P], _coords[5*Q + 1] - _coords[5*P + 1], // u_x, u_y
 179 //         _coords[5*R] - _coords[5*P], _coords[5*R + 1] - _coords[5*P + 1]  // v_x, v_y
 180 //       };
 181
 182 //     double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
 183     int sign = isTriangleInclinedToFacet( OXY );
 184
 185     if(sign == 0 )
 186       {
 187         LOG(2, " --- Triangle is perpendicular to z-plane - V = 0.0");
 188         return _volume = 0.0;
 189       }
 190
 191
 192     // normalize sign
 193     //sign = sign > 0.0 ? 1.0 : -1.0;
 194
 195     LOG(2, "-- Calculating intersection polygons ... ");
 196     calculateIntersectionAndProjectionPolygons();
 197
 198     double barycenter[3];
 199
 200     // calculate volume under A
 201     double volA = 0.0;
 202     if(_polygonA.size() > 2)
 203       {
 204         LOG(2, "---- Treating polygon A ... ");
 205 #if LOG_LEVEL > 0
 206         LOG(3, "  --- Final points in polygon A");
 207         for(const auto& pt:  _polygonA)
 208           LOG(3,vToStr(pt));
 209 #endif
 210         calculatePolygonBarycenter(A, barycenter);
 211         sortIntersectionPolygon(A, barycenter);
 212         volA = calculateVolumeUnderPolygon(A, barycenter);
 213         LOG(2, "Volume is " << sign * volA);
 214       }
 215
 216     double volB = 0.0;
 217     // if triangle is not in h = 0 plane, calculate volume under B
 218     if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ)) // second condition avoids double counting in case triangle fully included in h=0 facet
 219       {
 220         LOG(2, "---- Treating polygon B ... ");
 221 #if LOG_LEVEL > 0
 222         LOG(3, "  --- Final points in polygon B");
 223         for(const auto& pt:  _polygonB)
 224           LOG(3,vToStr(pt));
 225 #endif
 226         calculatePolygonBarycenter(B, barycenter);
 227         sortIntersectionPolygon(B, barycenter);
 228         volB = calculateVolumeUnderPolygon(B, barycenter);
 229         LOG(2, "Volume is " << sign * volB);
 230       }
 231
 232 #if LOG_LEVEL >= 2
 233     LOG(2, "############ Triangle :")
 234     dumpCoords();
 235     LOG(2, "vol A = " << volA);
 236     LOG(2, "vol B = " << volB);
 237     LOG(2, "TOTAL = " << sign*(volA+volB));
 238 #endif
 239
 240     return _volume = sign * (volA + volB);
 241
 242   }
 243
 244   /**
 245    * Calculates the volume of intersection between the triangle and the
 246    * unit tetrahedron.
 247    *
 248    * @return   volume of intersection of this triangle with unit tetrahedron,
 249    *            as described in Grandy
 250    *
 251    */
 252   double TransformedTriangle::calculateIntersectionSurface(TetraAffineTransform* tat)
 253   {
 254     // check first that we are not below z - plane
 255     if(isTriangleBelowTetraeder())
 256       {
 257         LOG(2, " --- Triangle is below tetraeder - V = 0.0");
 258         return 0.0;
 259       }
 260
 261     LOG(2, "-- Calculating intersection polygon ... ");
 262     calculateIntersectionPolygon();
 263
 264     _volume = 0.;
 265     if(_polygonA.size() > 2) {
 266       double barycenter[3];
 267       calculatePolygonBarycenter(A, barycenter);
 268       sortIntersectionPolygon(A, barycenter);
 269       const std::size_t nbPoints = _polygonA.size();
 270       for(std::size_t i = 0 ; i < nbPoints ; ++i)
 271         tat->reverseApply(_polygonA[i], _polygonA[i]);
 272       _volume = calculateSurfacePolygon();
 273     }
 274
 275     return _volume;
 276   }
 277
 278   // ----------------------------------------------------------------------------------
 279   // TransformedTriangle PROTECTED
 280   // ----------------------------------------------------------------------------------
 281
 282   /**
 283    * Calculates the intersection polygons A and B, performing the intersection tests
 284    * and storing the corresponding points in the vectors _polygonA and _polygonB.
 285    *
 286    * @post _polygonA contains the intersection polygon A and _polygonB contains the
 287    *       intersection polygon B.
 288    *
 289    */
 290   void TransformedTriangle::calculateIntersectionAndProjectionPolygons()
 291   {
 292 #if LOG_LEVEL >= 3
 293     std::cout << " @@@@@@@@ COORDS @@@@@@  " << std::endl;
 294     dumpCoords();
 295 #endif
 296
 297     assert(_polygonA.size() == 0);
 298     assert(_polygonB.size() == 0);
 299     // avoid reallocations in push_back() by pre-allocating enough memory
 300     // we should never have more than 20 points
 301     _polygonA.reserve(20);
 302     _polygonB.reserve(20);
 303     // -- surface intersections
 304     // surface - edge
 305     for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
 306       {
 307         if(testSurfaceEdgeIntersection(edge))
 308           {
 309             double* ptA = new double[3];
 310             calcIntersectionPtSurfaceEdge(edge, ptA);
 311             _polygonA.push_back(ptA);
 312             LOG(3,"Surface-edge (edge " << strTE(edge) << "): " << vToStr(ptA) << " added to A ");
 313             if(edge >= XY)
 314               {
 315                 double* ptB = new double[3];
 316                 copyVector3(ptA, ptB);
 317                 _polygonB.push_back(ptB);
 318                 LOG(3,"Surface-edge (edge " << strTE(edge) << "): " << vToStr(ptB) << " added to B ");
 319               }
 320
 321           }
 322       }
 323
 324     // surface - ray
 325     for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
 326       {
 327         if(testSurfaceRayIntersection(corner))
 328           {
 329             double* ptB = new double[3];
 330             copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
 331             _polygonB.push_back(ptB);
 332             LOG(3,"Surface-ray (corner " << strTC(corner) << "): " << vToStr(ptB) << " added to B");
 333           }
 334       }
 335
 336     // -- segment intersections
 337     for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
 338       {
 339
 340         bool isZero[NO_DP];
 341
 342         // check beforehand which double-products are zero.
 343         for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
 344           isZero[dp] = (calcStableC(seg, dp) == 0.0);
 345
 346         // segment - facet
 347         for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
 348           {
 349             // is this test worth it?
 350             const bool doTest =
 351                 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
 352                 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
 353                 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
 354
 355             if(doTest && testSegmentFacetIntersection(seg, facet))
 356               {
 357                 double* ptA = new double[3];
 358                 calcIntersectionPtSegmentFacet(seg, facet, ptA);
 359                 _polygonA.push_back(ptA);
 360                 LOG(3,"Segment-facet (facet " << strTF(facet) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A");
 361                 if(facet == XYZ)
 362                   {
 363                     double* ptB = new double[3];
 364                     copyVector3(ptA, ptB);
 365                     _polygonB.push_back(ptB);
 366                     LOG(3,"Segment-facet (facet " << strTF(facet) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B");
 367                   }
 368
 369               }
 370           }
 371
 372         // segment - edge
 373         for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
 374           {
 375             const DoubleProduct edge_dp = DoubleProduct(edge);
 376
 377             if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
 378               {
 379                 double* ptA = new double[3];
 380                 calcIntersectionPtSegmentEdge(seg, edge, ptA);
 381                 _polygonA.push_back(ptA);
 382                 LOG(3,"Segment-edge (edge " << strTE(edge) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A");
 383                 if(edge >= XY)
 384                   {
 385                     double* ptB = new double[3];
 386                     copyVector3(ptA, ptB);
 387                     _polygonB.push_back(ptB);
 388                     LOG(3,"Segment-edge (edge " << strTE(edge) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to B");
 389                   }
 390               }
 391           }
 392
 393         // segment - corner
 394         for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
 395           {
 396             const bool doTest =
 397                 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
 398                 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
 399                 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
 400
 401             if(doTest && testSegmentCornerIntersection(seg, corner))
 402               {
 403                 double* ptA = new double[3];
 404                 copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
 405                 _polygonA.push_back(ptA);
 406                 LOG(3,"Segment-corner (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A");
 407                 if(corner != O)
 408                   {
 409                     double* ptB = new double[3];
 410                     _polygonB.push_back(ptB);
 411                     copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
 412                     LOG(3,"Segment-corner (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B");
 413                   }
 414               }
 415           }
 416
 417         // segment - ray
 418         for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
 419           {
 420             if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner))
 421               {
 422                 double* ptB = new double[3];
 423                 copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
 424                 _polygonB.push_back(ptB);
 425                 LOG(3,"Segment-ray (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B");
 426               }
 427           }
 428
 429         // segment - halfstrip
 430         for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
 431           {
 432
 433 #if 0
 434             const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX
 435             const bool doTest =
 436                 !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] &&
 437                 !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]];
 438
 439
 440             if(doTest && testSegmentHalfstripIntersection(seg, edge))
 441 #endif
 442               if(testSegmentHalfstripIntersection(seg, edge))
 443                 {
 444                   double* ptB = new double[3];
 445                   calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
 446                   _polygonB.push_back(ptB);
 447                   LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
 448                 }
 449           }
 450       }
 451
 452     // inclusion tests
 453     for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
 454       {
 455         // { XYZ - inclusion only possible if in Tetrahedron?
 456         // tetrahedron
 457         if(testCornerInTetrahedron(corner))
 458           {
 459             double* ptA = new double[3];
 460             copyVector3(&_coords[5*corner], ptA);
 461             _polygonA.push_back(ptA);
 462             LOG(3,"Inclusion tetrahedron (corner " << strTriC(corner) << "): " << vToStr(ptA) << " added to A");
 463           }
 464
 465         // XYZ - plane
 466         if(testCornerOnXYZFacet(corner))
 467           {
 468             double* ptB = new double[3];
 469             copyVector3(&_coords[5*corner], ptB);
 470             _polygonB.push_back(ptB);
 471             LOG(3,"Inclusion XYZ-plane (corner " << strTriC(corner) << "): " << vToStr(ptB) << " added to B");
 472           }
 473
 474         // projection on XYZ - facet
 475         if(testCornerAboveXYZFacet(corner))
 476           {
 477             double* ptB = new double[3];
 478             copyVector3(&_coords[5*corner], ptB);
 479             ptB[2] = 1 - ptB[0] - ptB[1];   // lower z to project on XYZ
 480             assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
 481             _polygonB.push_back(ptB);
 482             LOG(3,"Projection XYZ-plane (corner " << strTriC(corner) << "): " << vToStr(ptB) << " added to B");
 483           }
 484
 485       }
 486
 487   }
 488
 489   /**
 490    * Calculates the intersection polygon A, performing the intersection tests
 491    * and storing the corresponding point in the vector _polygonA.
 492    *
 493    * @post _polygonA contains the intersection polygon A.
 494    *
 495    */
 496   void TransformedTriangle::calculateIntersectionPolygon()
 497   {
 498     assert(_polygonA.size() == 0);
 499     // avoid reallocations in push_back() by pre-allocating enough memory
 500     // we should never have more than 20 points
 501     _polygonA.reserve(20);
 502     // -- surface intersections
 503     // surface - edge
 504     for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
 505       {
 506         if(testSurfaceEdgeIntersection(edge))
 507           {
 508             double* ptA = new double[3];
 509             calcIntersectionPtSurfaceEdge(edge, ptA);
 510             _polygonA.push_back(ptA);
 511             LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A ");
 512           }
 513       }
 514
 515     // -- segment intersections
 516     for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
 517       {
 518
 519         bool isZero[NO_DP];
 520
 521         // check beforehand which double-products are zero
 522         // Test for "== 0.0" here is OK since doubleProduct has been fixed for rounding to zero already.
 523         for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
 524           isZero[dp] = (calcStableC(seg, dp) == 0.0);
 525
 526         // segment - facet
 527         for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
 528           {
 529             // is this test worth it?
 530             const bool doTest =
 531               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
 532               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
 533               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
 534
 535             if(doTest && testSegmentFacetIntersection(seg, facet))
 536               {
 537                 double* ptA = new double[3];
 538                 calcIntersectionPtSegmentFacet(seg, facet, ptA);
 539                 _polygonA.push_back(ptA);
 540                 LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
 541               }
 542           }
 543
 544         // segment - edge
 545         for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
 546           {
 547             const DoubleProduct edge_dp = DoubleProduct(edge);
 548
 549             if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
 550               {
 551                 double* ptA = new double[3];
 552                 calcIntersectionPtSegmentEdge(seg, edge, ptA);
 553                 _polygonA.push_back(ptA);
 554                 LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
 555               }
 556           }
 557
 558         // segment - corner
 559         for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
 560           {
 561             const bool doTest =
 562               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
 563               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
 564               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
 565
 566             if(doTest && testSegmentCornerIntersection(seg, corner))
 567               {
 568                 double* ptA = new double[3];
 569                 copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
 570                 _polygonA.push_back(ptA);
 571                 LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
 572               }
 573           }
 574
 575       }
 576
 577         // inclusion tests
 578         for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
 579           {
 580             // { XYZ - inclusion only possible if in Tetrahedron?
 581             // tetrahedron
 582             if(testCornerInTetrahedron(corner))
 583               {
 584                 double* ptA = new double[3];
 585                 copyVector3(&_coords[5*corner], ptA);
 586                 _polygonA.push_back(ptA);
 587                 LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
 588               }
 589
 590           }
 591
 592       }
 593
 594
 595     /**
 596      * Returns the surface of polygon A.
 597      *
 598      * @return the surface of polygon A.
 599      */
 600     double TransformedTriangle::calculateSurfacePolygon()
 601     {
 602       const std::size_t nbPoints = _polygonA.size();
 603       double pdt[3];
 604       double sum[3] = {0., 0., 0.};
 605
 606       for(std::size_t i = 0 ; i < nbPoints ; ++i)
 607         {
 608           const double *const ptCurr = _polygonA[i];  // pt "i"
 609           const double *const ptNext = _polygonA[(i + 1) % nbPoints]; // pt "i+1" (pt nbPoints == pt 0)
 610
 611           cross(ptCurr, ptNext, pdt);
 612           add(pdt, sum);
 613         }
 614
 615       const double surface = norm(sum) * 0.5;
 616       LOG(2,"Surface is " << surface);
 617       return surface;
 618     }
 619
 620     /**
 621      * Calculates the barycenters of the given intersection polygon.
 622      *
 623      * @pre  the intersection polygons have been calculated with calculateIntersectionAndProjectionPolygons()
 624      *
 625      * @param poly        one of the two intersection polygons
 626      * @param barycenter  array of three doubles where barycenter is stored
 627      *
 628      */
 629     void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter)
 630     {
 631       LOG(3,"--- Calculating polygon barycenter");
 632
 633       // get the polygon points
 634       std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
 635
 636       // calculate barycenter
 637       const std::size_t m = polygon.size();
 638
 639       for(int j = 0 ; j < 3 ; ++j)
 640         {
 641           barycenter[j] = 0.0;
 642         }
 643
 644       if(m != 0)
 645         {
 646           for(std::size_t i = 0 ; i < m ; ++i)
 647             {
 648               const double* pt = polygon[i];
 649               for(int j = 0 ; j < 3 ; ++j)
 650                 {
 651                   barycenter[j] += pt[j] / double(m);
 652                 }
 653             }
 654         }
 655       LOG(3,"Barycenter is " << vToStr(barycenter));
 656     }
 657
 658     /**
 659      * Sorts the given intersection polygon in circular order around its barycenter.
 660      * @pre  the intersection polygons have been calculated with calculateIntersectionAndProjectionPolygons()
 661      * @post the vertices in _polygonA and _polygonB are sorted in circular order around their
 662      *       respective barycenters
 663      *
 664      * @param poly        one of the two intersection polygons
 665      * @param barycenter  array of three doubles with the coordinates of the barycenter
 666      *
 667      */
 668     void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter)
 669     {
 670       LOG(3,"--- Sorting polygon ...");
 671
 672       using INTERP_KERNEL::ProjectedCentralCircularSortOrder;
 673       typedef ProjectedCentralCircularSortOrder SortOrder; // change is only necessary here and in constructor
 674       typedef SortOrder::CoordType CoordType;
 675
 676       // get the polygon points
 677       std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
 678
 679       if(polygon.size() == 0)
 680         return;
 681
 682       // determine type of sorting
 683       CoordType type = SortOrder::XY;
 684       if(poly == A && !isTriangleInclinedToFacet( OXY )) // B is on h = 0 plane -> ok
 685         {
 686           // NB : the following test is never true if we have eliminated the
 687           // triangles parallel to x == 0 and y == 0 in calculateIntersectionVolume().
 688           // We keep the test here anyway, to avoid interdependency.
 689
 690           // is triangle inclined to x == 0 ?
 691           if(isTriangleInclinedToFacet(OZX))
 692             {
 693               type = SortOrder::XZ;
 694             }
 695           else //if(isTriangleParallelToFacet(OYZ))
 696             {
 697               type = SortOrder::YZ;
 698             }
 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 (_coords[5*P + coord] == _coords[5*Q + coord]) && (_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