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