Back to index

salome-med  6.5.0
UnitTetraIntersectionBary.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 // File      : UnitTetraIntersectionBary.cxx
00021 // Created   : Tue Dec  9 16:48:49 2008
00022 // Author    : Edward AGAPOV (eap)
00023 //
00024 #include "UnitTetraIntersectionBary.hxx"
00025 
00026 #include "VectorUtils.hxx"
00027 #include "InterpolationUtils.hxx"
00028 #include "VolSurfFormulae.hxx"
00029 
00030 #define NB_TETRA_SIDES 4
00031 #define NB_TETRA_NODES 4
00032 
00033 //#define DMP_UNITTETRAINTERSECTIONBARY
00034 
00035 
00036 namespace INTERP_KERNEL
00037 {
00038   enum { _X=0, _Y, _Z };
00039 
00040   inline bool samePoint( const double* p1, const double* p2 )
00041   {
00042     return ( epsilonEqual( p1[0], p2[0]) &&
00043              epsilonEqual( p1[1], p2[1]) &&
00044              epsilonEqual( p1[2], p2[2]));
00045   }
00046 
00047   //================================================================================
00051   //================================================================================
00052 
00053   UnitTetraIntersectionBary::UnitTetraIntersectionBary(bool isTetraInversed)
00054     :TransformedTriangle(),_int_volume(0),_isTetraInversed( isTetraInversed )
00055   {
00056     //init();
00057   }
00058   //================================================================================
00062   //================================================================================
00063 
00064   void UnitTetraIntersectionBary::init(bool isTetraInversed)
00065   {
00066     _int_volume = 0;
00067     _isTetraInversed = isTetraInversed;
00068     _faces.clear();
00069     _polyNormals.clear();
00070   }
00071   
00072   //================================================================================
00077   //================================================================================
00078 
00079   void UnitTetraIntersectionBary::addSide(const TransformedTriangle& triangle)
00080   {
00081     _int_volume += triangle.getVolume();
00082 
00083     double triNormal[3], polyNormal[3];
00084     crossprod<3>( triangle.getCorner(P),triangle.getCorner(Q),triangle.getCorner(R), triNormal);
00085 
00086     const std::vector<double*> * pPolygonA = &triangle.getPolygonA();
00087     if ( pPolygonA->size() < 3 )
00088       {
00089         if ( !epsilonEqual( triNormal[_Z], 0 ))
00090           return; // not vertical triangle does not intersect the unit tetra
00091 
00092         // Vertical triangle. Use inherited methods of TransformedTriangle to
00093         // calculate intesection polygon
00094         *((TransformedTriangle*)this) = triangle; // copy triangle fields
00095         _polygonA.clear();
00096         _polygonB.clear();
00097         calculateIntersectionAndProjectionPolygons();
00098         if (this->_polygonA.size() < 3)
00099           return;
00100         calculatePolygonBarycenter(A, _barycenterA);
00101         sortIntersectionPolygon(A, _barycenterA);
00102         pPolygonA = & _polygonA;
00103       }
00104 
00105     // check if polygon orientation is same as the one of triangle
00106     std::vector<double*>::const_iterator p = pPolygonA->begin(), pEnd = pPolygonA->end();
00107 #ifdef DMP_UNITTETRAINTERSECTIONBARY
00108     std::cout.precision(18);
00109     std::cout << "**** int polygon() " << std::endl;
00110     while ( p != pEnd )
00111     {
00112       double* pp = *p++;
00113       std::cout << pEnd - p << ": ( " << pp[0] << ", " << pp[1] << ", " << pp[2] << " )" << std::endl;
00114     }
00115     p = pPolygonA->begin();
00116 #endif
00117     double* p1 = *p++;
00118     double* p2 = *p;
00119     while ( samePoint( p1, p2 ) && ++p != pEnd )
00120       p2 = *p;
00121     if ( p == pEnd )
00122       {
00123 #ifdef DMP_UNITTETRAINTERSECTIONBARY
00124         std::cout << "All points equal" << std::endl;
00125 #endif
00126         clearPolygons();
00127         return;
00128       }
00129     double* p3 = *p;
00130     while (( samePoint( p2, p3 ) || samePoint( p1, p3 )) && ++p != pEnd )
00131       p3 = *p;
00132     if ( p == pEnd )
00133       {
00134 #ifdef DMP_UNITTETRAINTERSECTIONBARY
00135         std::cout << "Only two points differ" << std::endl;
00136 #endif
00137         clearPolygons();
00138         return ;
00139       }
00140     crossprod<3>( p1, p2, p3, polyNormal );
00141     bool reverse = ( dotprod<3>( triNormal, polyNormal ) < 0.0 );
00142     if (_isTetraInversed) reverse = !reverse;
00143 
00144     // store polygon
00145     _faces.push_back( std::vector< double* > () );
00146     std::vector< double* >& faceCorner = _faces.back();
00147     faceCorner.resize( pPolygonA->size()/* + 1*/ );
00148 
00149     int i = 0;
00150     if ( reverse )
00151       {
00152         std::vector<double*>::const_reverse_iterator polyF = pPolygonA->rbegin(), polyEnd;
00153         for ( polyEnd = pPolygonA->rend(); polyF != polyEnd; ++i, ++polyF )
00154           if ( i==0 || !samePoint( *polyF, faceCorner[i-1] ))
00155             copyVector3( *polyF, faceCorner[i] = new double[3] );
00156           else
00157             --i;
00158         polyNormal[0] *= -1.;
00159         polyNormal[1] *= -1.;
00160         polyNormal[2] *= -1.;
00161       }
00162     else
00163       {
00164         std::vector<double*>::const_iterator polyF = pPolygonA->begin(), polyEnd;
00165         for ( polyEnd = pPolygonA->end(); polyF != polyEnd; ++i, ++polyF )
00166           if ( i==0 || !samePoint( *polyF, faceCorner[i-1] ))
00167             copyVector3( *polyF, faceCorner[i] = new double[3] );
00168           else
00169             --i;
00170       }
00171     if ( i < 3 )
00172       {
00173         clearPolygons(); // free memory of _polygonA
00174         _polygonA = faceCorner;
00175         _faces.pop_back();
00176       }
00177     else
00178       {
00179         if ( i < (int)pPolygonA->size() )
00180           faceCorner.resize( i );
00181 
00182         if ( _polyNormals.empty() )
00183           _polyNormals.reserve(4);
00184         _polyNormals.push_back( std::vector< double >( polyNormal, polyNormal+3 ));
00185       }
00186 
00187 #ifdef DMP_UNITTETRAINTERSECTIONBARY
00188     std::cout << "**** addSide() " << _faces.size() << std::endl;
00189     for ( int i = 0; i < faceCorner.size(); ++i )
00190       {
00191         double* p = faceCorner[i];
00192         std::cout << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << std::endl;
00193       }
00194     std::cout << "NORM: ( " << _polyNormals.back()[0] << ", " << _polyNormals.back()[1] << ", " << _polyNormals.back()[2] << " )" << std::endl;
00195 #endif
00196     clearPolygons();
00197   }
00198 
00199   //================================================================================
00203   //================================================================================
00204 
00205   bool UnitTetraIntersectionBary::getBary(double* baryCenter)
00206   {
00207     baryCenter[0] = baryCenter[1] = baryCenter[2] = -1.0;
00208     if ( addSideFaces() < NB_TETRA_SIDES )
00209       {
00210         // tetra is not intersected
00211         if ( fabs(_int_volume) > 1e-10 )
00212           {
00213             // tetra is fully inside the other cell
00214             baryCenter[0] = baryCenter[1] = baryCenter[2] = 0.25;
00215             _int_volume = 0.16666666666666666;
00216             return true;
00217           }
00218         return false;
00219       }
00220     // Algo:
00221     // - pick up one point P among the summits of the polyhedron
00222     // - for each face of the polyhedron which does not contain the point :
00223     //   - compute the barycenter of the volume obtained by forming the "pyramid" with
00224     //     the face as a base and point P as a summit
00225     //   - compute the volume of the "pyramid"
00226     // - Add up all barycenter positions weighting them with the volumes.
00227 
00228     baryCenter[0] = baryCenter[1] = baryCenter[2] = 0.;
00229 
00230     std::list< std::vector< double* > >::iterator f = _faces.begin(), fEnd = _faces.end();
00231     double * PP = f->at(0);
00232 
00233     for ( ++f; f != fEnd; ++f )
00234       {
00235         std::vector< double* >& polygon = *f;
00236         if ( polygon.empty() )
00237           continue;
00238 
00239         bool pBelongsToPoly = false;
00240         std::vector<double*>::iterator v = polygon.begin(), vEnd = polygon.end();
00241         for ( ; !pBelongsToPoly && v != vEnd; ++v )
00242           pBelongsToPoly = samePoint( PP, *v );
00243         if ( pBelongsToPoly )
00244           continue;
00245 
00246         // Compute the barycenter of the volume. Barycenter of pyramid is on line
00247         // ( barycenter of polygon -> PP ) with 1/4 of pyramid height from polygon.
00248 
00249         double bary[] = { 0, 0, 0 };
00250 
00251         // base polygon bary
00252         for ( v = polygon.begin(); v != vEnd ; ++v )
00253           {
00254             double* p = *v;
00255             bary[0] += p[0];
00256             bary[1] += p[1];
00257             bary[2] += p[2];
00258           }
00259         bary[0] /= (int)polygon.size();
00260         bary[1] /= (int)polygon.size();
00261         bary[2] /= (int)polygon.size();
00262 
00263         // pyramid volume
00264         double vol = 0;
00265         for ( int i = 0; i < (int)polygon.size(); ++i )
00266           {
00267             double* p1 = polygon[i];
00268             double* p2 = polygon[(i+1)%polygon.size()];
00269             vol += std::fabs( calculateVolumeForTetra( p1, p2, bary, PP ));
00270           }
00271 
00272         // put bary on the line ( barycenter of polygon -> PP ) and multiply by volume
00273         baryCenter[0] += ( bary[0] * 0.75 + PP[0] * 0.25 ) * vol;
00274         baryCenter[1] += ( bary[1] * 0.75 + PP[1] * 0.25 ) * vol;
00275         baryCenter[2] += ( bary[2] * 0.75 + PP[2] * 0.25 ) * vol;
00276       }
00277     if ( _int_volume < 0. )
00278       _int_volume = -_int_volume;
00279     baryCenter[0] /= _int_volume;
00280     baryCenter[1] /= _int_volume;
00281     baryCenter[2] /= _int_volume;
00282 
00283 #ifdef DMP_UNITTETRAINTERSECTIONBARY
00284     std::cout.precision(5);
00285     std::cout << "**** Barycenter " << baryCenter[0] <<", "<< baryCenter[1] <<", "<< baryCenter[2]
00286          << "\t **** Volume " << _int_volume << std::endl;
00287     std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
00288 #endif
00289     return true;
00290   }
00291 
00292   //================================================================================
00298   //================================================================================
00299 
00300   int UnitTetraIntersectionBary::addSideFaces()
00301   {
00302     int nbPolyhedraFaces = 0;
00303 
00304     if ( _faces.empty() )
00305       return nbPolyhedraFaces;
00306 
00307     // -------------------------------------------
00308     // Detect polygons laying on sides of a tetra
00309     // -------------------------------------------
00310 
00311     bool sideAdded[NB_TETRA_SIDES] = { false, false, false, false };
00312     int nbAddedSides = 0;
00313     std::list< std::vector< double* > >::iterator f = _faces.begin(), fEnd = _faces.end();
00314     for ( ; f != fEnd; ++f )
00315     {
00316       std::vector< double* >& polygon = *f;
00317       double coordSum[3] = {0,0,0};
00318       for ( int i = 0; i < (int)polygon.size(); ++i )
00319       {
00320         double* p = polygon[i];
00321         coordSum[0] += p[0];
00322         coordSum[1] += p[1];
00323         coordSum[2] += p[2];
00324       }
00325       for ( int j = 0; j < 3 && !sideAdded[j]; ++j )
00326       {
00327         if ( epsilonEqual( coordSum[j], 0.0 ))
00328           sideAdded[j] = ++nbAddedSides != 0 ;
00329       }
00330       if ( !sideAdded[3] &&
00331            ( epsilonEqual( (coordSum[0]+coordSum[1]+coordSum[2]) / (int)polygon.size(), 1. )))
00332         sideAdded[3] = ++nbAddedSides != 0 ;
00333     }
00334     if ( nbAddedSides == NB_TETRA_SIDES )
00335       return nbAddedSides;
00336 
00337     // ---------------------------------------------------------------------------------
00338     // Add segments of already added polygons to future polygonal faces on sides of tetra
00339     // ---------------------------------------------------------------------------------
00340 
00341     std::size_t nbIntersectPolygs = _faces.size();
00342 
00343     std::vector< double* > * sideFaces[ 4 ]; // future polygons on sides of tetra
00344     for ( int i = 0; i < NB_TETRA_SIDES; ++i )
00345     {
00346       sideFaces[ i ]=0;
00347       if ( !sideAdded[ i ] )
00348       {
00349         _faces.push_back( std::vector< double* > () );
00350         sideFaces[ i ] = &_faces.back();
00351       }
00352     }
00353     f = _faces.begin(), fEnd = _faces.end();
00354     for ( std::size_t iF = 0; iF < nbIntersectPolygs; ++f, ++iF ) // loop on added intersection polygons
00355     {
00356       std::vector< double* >& polygon = *f;
00357       for ( std::size_t i = 0; i < polygon.size(); ++i )
00358       {
00359         // segment ends
00360         double* p1 = polygon[i];
00361         double* p2 = polygon[(i+1)%polygon.size()];
00362         bool p1OnSide, p2OnSide;//, onZeroSide = false;
00363         for ( int j = 0; j < 3; ++j )
00364         {
00365           if ( !sideFaces[ j ] )
00366             continue;
00367           p1OnSide = epsilonEqual( p1[j], 0. );
00368           p2OnSide = epsilonEqual( p2[j], 0. );
00369           if ( p1OnSide && p2OnSide )
00370           {
00371             // segment p1-p2 is on j-th orthogonal side of tetra
00372             sideFaces[j]->push_back( new double[3] );
00373             copyVector3( p1, sideFaces[j]->back() );
00374             sideFaces[j]->push_back( new double[3] );
00375             copyVector3( p2, sideFaces[j]->back() );
00376             //break;
00377           }
00378         }
00379         // check if the segment p1-p2 is on the inclined side
00380         if ( sideFaces[3] &&
00381              epsilonEqual( p1[_X] + p1[_Y] + p1[_Z], 1.0 ) &&
00382              epsilonEqual( p2[_X] + p2[_Y] + p2[_Z], 1.0 ))
00383         {
00384           sideFaces[3]->push_back( new double[3] );
00385           copyVector3( p1, sideFaces[3]->back() );
00386           sideFaces[3]->push_back( new double[3] );
00387           copyVector3( p2, sideFaces[3]->back() );
00388         }
00389       }
00390     }
00391 #ifdef DMP_UNITTETRAINTERSECTIONBARY
00392     std::cout << "**** after Add segments to sides " << std::endl;
00393     for ( int i = 0; i < NB_TETRA_SIDES; ++i )
00394     {
00395       std::cout << "\t Side " << i << std::endl;
00396       if ( !sideFaces[i] )
00397       {
00398         std::cout << "\t cut by triagle" << std::endl;
00399       }
00400       else
00401       {
00402         std::vector< double* >& sideFace = *sideFaces[i];
00403         for ( int i = 0; i < sideFace.size(); ++i )
00404         {
00405           double* p = sideFace[i];
00406           std::cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << std::endl;
00407         }
00408       }
00409     }
00410 #endif
00411 
00412     // ---------------------------------------------------------------------------
00413     // Make closed polygons on tetra sides by adding not cut off corners of tetra
00414     // ---------------------------------------------------------------------------
00415 
00416     double origin[3] = { 0,0,0 };
00417 
00418     // find corners of tetra cut off by triangles of other tetra
00419     // ---------------------------------------------------------
00420 
00421     // corners are coded like this: index = 1*X + 2*Y + 3*Z
00422     // (0,0,0) -> index == 0; (0,0,1) -> index == 3
00423     int cutOffCorners[NB_TETRA_NODES] = { false, false, false, false };
00424     int passedCorners[NB_TETRA_NODES] = { false, false, false, false };
00425 
00426     // find cutOffCorners by normals of intersection polygons
00427     int nbCutOffCorners = 0;
00428     for ( int ic = 0; ic < NB_TETRA_NODES; ++ic )
00429     {
00430       f = _faces.begin(), fEnd = _faces.end();
00431       for ( std::size_t iF = 0; iF < nbIntersectPolygs; ++f, ++iF ) // loop on added intersection polygons
00432       {
00433         std::vector< double* >& polygon = *f;
00434 
00435         double corner2Poly[3] = { polygon[0][0], polygon[0][1], polygon[0][2] };
00436         if ( ic ) corner2Poly[ ic-1 ] -= 1.0;
00437 
00438         // _polyNormals are outside of a tetrahedron
00439         double dot = dotprod<3>( corner2Poly, &_polyNormals[iF][0] );
00440         if ( dot < -DEFAULT_ABS_TOL*DEFAULT_ABS_TOL )
00441         {
00442 #ifdef DMP_UNITTETRAINTERSECTIONBARY
00443           std::cout << "side " << iF+1 << ": cut " << ic << std::endl;
00444 #endif
00445           cutOffCorners[ ic ] = true;
00446           nbCutOffCorners++;
00447           break;
00448         }
00449       }
00450     }
00451 
00452     for ( int i = 0; i < 3; ++i ) // loop on orthogonal faces of the unit tetra
00453     {
00454       if ( !sideFaces[i] ) continue;
00455       std::vector< double* >& sideFace = *sideFaces[i];
00456 
00457       std::size_t nbPoints = sideFace.size();
00458       if ( nbPoints == 0 )
00459         continue; // not intersected face at all - no cut off corners can be detected
00460 
00461       int ind1 = (i+1)%3, ind2 = (i+2)%3; // indices of coords on i-th tetra side
00462 
00463       int nbCutOnSide = 0;
00464       bool isSegmentOnEdge=false;
00465       for ( std::size_t ip = 0; ip < nbPoints; ++ip )
00466       {
00467         std::size_t isSegmentEnd = ( ip % 2 );
00468 
00469         double* p = sideFace[ ip ];
00470         double* p2 = isSegmentEnd ? 0 : sideFace[ip+1];
00471 
00472         if ( !isSegmentEnd )
00473           isSegmentOnEdge = false; // initialize
00474 
00475         int cutOffIndex = -1, passThIndex = -1;// no cut off neither pass through
00476         int pCut[] = { 0,0,0 }, pPass[] = { 0,0,0 };
00477 
00478         if ( epsilonEqual( p[ind1], 0.))
00479         {
00480           // point is on orthogonal edge
00481           if ( !isSegmentEnd && epsilonEqual( p2[ind1], 0. ))
00482             isSegmentOnEdge = true;
00483 
00484           if ( !isSegmentOnEdge )
00485           { // segment ends are on different edges
00486             pCut[ind2] = (int)isSegmentEnd; // believe that cutting triangles are well oriented
00487             cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2];
00488           }
00489           if ( epsilonEqual( p[ind2], 0.) || epsilonEqual( p[ind2], 1.))
00490           {
00491             pPass[ind2] = ( p[ind2] < 0.5 ) ? 0 : 1;
00492             passThIndex = pPass[0] + 2*pPass[1] + 3*pPass[2];
00493           }
00494         }
00495         else if ( epsilonEqual( p[ind2], 0.))
00496         {
00497           // point is on orthogonal edge
00498           if ( !isSegmentEnd && epsilonEqual( p2[ind2], 0. ))
00499             isSegmentOnEdge = true;
00500           if ( !isSegmentEnd )
00501           {// segment ends are on different edges
00502             pCut[ind1] = 1-(int)isSegmentEnd;
00503             cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2];
00504           }
00505           if ( epsilonEqual( p[ind1], 0.) || epsilonEqual( p[ind1], 1.))
00506           {
00507             pPass[ind1] = ( p[ind1] < 0.5 ) ? 0 : 1;
00508             passThIndex = pPass[0] + 2*pPass[1] + 3*pPass[2];
00509           }
00510         }
00511         else if ( epsilonEqual(p[ind1] + p[ind2], 1.0 ))
00512         {
00513           // point is on inclined edge
00514           if ( !isSegmentEnd && epsilonEqual(p2[ind1] + p2[ind2], 1.0 ))
00515             isSegmentOnEdge = true;
00516           if ( !isSegmentOnEdge )
00517           { //segment ends are on different edges
00518             pCut[ind1] = (int)isSegmentEnd;
00519             pCut[ind2] = 1-(int)isSegmentEnd;
00520             cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2];
00521           }
00522         }
00523         else
00524         {
00525           continue;
00526         }
00527         // remember cut off and passed through points
00528         if ( passThIndex >= 0 )
00529         {
00530           passedCorners[ passThIndex ] = true;
00531           if ( cutOffCorners[ passThIndex ] )
00532           {
00533             nbCutOffCorners--;
00534             cutOffCorners[ passThIndex ] = false;
00535 #ifdef DMP_UNITTETRAINTERSECTIONBARY
00536             std::cout << "PASS THROUGH " << passThIndex << std::endl;
00537 #endif
00538           }
00539         }
00540         if ( cutOffIndex >= 0 )
00541         {
00542           nbCutOnSide++;
00543           if ( !passedCorners[ cutOffIndex ] && !cutOffCorners[ cutOffIndex ] )
00544           {
00545             nbCutOffCorners++;
00546             cutOffCorners[ cutOffIndex ] = true;
00547           }
00548         }
00549       } // loop on points on a unit tetra side
00550 
00551       if ( nbCutOnSide == 0 && nbPoints <= 2 )
00552         continue; // one segment laying on edge at most
00553 
00554       if ( nbCutOffCorners == NB_TETRA_NODES )
00555         break; // all tetra corners are cut off
00556 
00557       if ( /*nbCutOnSide <= 2 &&*/ nbPoints >= 6 )
00558       {
00559         // at least 3 segments - all corners of a side are cut off
00560         for (int cutIndex = 0; cutIndex < NB_TETRA_NODES; ++cutIndex )
00561           if ( cutIndex != i+1 && !passedCorners[ cutIndex ] && !cutOffCorners[ cutIndex ])
00562             cutOffCorners[ cutIndex ] = ++nbCutOffCorners != 0 ;
00563       }
00564 
00565     }
00566     // loop on orthogonal faces of tetra
00567 
00568     // check if all corners are cut off on the inclined tetra side
00569     if ( sideFaces[ XYZ ] && sideFaces[ XYZ ]->size() >= 6 )
00570     {
00571       for (int cutIndex = 1; cutIndex < NB_TETRA_NODES; ++cutIndex )
00572         if ( !passedCorners[ cutIndex ] && !cutOffCorners[ cutIndex ])
00573           cutOffCorners[ cutIndex ] = ++nbCutOffCorners != 0 ;
00574     }
00575 
00576     // Add to faces on tetra sides the corners not cut off by segments of intersection polygons
00577     // ----------------------------------------------------------------------------------
00578     if ( nbCutOffCorners > 0 )
00579     {
00580       for ( int i = 0; i < NB_TETRA_SIDES; ++i )
00581       {
00582         if ( !sideFaces[ i ] ) continue;
00583         std::vector< double* >& sideFace = *sideFaces[i];
00584 
00585         int excludeCorner = (i + 1) % NB_TETRA_NODES;
00586         for ( int ic = 0; ic < NB_TETRA_NODES; ++ic )
00587         {
00588           if ( !cutOffCorners[ ic ] && ic != excludeCorner )
00589           {
00590             sideFace.push_back( new double[3] );
00591             copyVector3( origin, sideFace.back() );
00592             if ( ic )
00593               sideFace.back()[ ic-1 ] = 1.0;
00594           }
00595         }
00596       }
00597     }
00598 
00599 #ifdef DMP_UNITTETRAINTERSECTIONBARY
00600     std::cout << "**** after Add corners to sides " << std::endl;
00601     for ( int i = 0; i < NB_TETRA_SIDES; ++i )
00602     {
00603       std::cout << "\t Side " << i << std::endl;
00604       if ( !sideFaces[i] ) {
00605         std::cout << "\t cut by triagle" << std::endl;
00606       }
00607       else 
00608       {
00609         std::vector< double* >& sideFace = *sideFaces[i];
00610         for ( int i = 0; i < sideFace.size(); ++i )
00611         {
00612           double* p = sideFace[i];
00613           std::cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << std::endl;
00614         }
00615       }
00616     }
00617     std::cout << "Cut off corners: ";
00618     if ( nbCutOffCorners == 0 )
00619       std::cout << "NO";
00620     else 
00621       for ( int ic = 0; ic < NB_TETRA_NODES; ++ic )
00622         std::cout << cutOffCorners[ ic ];
00623     std::cout << std::endl;
00624 #endif
00625     // ------------------------------------------------------------------------
00626     // Sort corners of filled up faces on tetra sides and exclude equal points
00627     // ------------------------------------------------------------------------
00628 
00629     std::size_t iF = 0;
00630     for ( f = _faces.begin(); f != fEnd; ++f, ++iF )
00631     {
00632       std::vector< double* >&  face = *f;
00633       if ( face.size() >= 3 )
00634       {
00635         clearPolygons(); // free memory of _polygonA
00636         _polygonA = face;
00637         face.clear();
00638         face.reserve( _polygonA.size() );
00639         if ( iF >= nbIntersectPolygs )
00640         { // sort points of side faces
00641           calculatePolygonBarycenter( A, _barycenterA );
00642           setTriangleOnSide( (int)(iF-nbIntersectPolygs) );
00643           sortIntersectionPolygon( A, _barycenterA );
00644         }
00645         // exclude equal points
00646         std::vector< double* >::iterator v = _polygonA.begin(), vEnd = _polygonA.end();
00647         face.push_back( *v );
00648         *v = 0;
00649         for ( ++v; v != vEnd; ++v )
00650         {
00651           double* pPrev = face.back();
00652           double* p     = *v;
00653           if ( !samePoint( p, pPrev ))
00654           {
00655             face.push_back( p );
00656             *v = 0;
00657           }
00658         }
00659       }
00660       if ( face.size() < 3 )
00661       { // size could decrease
00662         clearPolygons(); // free memory of _polygonA
00663         _polygonA = face;
00664         face.clear();
00665       }
00666       else
00667       {
00668         nbPolyhedraFaces++;
00669       }
00670     }
00671 #ifdef DMP_UNITTETRAINTERSECTIONBARY
00672     std::cout << "**** after HEALING all faces " << std::endl;
00673     for (iF=0, f = _faces.begin(); f != fEnd; ++f, ++iF )
00674     {
00675       std::cout << "\t Side " << iF << std::endl;
00676       std::vector< double* >& sideFace = *f;
00677       for ( int i = 0; i < sideFace.size(); ++i )
00678       {
00679         double* p = sideFace[i];
00680         std::cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << std::endl;
00681       }
00682     }
00683 #endif
00684     return nbPolyhedraFaces;
00685   }
00686 
00687   //================================================================================
00693   //================================================================================
00694 
00695   void UnitTetraIntersectionBary::setTriangleOnSide(int iSide)
00696   {
00697     if ( iSide >= 3 )
00698       iSide = 0;
00699     for(int i = 0 ; i < 3 ; ++i) 
00700       {
00701         _coords[5*i] = _coords[5*i + 1] = _coords[5*i + 2] = 0.;
00702         if ( i != iSide )
00703           _coords[5*i + i] = 1.;
00704       }
00705   }
00706 
00707   //================================================================================
00711   //================================================================================
00712 
00713   void UnitTetraIntersectionBary::clearPolygons(bool andFaces)
00714   {
00715     for(std::vector<double*>::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it)
00716       {  delete[] *it;
00717         *it = 0; 
00718       }
00719     for(std::vector<double*>::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it)
00720       { 
00721         delete[] *it; 
00722         *it = 0; 
00723       }
00724 
00725     _polygonA.clear();
00726     _polygonB.clear();
00727 
00728     if ( andFaces )
00729       {
00730         std::list< std::vector< double* > >::iterator f = this->_faces.begin(), fEnd = this->_faces.end();
00731         for ( ; f != fEnd; ++f )
00732           {
00733             std::vector< double* >& polygon = *f;
00734             for(std::vector<double*>::iterator it = polygon.begin() ; it != polygon.end() ; ++it)
00735               { 
00736                 delete[] *it;
00737                 *it = 0;
00738               }
00739           }
00740         this->_faces.clear();
00741       }
00742   }
00743 
00744   //================================================================================
00748   //================================================================================
00749 
00750   UnitTetraIntersectionBary::~UnitTetraIntersectionBary()
00751   {
00752     clearPolygons(/*andFaces=*/true );
00753   }
00754 
00755 }