From: eap Date: Tue, 29 Sep 2009 06:55:06 +0000 (+0000) Subject: 0020440: [CEA 349] P1P0 barycentric interpolators X-Git-Tag: V5_1_main_FINAL~339 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=0b34b97e48ef51afe0382d7b86648e6bdb33f417;p=tools%2Fmedcoupling.git 0020440: [CEA 349] P1P0 barycentric interpolators Fixing bugs - UnitTetraIntersectionBary(); + UnitTetraIntersectionBary(bool isTetraInversed=false); - void init(); + void init(bool isTetraInversed=false); + std::vector< std::vector< double > > _polyNormals; --- diff --git a/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx b/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx index 73efea7df..7c5d5735b 100644 --- a/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx +++ b/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx @@ -29,6 +29,8 @@ #define NB_TETRA_SIDES 4 #define NB_TETRA_NODES 4 +//#define DMP_UNITTETRAINTERSECTIONBARY + using namespace std; namespace INTERP_KERNEL @@ -37,7 +39,8 @@ namespace INTERP_KERNEL inline bool samePoint( const double* p1, const double* p2 ) { - return ( p1[0]==p2[0] && p1[1]==p2[1] && p1[2]==p2[2]); + //return ( p1[0]==p2[0] && p1[1]==p2[1] && p1[2]==p2[2]); + return ( epsilonEqual( p1[0], p2[0]) && epsilonEqual( p1[1], p2[1]) && epsilonEqual( p1[2], p2[2])); } //================================================================================ @@ -46,9 +49,9 @@ namespace INTERP_KERNEL */ //================================================================================ - UnitTetraIntersectionBary::UnitTetraIntersectionBary():TransformedTriangle() + UnitTetraIntersectionBary::UnitTetraIntersectionBary(bool isTetraInversed) + :TransformedTriangle(),_int_volume(0),_isTetraInversed( isTetraInversed ) { - _int_volume = 0; //init(); } //================================================================================ @@ -57,10 +60,12 @@ namespace INTERP_KERNEL */ //================================================================================ - void UnitTetraIntersectionBary::init() + void UnitTetraIntersectionBary::init(bool isTetraInversed) { _int_volume = 0; + _isTetraInversed = isTetraInversed; _faces.clear(); + _polyNormals.clear(); } //================================================================================ @@ -117,6 +122,7 @@ namespace INTERP_KERNEL } crossprod<3>( p1, p2, p3, polyNormal ); bool reverse = ( dotprod<3>( triNormal, polyNormal ) < 0.0 ); + if (_isTetraInversed) reverse = !reverse; // store polygon _faces.push_back( vector< double* > () ); @@ -132,6 +138,9 @@ namespace INTERP_KERNEL copyVector3( *polyF, faceCorner[i] = new double[3] ); else --i; + polyNormal[0] *= -1.; + polyNormal[1] *= -1.; + polyNormal[2] *= -1.; } else { @@ -152,9 +161,13 @@ namespace INTERP_KERNEL { if ( i < (int)pPolygonA->size() ) faceCorner.resize( i ); + + if ( _polyNormals.empty() ) + _polyNormals.reserve(4); + _polyNormals.push_back( vector< double >( polyNormal, polyNormal+3 )); } -#ifdef DMP_ADDSIDE +#ifdef DMP_UNITTETRAINTERSECTIONBARY cout << "**** addSide() " << _faces.size() << endl; for ( int i = 0; i < faceCorner.size(); ++i ) { @@ -249,6 +262,10 @@ namespace INTERP_KERNEL baryCenter[1] /= _int_volume; baryCenter[2] /= _int_volume; +#ifdef DMP_UNITTETRAINTERSECTIONBARY + cout << "**** Barycenter " << baryCenter[0] <<", "<< baryCenter[1] <<", "<< baryCenter[2] + << "\t **** Volume " << _int_volume << endl; +#endif return true; } @@ -275,25 +292,25 @@ namespace INTERP_KERNEL int nbAddedSides = 0; list< vector< double* > >::iterator f = _faces.begin(), fEnd = _faces.end(); for ( ; f != fEnd; ++f ) + { + vector< double* >& polygon = *f; + double coordSum[3] = {0,0,0}; + for ( int i = 0; i < (int)polygon.size(); ++i ) { - vector< double* >& polygon = *f; - double coordSum[3] = {0,0,0}; - for ( int i = 0; i < (int)polygon.size(); ++i ) - { - double* p = polygon[i]; - coordSum[0] += p[0]; - coordSum[1] += p[1]; - coordSum[2] += p[2]; - } - for ( int j = 0; j < 3 && !sideAdded[j]; ++j ) - { - if ( coordSum[j] == 0 ) - sideAdded[j] = bool( ++nbAddedSides ); - } - if ( !sideAdded[3] && - ( epsilonEqual( (coordSum[0]+coordSum[1]+coordSum[2]) / polygon.size(), 1. ))) - sideAdded[3] = bool( ++nbAddedSides ); + double* p = polygon[i]; + coordSum[0] += p[0]; + coordSum[1] += p[1]; + coordSum[2] += p[2]; } + for ( int j = 0; j < 3 && !sideAdded[j]; ++j ) + { + if ( coordSum[j] == 0 ) + sideAdded[j] = bool( ++nbAddedSides ); + } + if ( !sideAdded[3] && + ( epsilonEqual( (coordSum[0]+coordSum[1]+coordSum[2]) / polygon.size(), 1. ))) + sideAdded[3] = bool( ++nbAddedSides ); + } if ( nbAddedSides == NB_TETRA_SIDES ) return nbAddedSides; @@ -305,71 +322,71 @@ namespace INTERP_KERNEL vector< double* > * sideFaces[ 4 ]; // future polygons on sides of tetra for ( int i = 0; i < NB_TETRA_SIDES; ++i ) + { + sideFaces[ i ]=0; + if ( !sideAdded[ i ] ) { - sideFaces[ i ]=0; - if ( !sideAdded[ i ] ) - { - _faces.push_back( vector< double* > () ); - sideFaces[ i ] = &_faces.back(); - } + _faces.push_back( vector< double* > () ); + sideFaces[ i ] = &_faces.back(); } + } f = _faces.begin(), fEnd = _faces.end(); for ( int iF = 0; iF < nbIntersectPolygs; ++f, ++iF ) // loop on added intersection polygons + { + vector< double* >& polygon = *f; + for ( int i = 0; i < (int)polygon.size(); ++i ) { - vector< double* >& polygon = *f; - for ( int i = 0; i < (int)polygon.size(); ++i ) + // segment ends + double* p1 = polygon[i]; + double* p2 = polygon[(i+1)%polygon.size()]; + bool p1OnSide, p2OnSide;//, onZeroSide = false; + for ( int j = 0; j < 3; ++j ) + { + if ( !sideFaces[ j ] ) + continue; + p1OnSide = ( p1[j] == 0. ); + p2OnSide = ( p2[j] == 0. ); + if ( p1OnSide && p2OnSide ) { - // segment ends - double* p1 = polygon[i]; - double* p2 = polygon[(i+1)%polygon.size()]; - bool p1OnSide, p2OnSide;//, onZeroSide = false; - for ( int j = 0; j < 3; ++j ) - { - if ( !sideFaces[ j ] ) - continue; - p1OnSide = ( p1[j] == 0 ); - p2OnSide = ( p2[j] == 0 ); - if ( p1OnSide && p2OnSide ) - { - // segment p1-p2 is on j-th orthogonal side of tetra - sideFaces[j]->push_back( new double[3] ); - copyVector3( p1, sideFaces[j]->back() ); - sideFaces[j]->push_back( new double[3] ); - copyVector3( p2, sideFaces[j]->back() ); - //break; - } - } - // check if the segment p1-p2 is on the inclined side - if ( sideFaces[3] && - epsilonEqual( p1[_X] + p1[_Y] + p1[_Z], 1.0 ) && - epsilonEqual( p2[_X] + p2[_Y] + p2[_Z], 1.0 )) - { - sideFaces[3]->push_back( new double[3] ); - copyVector3( p1, sideFaces[3]->back() ); - sideFaces[3]->push_back( new double[3] ); - copyVector3( p2, sideFaces[3]->back() ); - } + // segment p1-p2 is on j-th orthogonal side of tetra + sideFaces[j]->push_back( new double[3] ); + copyVector3( p1, sideFaces[j]->back() ); + sideFaces[j]->push_back( new double[3] ); + copyVector3( p2, sideFaces[j]->back() ); + //break; } + } + // check if the segment p1-p2 is on the inclined side + if ( sideFaces[3] && + epsilonEqual( p1[_X] + p1[_Y] + p1[_Z], 1.0 ) && + epsilonEqual( p2[_X] + p2[_Y] + p2[_Z], 1.0 )) + { + sideFaces[3]->push_back( new double[3] ); + copyVector3( p1, sideFaces[3]->back() ); + sideFaces[3]->push_back( new double[3] ); + copyVector3( p2, sideFaces[3]->back() ); + } } -#ifdef DMP_ADDSIDEFACES + } +#ifdef DMP_UNITTETRAINTERSECTIONBARY cout << "**** after Add segments to sides " << endl; for ( int i = 0; i < NB_TETRA_SIDES; ++i ) + { + cout << "\t Side " << i << endl; + if ( !sideFaces[i] ) { - cout << "\t Side " << i << endl; - if ( !sideFaces[i] ) - { - cout << "\t cut by triagle" << endl; - } - else - { - vector< double* >& sideFace = *sideFaces[i]; - for ( int i = 0; i < sideFace.size(); ++i ) - { - double* p = sideFace[i]; - cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl; - } - } + cout << "\t cut by triagle" << endl; + } + else + { + vector< double* >& sideFace = *sideFaces[i]; + for ( int i = 0; i < sideFace.size(); ++i ) + { + double* p = sideFace[i]; + cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl; + } } + } #endif // --------------------------------------------------------------------------- @@ -383,170 +400,200 @@ namespace INTERP_KERNEL // corners are coded like this: index = 1*X + 2*Y + 3*Z // (0,0,0) -> index == 0; (0,0,1) -> index == 3 - vector< int > cutOffCorners(NB_TETRA_NODES, false), passedCorners(NB_TETRA_NODES, false); + int cutOffCorners[NB_TETRA_NODES] = { false, false, false, false }; + int passedCorners[NB_TETRA_NODES] = { false, false, false, false }; + // find cutOffCorners by normals of intersection polygons int nbCutOffCorners = 0; - for ( int i = 0; i < 3; ++i ) // loop on orthogonal faces of the unit tetra + for ( int ic = 0; ic < NB_TETRA_NODES; ++ic ) + { + //if ( !passedCorners[ ic ] && !cutOffCorners[ ic ]) { - if ( !sideFaces[i] ) continue; - vector< double* >& sideFace = *sideFaces[i]; +// double corner[3] = { 0., 0., 0. }; +// if ( ic ) corner[ ic-1 ] = 1.0; - int nbPoints = sideFace.size(); - if ( nbPoints == 0 ) - continue; // not intersected face at all - no cut off corners can be detected + f = _faces.begin(), fEnd = _faces.end(); + for ( int iF = 0; iF < nbIntersectPolygs; ++f, ++iF ) // loop on added intersection polygons + { + vector< double* >& polygon = *f; - int ind1 = (i+1)%3, ind2 = (i+2)%3; // indices of coords on i-th tetra side + double corner2Poly[3] = { polygon[0][0], polygon[0][1], polygon[0][2] }; + if ( ic ) corner2Poly[ ic-1 ] -= 1.0; - int nbCutOnSide = 0; - bool isSegmentOnEdge; - for ( int ip = 0; ip < nbPoints; ++ip ) + // _polyNormals are outside of a tetrahedron + double dot = dotprod<3>( corner2Poly, &_polyNormals[iF][0] ); + if ( dot < -DEFAULT_ABS_TOL ) { - int isSegmentEnd = ( ip % 2 ); - - double* p = sideFace[ ip ]; - double* p2 = isSegmentEnd ? 0 : sideFace[ip+1]; - - if ( !isSegmentEnd ) - isSegmentOnEdge = false; // initialize - - int cutOffIndex = -1, passThIndex = -1;// no cut off neither pass through - int pCut[] = { 0,0,0 }, pPass[] = { 0,0,0 }; - - if ( p[ind1]==0.) - { - // point is in on orthogonal edge - if ( !isSegmentEnd && p2[ind1]==0 ) - isSegmentOnEdge = true; - - if ( !isSegmentOnEdge ) - { // segment ends are on different edges - pCut[ind2] = isSegmentEnd; // believe that cutting triangles are well oriented - cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2]; - } - if ( p[ind2]==0. || p[ind2]==1.) - { - pPass[ind2] = int(p[ind2]); - passThIndex = pPass[0] + 2*pPass[1] + 3*pPass[2]; - } - } - else if ( p[ind2]==0.) - { - // point is on orthogonal edge - if ( !isSegmentEnd && p2[ind2]==0 ) - isSegmentOnEdge = true; - if ( !isSegmentEnd ) - {// segment ends are on different edges - pCut[ind1] = 1-isSegmentEnd; - cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2]; - } - if ( p[ind1]==0. || p[ind1]==1.) - { - pPass[ind1] = int(p[ind1]); - passThIndex = pPass[0] + 2*pPass[1] + 3*pPass[2]; - } - } - else if ( epsilonEqual(p[ind1] + p[ind2], 1.0 )) - { - // point is on inclined edge - if ( !isSegmentEnd && epsilonEqual(p2[ind1] + p2[ind2], 1.0 )) - isSegmentOnEdge = true; - if ( !isSegmentOnEdge ) - { //segment ends are on different edges - pCut[ind1] = isSegmentEnd; - pCut[ind2] = 1-isSegmentEnd; - cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2]; - } - } - else - { - continue; - } - // remember cut off and passed through points - if ( passThIndex >= 0. ) - { - passedCorners[ passThIndex ] = true; - if ( cutOffCorners[ passThIndex ] ) - { - nbCutOffCorners--; - cutOffCorners[ passThIndex ] = false; - } - } - if ( cutOffIndex >= 0. ) - { - nbCutOnSide++; - if ( !passedCorners[ cutOffIndex ] && !cutOffCorners[ cutOffIndex ] ) - { - nbCutOffCorners++; - cutOffCorners[ cutOffIndex ] = true; - } - } - } // loop on points on a unit tetra side + cutOffCorners[ ic ] = true; + nbCutOffCorners++; + break; + } + } + } + } + + for ( int i = 0; i < 3; ++i ) // loop on orthogonal faces of the unit tetra + { + if ( !sideFaces[i] ) continue; + vector< double* >& sideFace = *sideFaces[i]; + + int nbPoints = sideFace.size(); + if ( nbPoints == 0 ) + continue; // not intersected face at all - no cut off corners can be detected + + int ind1 = (i+1)%3, ind2 = (i+2)%3; // indices of coords on i-th tetra side + + int nbCutOnSide = 0; + bool isSegmentOnEdge; + for ( int ip = 0; ip < nbPoints; ++ip ) + { + int isSegmentEnd = ( ip % 2 ); + + double* p = sideFace[ ip ]; + double* p2 = isSegmentEnd ? 0 : sideFace[ip+1]; - if ( nbCutOnSide == 0 && nbPoints <= 2 ) - continue; // one segment laying on edge at most + if ( !isSegmentEnd ) + isSegmentOnEdge = false; // initialize - if ( nbCutOffCorners == NB_TETRA_NODES ) - break; // all tetra corners are cut off + int cutOffIndex = -1, passThIndex = -1;// no cut off neither pass through + int pCut[] = { 0,0,0 }, pPass[] = { 0,0,0 }; - if ( /*nbCutOnSide <= 2 &&*/ nbPoints >= 6 ) + if ( p[ind1]==0.) + { + // point is on orthogonal edge + if ( !isSegmentEnd && p2[ind1]==0. ) + isSegmentOnEdge = true; + + if ( !isSegmentOnEdge ) + { // segment ends are on different edges + pCut[ind2] = isSegmentEnd; // believe that cutting triangles are well oriented + cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2]; + } + if ( p[ind2]==0. || p[ind2]==1.) + { + pPass[ind2] = int(p[ind2]); + passThIndex = pPass[0] + 2*pPass[1] + 3*pPass[2]; + } + } + else if ( p[ind2]==0.) + { + // point is on orthogonal edge + if ( !isSegmentEnd && p2[ind2]==0. ) + isSegmentOnEdge = true; + if ( !isSegmentEnd ) + {// segment ends are on different edges + pCut[ind1] = 1-isSegmentEnd; + cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2]; + } + if ( p[ind1]==0. || p[ind1]==1.) + { + pPass[ind1] = int(p[ind1]); + passThIndex = pPass[0] + 2*pPass[1] + 3*pPass[2]; + } + } + else if ( epsilonEqual(p[ind1] + p[ind2], 1.0 )) + { + // point is on inclined edge + if ( !isSegmentEnd && epsilonEqual(p2[ind1] + p2[ind2], 1.0 )) + isSegmentOnEdge = true; + if ( !isSegmentOnEdge ) + { //segment ends are on different edges + pCut[ind1] = isSegmentEnd; + pCut[ind2] = 1-isSegmentEnd; + cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2]; + } + } + else + { + continue; + } + // remember cut off and passed through points + if ( passThIndex >= 0. ) + { + passedCorners[ passThIndex ] = true; + if ( cutOffCorners[ passThIndex ] ) + { + nbCutOffCorners--; + cutOffCorners[ passThIndex ] = false; + } + } + if ( cutOffIndex >= 0. ) + { + nbCutOnSide++; + if ( !passedCorners[ cutOffIndex ] && !cutOffCorners[ cutOffIndex ] ) { - // at least 3 segments - all corners of a side are cut off - for (int cutIndex = 0; cutIndex < NB_TETRA_NODES; ++cutIndex ) - if ( cutIndex != i+1 && !passedCorners[ cutIndex ] && !cutOffCorners[ cutIndex ]) - cutOffCorners[ cutIndex ] = bool ( ++nbCutOffCorners ); + nbCutOffCorners++; + cutOffCorners[ cutOffIndex ] = true; } + } + } // loop on points on a unit tetra side - } // loop on orthogonal faces of tetra + if ( nbCutOnSide == 0 && nbPoints <= 2 ) + continue; // one segment laying on edge at most - // check if all corners are cut off on the inclined tetra side - if ( sideFaces[ XYZ ] && sideFaces[ XYZ ]->size() >= 6 ) + if ( nbCutOffCorners == NB_TETRA_NODES ) + break; // all tetra corners are cut off + + if ( /*nbCutOnSide <= 2 &&*/ nbPoints >= 6 ) { - for (int cutIndex = 1; cutIndex < NB_TETRA_NODES; ++cutIndex ) - if ( !passedCorners[ cutIndex ] && !cutOffCorners[ cutIndex ]) + // at least 3 segments - all corners of a side are cut off + for (int cutIndex = 0; cutIndex < NB_TETRA_NODES; ++cutIndex ) + if ( cutIndex != i+1 && !passedCorners[ cutIndex ] && !cutOffCorners[ cutIndex ]) cutOffCorners[ cutIndex ] = bool ( ++nbCutOffCorners ); } + } + // loop on orthogonal faces of tetra + + // check if all corners are cut off on the inclined tetra side + if ( sideFaces[ XYZ ] && sideFaces[ XYZ ]->size() >= 6 ) + { + for (int cutIndex = 1; cutIndex < NB_TETRA_NODES; ++cutIndex ) + if ( !passedCorners[ cutIndex ] && !cutOffCorners[ cutIndex ]) + cutOffCorners[ cutIndex ] = bool ( ++nbCutOffCorners ); + } + // Add to faces on tetra sides the corners not cut off by segments of intersection polygons // ---------------------------------------------------------------------------------- if ( nbCutOffCorners > 0 ) + { + for ( int i = 0; i < NB_TETRA_SIDES; ++i ) { - for ( int i = 0; i < NB_TETRA_SIDES; ++i ) + if ( !sideFaces[ i ] ) continue; + vector< double* >& sideFace = *sideFaces[i]; + + int excludeCorner = (i + 1) % NB_TETRA_NODES; + for ( int ic = 0; ic < NB_TETRA_NODES; ++ic ) + { + if ( !cutOffCorners[ ic ] && ic != excludeCorner ) { - if ( !sideFaces[ i ] ) continue; - vector< double* >& sideFace = *sideFaces[i]; - - int excludeCorner = (i + 1) % NB_TETRA_NODES; - for ( int ic = 0; ic < NB_TETRA_NODES; ++ic ) - { - if ( !cutOffCorners[ ic ] && ic != excludeCorner ) - { - sideFace.push_back( new double[3] ); - copyVector3( origin, sideFace.back() ); - if ( ic ) - sideFace.back()[ ic-1 ] = 1.0; - } - } + sideFace.push_back( new double[3] ); + copyVector3( origin, sideFace.back() ); + if ( ic ) + sideFace.back()[ ic-1 ] = 1.0; } + } } + } -#ifdef DMP_ADDSIDEFACES +#ifdef DMP_UNITTETRAINTERSECTIONBARY cout << "**** after Add corners to sides " << endl; for ( int i = 0; i < NB_TETRA_SIDES; ++i ) + { + cout << "\t Side " << i << endl; + if ( !sideFaces[i] ) { + cout << "\t cut by triagle" << endl; + } + else { - cout << "\t Side " << i << endl; - if ( !sideFaces[i] ) { - cout << "\t cut by triagle" << endl; + vector< double* >& sideFace = *sideFaces[i]; + for ( int i = 0; i < sideFace.size(); ++i ) + { + double* p = sideFace[i]; + cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl; } - else - { - vector< double* >& sideFace = *sideFaces[i]; - for ( int i = 0; i < sideFace.size(); ++i ) - { - double* p = sideFace[i]; - cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl; - } - } } + } cout << "Cut off corners: "; if ( nbCutOffCorners == 0 ) cout << "NO"; @@ -561,58 +608,58 @@ namespace INTERP_KERNEL int iF = 0; for ( f = _faces.begin(); f != fEnd; ++f, ++iF ) + { + vector< double* >& face = *f; + if ( face.size() >= 3 ) { - vector< double* >& face = *f; - if ( face.size() >= 3 ) + clearPolygons(); // free memory of _polygonA + _polygonA = face; + face.clear(); + face.reserve( _polygonA.size() ); + if ( iF >= nbIntersectPolygs ) + { // sort points of side faces + calculatePolygonBarycenter( A, _barycenterA ); + setTriangleOnSide( iF - nbIntersectPolygs ); + sortIntersectionPolygon( A, _barycenterA ); + } + // exclude equal points + vector< double* >::iterator v = _polygonA.begin(), vEnd = _polygonA.end(); + face.push_back( *v ); + *v = 0; + for ( ++v; v != vEnd; ++v ) + { + double* pPrev = face.back(); + double* p = *v; + if ( !samePoint( p, pPrev )) { - clearPolygons(); // free memory of _polygonA - _polygonA = face; - face.clear(); - face.reserve( _polygonA.size() ); - if ( iF >= nbIntersectPolygs ) - { // sort points of side faces - calculatePolygonBarycenter( A, _barycenterA ); - setTriangleOnSide( iF - nbIntersectPolygs ); - sortIntersectionPolygon( A, _barycenterA ); - } - // exclude equal points - vector< double* >::iterator v = _polygonA.begin(), vEnd = _polygonA.end(); - face.push_back( *v ); + face.push_back( p ); *v = 0; - for ( ++v; v != vEnd; ++v ) - { - double* pPrev = face.back(); - double* p = *v; - if ( !samePoint( p, pPrev )) - { - face.push_back( p ); - *v = 0; - } - } - } - if ( face.size() < 3 ) - { // size could decrease - clearPolygons(); // free memory of _polygonA - _polygonA = face; - face.clear(); - } - else - { - nbPolyhedraFaces++; } + } + } + if ( face.size() < 3 ) + { // size could decrease + clearPolygons(); // free memory of _polygonA + _polygonA = face; + face.clear(); + } + else + { + nbPolyhedraFaces++; } -#ifdef DMP_ADDSIDEFACES + } +#ifdef DMP_UNITTETRAINTERSECTIONBARY cout << "**** after HEALING all faces " << endl; for (iF=0, f = _faces.begin(); f != fEnd; ++f, ++iF ) + { + cout << "\t Side " << iF << endl; + vector< double* >& sideFace = *f; + for ( int i = 0; i < sideFace.size(); ++i ) { - cout << "\t Side " << iF << endl; - vector< double* >& sideFace = *f; - for ( int i = 0; i < sideFace.size(); ++i ) - { - double* p = sideFace[i]; - cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl; - } + double* p = sideFace[i]; + cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl; } + } #endif return nbPolyhedraFaces; } diff --git a/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx b/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx index 41aa72148..f03cab4a5 100644 --- a/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx +++ b/src/INTERP_KERNEL/UnitTetraIntersectionBary.hxx @@ -33,9 +33,9 @@ namespace INTERP_KERNEL class INTERPKERNEL_EXPORT UnitTetraIntersectionBary : protected TransformedTriangle { public: - UnitTetraIntersectionBary(); + UnitTetraIntersectionBary(bool isTetraInversed=false); - void init(); + void init(bool isTetraInversed=false); /*! * \brief Stores a part of triangle common with the unit tetrahedron * \param triangle - triangle side of other cell, whose calculateIntersectionVolume() @@ -68,7 +68,10 @@ namespace INTERP_KERNEL double _int_volume; /// faces of intersection polyhedron - std::list< std::vector< double* > > _faces; + std::list< std::vector< double* > > _faces; + std::vector< std::vector< double > > _polyNormals; + + bool _isTetraInversed; }; }