src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx

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