src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   2 //
   3 //  This library is free software; you can redistribute it and/or
   4 //  modify it under the terms of the GNU Lesser General Public
   5 //  License as published by the Free Software Foundation; either
   6 //  version 2.1 of the License.
   7 //
   8 //  This library is distributed in the hope that it will be useful,
   9 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 //  Lesser General Public License for more details.
  12 //
  13 //  You should have received a copy of the GNU Lesser General Public
  14 //  License along with this library; if not, write to the Free Software
  15 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // File      : UnitTetraIntersectionBary.cxx
  20 // Created   : Tue Dec  9 16:48:49 2008
  21 // Author    : Edward AGAPOV (eap)
  22
  23 #include "UnitTetraIntersectionBary.hxx"
  24
  25 #include "VectorUtils.hxx"
  26 #include "InterpolationUtils.hxx"
  27 #include "VolSurfFormulae.hxx"
  28
  29 #define NB_TETRA_SIDES 4
  30 #define NB_TETRA_NODES 4
  31
  32 using namespace std;
  33
  34 namespace INTERP_KERNEL
  35 {
  36   enum { _X=0, _Y, _Z };
  37
  38   inline bool samePoint( const double* p1, const double* p2 )
  39   {
  40     return ( p1[0]==p2[0] && p1[1]==p2[1] && p1[2]==p2[2]);
  41   }
  42
  43   //================================================================================
  44   /*!
  45    * \brief Creates a ready-to-use tool
  46    */
  47   //================================================================================
  48
  49   UnitTetraIntersectionBary::UnitTetraIntersectionBary():TransformedTriangle()
  50   {
  51     _int_volume = 0;
  52     //init();
  53   }
  54   //================================================================================
  55   /*!
  56    * \brief Initializes fields
  57    */
  58   //================================================================================
  59
  60   void UnitTetraIntersectionBary::init()
  61   {
  62     _int_volume = 0;
  63     _faces.clear();
  64   }
  65
  66   //================================================================================
  67   /*!
  68    * \brief Stores a part of triangle common with the unit tetrahedron
  69    *  \param triangle - triangle side of other cell
  70    */
  71   //================================================================================
  72
  73   void UnitTetraIntersectionBary::addSide(const TransformedTriangle& triangle)
  74   {
  75     _int_volume += triangle.getVolume();
  76
  77     double triNormal[3], polyNormal[3];
  78     crossprod<3>( triangle.getCorner(P),triangle.getCorner(Q),triangle.getCorner(R), triNormal);
  79
  80     const std::vector<double*> * pPolygonA = &triangle.getPolygonA();
  81     if ( pPolygonA->size() < 3 )
  82       {
  83         if ( !epsilonEqual( triNormal[_Z], 0 ))
  84           return; // not vertical triangle does not intersect the unit tetra
  85
  86         // Vertical triangle. Use inherited methods of TransformedTriangle to
  87         // calculate intesection polygon
  88         *((TransformedTriangle*)this) = triangle; // copy triangle fields
  89         _polygonA.clear();
  90         _polygonB.clear();
  91         calculateIntersectionPolygons();
  92         if (this->_polygonA.size() < 3)
  93           return;
  94         calculatePolygonBarycenter(A, _barycenterA);
  95         sortIntersectionPolygon(A, _barycenterA);
  96         pPolygonA = & _polygonA;
  97       }
  98
  99     // check if polygon orientation is same as the one of triangle
 100     vector<double*>::const_iterator p = pPolygonA->begin(), pEnd = pPolygonA->end();
 101     double* p1 = *p++;
 102     double* p2 = *p;
 103     while ( samePoint( p1, p2 ) && ++p != pEnd )
 104       p2 = *p;
 105     if ( p == pEnd )
 106       {
 107         clearPolygons();
 108         return;
 109       }
 110     double* p3 = *p;
 111     while ( samePoint( p2, p3 ) && ++p != pEnd )
 112       p3 = *p;
 113     if ( p == pEnd )
 114       {
 115         clearPolygons();
 116         return ;
 117       }
 118     crossprod<3>( p1, p2, p3, polyNormal );
 119     bool reverse = ( dotprod<3>( triNormal, polyNormal ) < 0.0 );
 120
 121     // store polygon
 122     _faces.push_back( vector< double* > () );
 123     vector< double* >& faceCorner = _faces.back();
 124     faceCorner.resize( pPolygonA->size()/* + 1*/ );
 125
 126     int i = 0;
 127     if ( reverse )
 128       {
 129         vector<double*>::const_reverse_iterator polyF = pPolygonA->rbegin(), polyEnd;
 130         for ( polyEnd = pPolygonA->rend(); polyF != polyEnd; ++i, ++polyF )
 131           if ( i==0 || !samePoint( *polyF, faceCorner[i-1] ))
 132             copyVector3( *polyF, faceCorner[i] = new double[3] );
 133           else
 134             --i;
 135       }
 136     else
 137       {
 138         vector<double*>::const_iterator polyF = pPolygonA->begin(), polyEnd;
 139         for ( polyEnd = pPolygonA->end(); polyF != polyEnd; ++i, ++polyF )
 140           if ( i==0 || !samePoint( *polyF, faceCorner[i-1] ))
 141             copyVector3( *polyF, faceCorner[i] = new double[3] );
 142           else
 143             --i;
 144       }
 145     if ( i < 3 )
 146       {
 147         clearPolygons(); // free memory of _polygonA
 148         _polygonA = faceCorner;
 149         _faces.pop_back();
 150       }
 151     else
 152       {
 153         if ( i < (int)pPolygonA->size() )
 154           faceCorner.resize( i );
 155       }
 156
 157 #ifdef DMP_ADDSIDE
 158     cout << "**** addSide() " << _faces.size() << endl;
 159     for ( int i = 0; i < faceCorner.size(); ++i )
 160       {
 161         double* p = faceCorner[i];
 162         cout << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl;
 163       }
 164 #endif
 165     clearPolygons();
 166   }
 167
 168   //================================================================================
 169   /*!
 170    * \brief Computes and returns coordinates of barycentre
 171    */
 172   //================================================================================
 173
 174   bool UnitTetraIntersectionBary::getBary(double* baryCenter)
 175   {
 176     baryCenter[0] = baryCenter[1] = baryCenter[2] = -1.0;
 177     if ( addSideFaces() < NB_TETRA_SIDES )
 178       {
 179         // tetra is not intersected
 180         if ( _int_volume != 0.0 )
 181           {
 182             // tetra is fully inside the other cell
 183             baryCenter[0] = baryCenter[1] = baryCenter[2] = 0.25;
 184             _int_volume = 0.16666666666666666;
 185             return true;
 186           }
 187         return false;
 188       }
 189     // Algo:
 190     // - pick up one point P among the summits of the polyhedron
 191     // - for each face of the polyhedron which does not contain the point :
 192     //   - compute the barycenter of the volume obtained by forming the "pyramid" with
 193     //     the face as a base and point P as a summit
 194     //   - compute the volume of the "pyramid"
 195     // - Add up all barycenter positions weighting them with the volumes.
 196
 197     baryCenter[0] = baryCenter[1] = baryCenter[2] = 0.;
 198
 199     list< vector< double* > >::iterator f = _faces.begin(), fEnd = _faces.end();
 200     double * P = f->at(0);
 201
 202     for ( ++f; f != fEnd; ++f )
 203       {
 204         vector< double* >& polygon = *f;
 205         if ( polygon.empty() )
 206           continue;
 207
 208         bool pBelongsToPoly = false;
 209         vector<double*>::iterator v = polygon.begin(), vEnd = polygon.end();
 210         for ( ; !pBelongsToPoly && v != vEnd; ++v )
 211           pBelongsToPoly = samePoint( P, *v );
 212         if ( pBelongsToPoly )
 213           continue;
 214
 215         // Compute the barycenter of the volume. Barycenter of pyramid is on line
 216         // ( barycenter of polygon -> P ) with 1/4 of pyramid height from polygon.
 217
 218         double bary[] = { 0, 0, 0 };
 219
 220         // base polygon bary
 221         for ( v = polygon.begin(); v != vEnd ; ++v )
 222           {
 223             double* p = *v;
 224             bary[0] += p[0];
 225             bary[1] += p[1];
 226             bary[2] += p[2];
 227           }
 228         bary[0] /= polygon.size();
 229         bary[1] /= polygon.size();
 230         bary[2] /= polygon.size();
 231
 232         // pyramid volume
 233         double vol = 0;
 234         for ( int i = 0; i < (int)polygon.size(); ++i )
 235           {
 236             double* p1 = polygon[i];
 237             double* p2 = polygon[(i+1)%polygon.size()];
 238             vol += std::fabs( calculateVolumeForTetra( p1, p2, bary, P ));
 239           }
 240
 241         // put bary on the line ( barycenter of polygon -> P ) and multiply by volume
 242         baryCenter[0] += ( bary[0] * 0.75 + P[0] * 0.25 ) * vol;
 243         baryCenter[1] += ( bary[1] * 0.75 + P[1] * 0.25 ) * vol;
 244         baryCenter[2] += ( bary[2] * 0.75 + P[2] * 0.25 ) * vol;
 245       }
 246     if ( _int_volume < 0. )
 247       _int_volume = -_int_volume;
 248     baryCenter[0] /= _int_volume;
 249     baryCenter[1] /= _int_volume;
 250     baryCenter[2] /= _int_volume;
 251
 252     return true;
 253   }
 254
 255   //================================================================================
 256   /*!
 257    * \brief Add faces of the intersection polyhedron formed on faces of the
 258    *        unit tetrahedron by sides of already added faces
 259    *  \retval int - number of faces of intersection polyhedron
 260    */
 261   //================================================================================
 262
 263   int UnitTetraIntersectionBary::addSideFaces()
 264   {
 265     int nbPolyhedraFaces = 0;
 266
 267     if ( _faces.empty() )
 268       return nbPolyhedraFaces;
 269
 270     // -------------------------------------------
 271     // Detect polygons laying on sides of a tetra
 272     // -------------------------------------------
 273
 274     bool sideAdded[NB_TETRA_SIDES] = { false, false, false, false };
 275     int nbAddedSides = 0;
 276     list< vector< double* > >::iterator f = _faces.begin(), fEnd = _faces.end();
 277     for ( ; f != fEnd; ++f )
 278       {
 279         vector< double* >& polygon = *f;
 280         double coordSum[3] = {0,0,0};
 281         for ( int i = 0; i < (int)polygon.size(); ++i )
 282           {
 283             double* p = polygon[i];
 284             coordSum[0] += p[0];
 285             coordSum[1] += p[1];
 286             coordSum[2] += p[2];
 287           }
 288         for ( int j = 0; j < 3 && !sideAdded[j]; ++j )
 289           {
 290             if ( coordSum[j] == 0 )
 291               sideAdded[j] = bool( ++nbAddedSides );
 292           }
 293         if ( !sideAdded[3] &&
 294              ( epsilonEqual( (coordSum[0]+coordSum[1]+coordSum[2]) / polygon.size(), 1. )))
 295           sideAdded[3] = bool( ++nbAddedSides );
 296       }
 297     if ( nbAddedSides == NB_TETRA_SIDES )
 298       return nbAddedSides;
 299
 300     // ---------------------------------------------------------------------------------
 301     // Add segments of already added polygons to future polygonal faces on sides of tetra
 302     // ---------------------------------------------------------------------------------
 303
 304     int nbIntersectPolygs = _faces.size();
 305
 306     vector< double* > * sideFaces[ 4 ]; // future polygons on sides of tetra
 307     for ( int i = 0; i < NB_TETRA_SIDES; ++i )
 308       {
 309         sideFaces[ i ]=0;
 310         if ( !sideAdded[ i ] )
 311           {
 312             _faces.push_back( vector< double* > () );
 313             sideFaces[ i ] = &_faces.back();
 314           }
 315       }
 316     f = _faces.begin(), fEnd = _faces.end();
 317     for ( int iF = 0; iF < nbIntersectPolygs; ++f, ++iF ) // loop on added intersection polygons
 318       {
 319         vector< double* >& polygon = *f;
 320         for ( int i = 0; i < (int)polygon.size(); ++i )
 321           {
 322             // segment ends
 323             double* p1 = polygon[i];
 324             double* p2 = polygon[(i+1)%polygon.size()];
 325             bool p1OnSide, p2OnSide;//, onZeroSide = false;
 326             for ( int j = 0; j < 3; ++j )
 327               {
 328                 if ( !sideFaces[ j ] )
 329                   continue;
 330                 p1OnSide = ( p1[j] == 0 );
 331                 p2OnSide = ( p2[j] == 0 );
 332                 if ( p1OnSide && p2OnSide )
 333                   {
 334                     // segment p1-p2 is on j-th orthogonal side of tetra
 335                     sideFaces[j]->push_back( new double[3] );
 336                     copyVector3( p1, sideFaces[j]->back() );
 337                     sideFaces[j]->push_back( new double[3] );
 338                     copyVector3( p2, sideFaces[j]->back() );
 339                     //break;
 340                   }
 341               }
 342             // check if the segment p1-p2 is on the inclined side
 343             if ( sideFaces[3] &&
 344                  epsilonEqual( p1[_X] + p1[_Y] + p1[_Z], 1.0 ) &&
 345                  epsilonEqual( p2[_X] + p2[_Y] + p2[_Z], 1.0 ))
 346               {
 347                 sideFaces[3]->push_back( new double[3] );
 348                 copyVector3( p1, sideFaces[3]->back() );
 349                 sideFaces[3]->push_back( new double[3] );
 350                 copyVector3( p2, sideFaces[3]->back() );
 351               }
 352           }
 353       }
 354 #ifdef DMP_ADDSIDEFACES
 355     cout << "**** after Add segments to sides " << endl;
 356     for ( int i = 0; i < NB_TETRA_SIDES; ++i )
 357       {
 358         cout << "\t Side " << i << endl;
 359         if ( !sideFaces[i] )
 360           {
 361             cout << "\t cut by triagle" << endl;
 362           }
 363         else
 364           {
 365             vector< double* >& sideFace = *sideFaces[i];
 366             for ( int i = 0; i < sideFace.size(); ++i )
 367               {
 368                 double* p = sideFace[i];
 369                 cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl;
 370               }
 371           }
 372       }
 373 #endif
 374
 375     // ---------------------------------------------------------------------------
 376     // Make closed polygons on tetra sides by adding not cut off corners of tetra
 377     // ---------------------------------------------------------------------------
 378
 379     double origin[3] = { 0,0,0 };
 380
 381     // find corners of tetra cut off by triangles of other tetra
 382     // ---------------------------------------------------------
 383
 384     // corners are coded like this: index = 1*X + 2*Y + 3*Z
 385     // (0,0,0) -> index == 0; (0,0,1) -> index == 3
 386     vector< int > cutOffCorners(NB_TETRA_NODES, false), passedCorners(NB_TETRA_NODES, false);
 387
 388     int nbCutOffCorners = 0;
 389     for ( int i = 0; i < 3; ++i ) // loop on orthogonal faces of the unit tetra
 390       {
 391         if ( !sideFaces[i] ) continue;
 392         vector< double* >& sideFace = *sideFaces[i];
 393
 394         int nbPoints = sideFace.size();
 395         if ( nbPoints == 0 )
 396           continue; // not intersected face at all - no cut off corners can be detected
 397
 398         int ind1 = (i+1)%3, ind2 = (i+2)%3; // indices of coords on i-th tetra side
 399
 400         int nbCutOnSide = 0;
 401         bool isSegmentOnEdge;
 402         for ( int ip = 0; ip < nbPoints; ++ip )
 403           {
 404             int isSegmentEnd = ( ip % 2 );
 405
 406             double* p = sideFace[ ip ];
 407             double* p2 = isSegmentEnd ? 0 : sideFace[ip+1];
 408
 409             if ( !isSegmentEnd )
 410               isSegmentOnEdge = false; // initialize
 411
 412             int cutOffIndex = -1, passThIndex = -1;// no cut off neither pass through
 413             int pCut[] = { 0,0,0 }, pPass[] = { 0,0,0 };
 414
 415             if ( p[ind1]==0.)
 416               {
 417                 // point is in on orthogonal edge
 418                 if ( !isSegmentEnd && p2[ind1]==0 )
 419                   isSegmentOnEdge = true;
 420
 421                 if ( !isSegmentOnEdge )
 422                   { // segment ends are on different edges
 423                     pCut[ind2] = isSegmentEnd; // believe that cutting triangles are well oriented
 424                     cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2];
 425                   }
 426                 if ( p[ind2]==0. || p[ind2]==1.)
 427                   {
 428                     pPass[ind2] = int(p[ind2]);
 429                     passThIndex = pPass[0] + 2*pPass[1] + 3*pPass[2];
 430                   }
 431               }
 432             else if ( p[ind2]==0.)
 433               {
 434                 // point is on orthogonal edge
 435                 if ( !isSegmentEnd && p2[ind2]==0 )
 436                   isSegmentOnEdge = true;
 437                 if ( !isSegmentEnd )
 438                   {// segment ends are on different edges
 439                     pCut[ind1] = 1-isSegmentEnd;
 440                     cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2];
 441                   }
 442                 if ( p[ind1]==0. || p[ind1]==1.)
 443                   {
 444                     pPass[ind1] = int(p[ind1]);
 445                     passThIndex = pPass[0] + 2*pPass[1] + 3*pPass[2];
 446                   }
 447               }
 448             else if ( epsilonEqual(p[ind1] + p[ind2], 1.0 ))
 449               {
 450                 // point is on inclined edge
 451                 if ( !isSegmentEnd && epsilonEqual(p2[ind1] + p2[ind2], 1.0 ))
 452                   isSegmentOnEdge = true;
 453                 if ( !isSegmentOnEdge )
 454                   { //segment ends are on different edges
 455                     pCut[ind1] = isSegmentEnd;
 456                     pCut[ind2] = 1-isSegmentEnd;
 457                     cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2];
 458                   }
 459               }
 460             else
 461               {
 462                 continue;
 463               }
 464             // remember cut off and passed through points
 465             if ( passThIndex >= 0. )
 466               {
 467                 passedCorners[ passThIndex ] = true;
 468                 if ( cutOffCorners[ passThIndex ] )
 469                   {
 470                     nbCutOffCorners--;
 471                     cutOffCorners[ passThIndex ] = false;
 472                   }
 473               }
 474             if ( cutOffIndex >= 0. )
 475               {
 476                 nbCutOnSide++;
 477                 if ( !passedCorners[ cutOffIndex ] && !cutOffCorners[ cutOffIndex ] )
 478                   {
 479                     nbCutOffCorners++;
 480                     cutOffCorners[ cutOffIndex ] = true;
 481                   }
 482               }
 483           } // loop on points on a unit tetra side
 484
 485         if ( nbCutOnSide == 0 && nbPoints <= 2 )
 486           continue; // one segment laying on edge at most
 487
 488         if ( nbCutOffCorners == NB_TETRA_NODES )
 489           break; // all tetra corners are cut off
 490
 491         if ( /*nbCutOnSide <= 2 &&*/ nbPoints >= 6 )
 492           {
 493             // at least 3 segments - all corners of a side are cut off
 494             for (int cutIndex = 0; cutIndex < NB_TETRA_NODES; ++cutIndex )
 495               if ( cutIndex != i+1 && !passedCorners[ cutIndex ] && !cutOffCorners[ cutIndex ])
 496                 cutOffCorners[ cutIndex ] = bool ( ++nbCutOffCorners );
 497           }
 498
 499       } // loop on orthogonal faces of tetra
 500
 501     // check if all corners are cut off on the inclined tetra side
 502     if ( sideFaces[ XYZ ] && sideFaces[ XYZ ]->size() >= 6 )
 503       {
 504         for (int cutIndex = 1; cutIndex < NB_TETRA_NODES; ++cutIndex )
 505           if ( !passedCorners[ cutIndex ] && !cutOffCorners[ cutIndex ])
 506             cutOffCorners[ cutIndex ] = bool ( ++nbCutOffCorners );
 507       }
 508
 509     // Add to faces on tetra sides the corners not cut off by segments of intersection polygons
 510     // ----------------------------------------------------------------------------------
 511     if ( nbCutOffCorners > 0 )
 512       {
 513         for ( int i = 0; i < NB_TETRA_SIDES; ++i )
 514           {
 515             if ( !sideFaces[ i ] ) continue;
 516             vector< double* >& sideFace = *sideFaces[i];
 517
 518             int excludeCorner = (i + 1) % NB_TETRA_NODES;
 519             for ( int ic = 0; ic < NB_TETRA_NODES; ++ic )
 520               {
 521                 if ( !cutOffCorners[ ic ] && ic != excludeCorner )
 522                   {
 523                     sideFace.push_back( new double[3] );
 524                     copyVector3( origin, sideFace.back() );
 525                     if ( ic )
 526                       sideFace.back()[ ic-1 ] = 1.0;
 527                   }
 528               }
 529           }
 530       }
 531
 532 #ifdef DMP_ADDSIDEFACES
 533     cout << "**** after Add corners to sides " << endl;
 534     for ( int i = 0; i < NB_TETRA_SIDES; ++i )
 535       {
 536         cout << "\t Side " << i << endl;
 537         if ( !sideFaces[i] ) {
 538           cout << "\t cut by triagle" << endl;
 539         }
 540         else
 541           {
 542             vector< double* >& sideFace = *sideFaces[i];
 543             for ( int i = 0; i < sideFace.size(); ++i )
 544               {
 545                 double* p = sideFace[i];
 546                 cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl;
 547               }
 548           }
 549       }
 550     cout << "Cut off corners: ";
 551     if ( nbCutOffCorners == 0 )
 552       cout << "NO";
 553     else
 554       for ( int ic = 0; ic < NB_TETRA_NODES; ++ic )
 555         cout << cutOffCorners[ ic ];
 556     cout << endl;
 557 #endif
 558     // ------------------------------------------------------------------------
 559     // Sort corners of filled up faces on tetra sides and exclude equal points
 560     // ------------------------------------------------------------------------
 561
 562     int iF = 0;
 563     for ( f = _faces.begin(); f != fEnd; ++f, ++iF )
 564       {
 565         vector< double* >&  face = *f;
 566         if ( face.size() >= 3 )
 567           {
 568             clearPolygons(); // free memory of _polygonA
 569             _polygonA = face;
 570             face.clear();
 571             face.reserve( _polygonA.size() );
 572             if ( iF >= nbIntersectPolygs )
 573               { // sort points of side faces
 574                 calculatePolygonBarycenter( A, _barycenterA );
 575                 setTriangleOnSide( iF - nbIntersectPolygs );
 576                 sortIntersectionPolygon( A, _barycenterA );
 577               }
 578             // exclude equal points
 579             vector< double* >::iterator v = _polygonA.begin(), vEnd = _polygonA.end();
 580             face.push_back( *v );
 581             *v = 0;
 582             for ( ++v; v != vEnd; ++v )
 583               {
 584                 double* pPrev = face.back();
 585                 double* p     = *v;
 586                 if ( !samePoint( p, pPrev ))
 587                   {
 588                     face.push_back( p );
 589                     *v = 0;
 590                   }
 591               }
 592           }
 593         if ( face.size() < 3 )
 594           { // size could decrease
 595             clearPolygons(); // free memory of _polygonA
 596             _polygonA = face;
 597             face.clear();
 598           }
 599         else
 600           {
 601             nbPolyhedraFaces++;
 602           }
 603       }
 604 #ifdef DMP_ADDSIDEFACES
 605     cout << "**** after HEALING all faces " << endl;
 606     for (iF=0, f = _faces.begin(); f != fEnd; ++f, ++iF )
 607       {
 608         cout << "\t Side " << iF << endl;
 609         vector< double* >& sideFace = *f;
 610         for ( int i = 0; i < sideFace.size(); ++i )
 611           {
 612             double* p = sideFace[i];
 613             cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl;
 614           }
 615       }
 616 #endif
 617     return nbPolyhedraFaces;
 618   }
 619
 620   //================================================================================
 621   /*!
 622    * \brief set corners of inherited TransformedTriangle as corners of i-th side of
 623    * the Unit tetra. It is necessary to sort points of faces on sides of the unit
 624    * tetrahedron using sortIntersectionPolygon(A)
 625    */
 626   //================================================================================
 627
 628   void UnitTetraIntersectionBary::setTriangleOnSide(int iSide)
 629   {
 630     if ( iSide >= 3 )
 631       iSide = 0;
 632     for(int i = 0 ; i < 3 ; ++i)
 633       {
 634         _coords[5*i] = _coords[5*i + 1] = _coords[5*i + 2] = 0.;
 635         if ( i != iSide )
 636           _coords[5*i + i] = 1.;
 637       }
 638   }
 639
 640   //================================================================================
 641   /*!
 642    * \brief Free memory of polygons
 643    */
 644   //================================================================================
 645
 646   void UnitTetraIntersectionBary::clearPolygons(bool andFaces)
 647   {
 648     for(vector<double*>::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it)
 649       {  delete[] *it;
 650         *it = 0;
 651       }
 652     for(vector<double*>::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it)
 653       {
 654         delete[] *it;
 655         *it = 0;
 656       }
 657
 658     _polygonA.clear();
 659     _polygonB.clear();
 660
 661     if ( andFaces )
 662       {
 663         list< vector< double* > >::iterator f = this->_faces.begin(), fEnd = this->_faces.end();
 664         for ( ; f != fEnd; ++f )
 665           {
 666             vector< double* >& polygon = *f;
 667             for(vector<double*>::iterator it = polygon.begin() ; it != polygon.end() ; ++it)
 668               {
 669                 delete[] *it;
 670                 *it = 0;
 671               }
 672           }
 673         this->_faces.clear();
 674       }
 675   }
 676
 677   //================================================================================
 678   /*!
 679    * \brief Destructor clears coordinates of faces
 680    */
 681   //================================================================================
 682
 683   UnitTetraIntersectionBary::~UnitTetraIntersectionBary()
 684   {
 685     clearPolygons(/*andFaces=*/true );
 686   }
 687
 688 }