1 // Copyright (C) 2007-2008 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // File : UnitTetraIntersectionBary.cxx
20 // Created : Tue Dec 9 16:48:49 2008
21 // Author : Edward AGAPOV (eap)
23 #include "UnitTetraIntersectionBary.hxx"
25 #include "VectorUtils.hxx"
26 #include "InterpolationUtils.hxx"
27 #include "VolSurfFormulae.hxx"
29 #define NB_TETRA_SIDES 4
30 #define NB_TETRA_NODES 4
34 namespace INTERP_KERNEL
36 enum { _X=0, _Y, _Z };
38 inline bool samePoint( const double* p1, const double* p2 )
40 return ( p1[0]==p2[0] && p1[1]==p2[1] && p1[2]==p2[2]);
43 //================================================================================
45 * \brief Creates a ready-to-use tool
47 //================================================================================
49 UnitTetraIntersectionBary::UnitTetraIntersectionBary():TransformedTriangle()
54 //================================================================================
56 * \brief Initializes fields
58 //================================================================================
60 void UnitTetraIntersectionBary::init()
66 //================================================================================
68 * \brief Stores a part of triangle common with the unit tetrahedron
69 * \param triangle - triangle side of other cell
71 //================================================================================
73 void UnitTetraIntersectionBary::addSide(const TransformedTriangle& triangle)
75 _int_volume += triangle.getVolume();
77 double triNormal[3], polyNormal[3];
78 crossprod<3>( triangle.getCorner(P),triangle.getCorner(Q),triangle.getCorner(R), triNormal);
80 const std::vector<double*> * pPolygonA = &triangle.getPolygonA();
81 if ( pPolygonA->size() < 3 )
83 if ( !epsilonEqual( triNormal[_Z], 0 ))
84 return; // not vertical triangle does not intersect the unit tetra
86 // Vertical triangle. Use inherited methods of TransformedTriangle to
87 // calculate intesection polygon
88 *((TransformedTriangle*)this) = triangle; // copy triangle fields
91 calculateIntersectionPolygons();
92 if (this->_polygonA.size() < 3)
94 calculatePolygonBarycenter(A, _barycenterA);
95 sortIntersectionPolygon(A, _barycenterA);
96 pPolygonA = & _polygonA;
99 // check if polygon orientation is same as the one of triangle
100 vector<double*>::const_iterator p = pPolygonA->begin(), pEnd = pPolygonA->end();
103 while ( samePoint( p1, p2 ) && ++p != pEnd )
111 while ( samePoint( p2, p3 ) && ++p != pEnd )
118 crossprod<3>( p1, p2, p3, polyNormal );
119 bool reverse = ( dotprod<3>( triNormal, polyNormal ) < 0.0 );
122 _faces.push_back( vector< double* > () );
123 vector< double* >& faceCorner = _faces.back();
124 faceCorner.resize( pPolygonA->size()/* + 1*/ );
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] );
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] );
147 clearPolygons(); // free memory of _polygonA
148 _polygonA = faceCorner;
153 if ( i < (int)pPolygonA->size() )
154 faceCorner.resize( i );
158 cout << "**** addSide() " << _faces.size() << endl;
159 for ( int i = 0; i < faceCorner.size(); ++i )
161 double* p = faceCorner[i];
162 cout << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl;
168 //================================================================================
170 * \brief Computes and returns coordinates of barycentre
172 //================================================================================
174 bool UnitTetraIntersectionBary::getBary(double* baryCenter)
176 baryCenter[0] = baryCenter[1] = baryCenter[2] = -1.0;
177 if ( addSideFaces() < NB_TETRA_SIDES )
179 // tetra is not intersected
180 if ( _int_volume != 0.0 )
182 // tetra is fully inside the other cell
183 baryCenter[0] = baryCenter[1] = baryCenter[2] = 0.25;
184 _int_volume = 0.16666666666666666;
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.
197 baryCenter[0] = baryCenter[1] = baryCenter[2] = 0.;
199 list< vector< double* > >::iterator f = _faces.begin(), fEnd = _faces.end();
200 double * P = f->at(0);
202 for ( ++f; f != fEnd; ++f )
204 vector< double* >& polygon = *f;
205 if ( polygon.empty() )
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 )
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.
218 double bary[] = { 0, 0, 0 };
221 for ( v = polygon.begin(); v != vEnd ; ++v )
228 bary[0] /= polygon.size();
229 bary[1] /= polygon.size();
230 bary[2] /= polygon.size();
234 for ( int i = 0; i < (int)polygon.size(); ++i )
236 double* p1 = polygon[i];
237 double* p2 = polygon[(i+1)%polygon.size()];
238 vol += std::fabs( calculateVolumeForTetra( p1, p2, bary, P ));
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;
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;
255 //================================================================================
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
261 //================================================================================
263 int UnitTetraIntersectionBary::addSideFaces()
265 int nbPolyhedraFaces = 0;
267 if ( _faces.empty() )
268 return nbPolyhedraFaces;
270 // -------------------------------------------
271 // Detect polygons laying on sides of a tetra
272 // -------------------------------------------
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 )
279 vector< double* >& polygon = *f;
280 double coordSum[3] = {0,0,0};
281 for ( int i = 0; i < (int)polygon.size(); ++i )
283 double* p = polygon[i];
288 for ( int j = 0; j < 3 && !sideAdded[j]; ++j )
290 if ( coordSum[j] == 0 )
291 sideAdded[j] = bool( ++nbAddedSides );
293 if ( !sideAdded[3] &&
294 ( epsilonEqual( (coordSum[0]+coordSum[1]+coordSum[2]) / polygon.size(), 1. )))
295 sideAdded[3] = bool( ++nbAddedSides );
297 if ( nbAddedSides == NB_TETRA_SIDES )
300 // ---------------------------------------------------------------------------------
301 // Add segments of already added polygons to future polygonal faces on sides of tetra
302 // ---------------------------------------------------------------------------------
304 int nbIntersectPolygs = _faces.size();
306 vector< double* > * sideFaces[ 4 ]; // future polygons on sides of tetra
307 for ( int i = 0; i < NB_TETRA_SIDES; ++i )
310 if ( !sideAdded[ i ] )
312 _faces.push_back( vector< double* > () );
313 sideFaces[ i ] = &_faces.back();
316 f = _faces.begin(), fEnd = _faces.end();
317 for ( int iF = 0; iF < nbIntersectPolygs; ++f, ++iF ) // loop on added intersection polygons
319 vector< double* >& polygon = *f;
320 for ( int i = 0; i < (int)polygon.size(); ++i )
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 )
328 if ( !sideFaces[ j ] )
330 p1OnSide = ( p1[j] == 0 );
331 p2OnSide = ( p2[j] == 0 );
332 if ( p1OnSide && p2OnSide )
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() );
342 // check if the segment p1-p2 is on the inclined side
344 epsilonEqual( p1[_X] + p1[_Y] + p1[_Z], 1.0 ) &&
345 epsilonEqual( p2[_X] + p2[_Y] + p2[_Z], 1.0 ))
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() );
354 #ifdef DMP_ADDSIDEFACES
355 cout << "**** after Add segments to sides " << endl;
356 for ( int i = 0; i < NB_TETRA_SIDES; ++i )
358 cout << "\t Side " << i << endl;
361 cout << "\t cut by triagle" << endl;
365 vector< double* >& sideFace = *sideFaces[i];
366 for ( int i = 0; i < sideFace.size(); ++i )
368 double* p = sideFace[i];
369 cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl;
375 // ---------------------------------------------------------------------------
376 // Make closed polygons on tetra sides by adding not cut off corners of tetra
377 // ---------------------------------------------------------------------------
379 double origin[3] = { 0,0,0 };
381 // find corners of tetra cut off by triangles of other tetra
382 // ---------------------------------------------------------
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);
388 int nbCutOffCorners = 0;
389 for ( int i = 0; i < 3; ++i ) // loop on orthogonal faces of the unit tetra
391 if ( !sideFaces[i] ) continue;
392 vector< double* >& sideFace = *sideFaces[i];
394 int nbPoints = sideFace.size();
396 continue; // not intersected face at all - no cut off corners can be detected
398 int ind1 = (i+1)%3, ind2 = (i+2)%3; // indices of coords on i-th tetra side
401 bool isSegmentOnEdge;
402 for ( int ip = 0; ip < nbPoints; ++ip )
404 int isSegmentEnd = ( ip % 2 );
406 double* p = sideFace[ ip ];
407 double* p2 = isSegmentEnd ? 0 : sideFace[ip+1];
410 isSegmentOnEdge = false; // initialize
412 int cutOffIndex = -1, passThIndex = -1;// no cut off neither pass through
413 int pCut[] = { 0,0,0 }, pPass[] = { 0,0,0 };
417 // point is in on orthogonal edge
418 if ( !isSegmentEnd && p2[ind1]==0 )
419 isSegmentOnEdge = true;
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];
426 if ( p[ind2]==0. || p[ind2]==1.)
428 pPass[ind2] = int(p[ind2]);
429 passThIndex = pPass[0] + 2*pPass[1] + 3*pPass[2];
432 else if ( p[ind2]==0.)
434 // point is on orthogonal edge
435 if ( !isSegmentEnd && p2[ind2]==0 )
436 isSegmentOnEdge = true;
438 {// segment ends are on different edges
439 pCut[ind1] = 1-isSegmentEnd;
440 cutOffIndex = pCut[0] + 2*pCut[1] + 3*pCut[2];
442 if ( p[ind1]==0. || p[ind1]==1.)
444 pPass[ind1] = int(p[ind1]);
445 passThIndex = pPass[0] + 2*pPass[1] + 3*pPass[2];
448 else if ( epsilonEqual(p[ind1] + p[ind2], 1.0 ))
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];
464 // remember cut off and passed through points
465 if ( passThIndex >= 0. )
467 passedCorners[ passThIndex ] = true;
468 if ( cutOffCorners[ passThIndex ] )
471 cutOffCorners[ passThIndex ] = false;
474 if ( cutOffIndex >= 0. )
477 if ( !passedCorners[ cutOffIndex ] && !cutOffCorners[ cutOffIndex ] )
480 cutOffCorners[ cutOffIndex ] = true;
483 } // loop on points on a unit tetra side
485 if ( nbCutOnSide == 0 && nbPoints <= 2 )
486 continue; // one segment laying on edge at most
488 if ( nbCutOffCorners == NB_TETRA_NODES )
489 break; // all tetra corners are cut off
491 if ( /*nbCutOnSide <= 2 &&*/ nbPoints >= 6 )
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 );
499 } // loop on orthogonal faces of tetra
501 // check if all corners are cut off on the inclined tetra side
502 if ( sideFaces[ XYZ ] && sideFaces[ XYZ ]->size() >= 6 )
504 for (int cutIndex = 1; cutIndex < NB_TETRA_NODES; ++cutIndex )
505 if ( !passedCorners[ cutIndex ] && !cutOffCorners[ cutIndex ])
506 cutOffCorners[ cutIndex ] = bool ( ++nbCutOffCorners );
509 // Add to faces on tetra sides the corners not cut off by segments of intersection polygons
510 // ----------------------------------------------------------------------------------
511 if ( nbCutOffCorners > 0 )
513 for ( int i = 0; i < NB_TETRA_SIDES; ++i )
515 if ( !sideFaces[ i ] ) continue;
516 vector< double* >& sideFace = *sideFaces[i];
518 int excludeCorner = (i + 1) % NB_TETRA_NODES;
519 for ( int ic = 0; ic < NB_TETRA_NODES; ++ic )
521 if ( !cutOffCorners[ ic ] && ic != excludeCorner )
523 sideFace.push_back( new double[3] );
524 copyVector3( origin, sideFace.back() );
526 sideFace.back()[ ic-1 ] = 1.0;
532 #ifdef DMP_ADDSIDEFACES
533 cout << "**** after Add corners to sides " << endl;
534 for ( int i = 0; i < NB_TETRA_SIDES; ++i )
536 cout << "\t Side " << i << endl;
537 if ( !sideFaces[i] ) {
538 cout << "\t cut by triagle" << endl;
542 vector< double* >& sideFace = *sideFaces[i];
543 for ( int i = 0; i < sideFace.size(); ++i )
545 double* p = sideFace[i];
546 cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl;
550 cout << "Cut off corners: ";
551 if ( nbCutOffCorners == 0 )
554 for ( int ic = 0; ic < NB_TETRA_NODES; ++ic )
555 cout << cutOffCorners[ ic ];
558 // ------------------------------------------------------------------------
559 // Sort corners of filled up faces on tetra sides and exclude equal points
560 // ------------------------------------------------------------------------
563 for ( f = _faces.begin(); f != fEnd; ++f, ++iF )
565 vector< double* >& face = *f;
566 if ( face.size() >= 3 )
568 clearPolygons(); // free memory of _polygonA
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 );
578 // exclude equal points
579 vector< double* >::iterator v = _polygonA.begin(), vEnd = _polygonA.end();
580 face.push_back( *v );
582 for ( ++v; v != vEnd; ++v )
584 double* pPrev = face.back();
586 if ( !samePoint( p, pPrev ))
593 if ( face.size() < 3 )
594 { // size could decrease
595 clearPolygons(); // free memory of _polygonA
604 #ifdef DMP_ADDSIDEFACES
605 cout << "**** after HEALING all faces " << endl;
606 for (iF=0, f = _faces.begin(); f != fEnd; ++f, ++iF )
608 cout << "\t Side " << iF << endl;
609 vector< double* >& sideFace = *f;
610 for ( int i = 0; i < sideFace.size(); ++i )
612 double* p = sideFace[i];
613 cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << endl;
617 return nbPolyhedraFaces;
620 //================================================================================
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)
626 //================================================================================
628 void UnitTetraIntersectionBary::setTriangleOnSide(int iSide)
632 for(int i = 0 ; i < 3 ; ++i)
634 _coords[5*i] = _coords[5*i + 1] = _coords[5*i + 2] = 0.;
636 _coords[5*i + i] = 1.;
640 //================================================================================
642 * \brief Free memory of polygons
644 //================================================================================
646 void UnitTetraIntersectionBary::clearPolygons(bool andFaces)
648 for(vector<double*>::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it)
652 for(vector<double*>::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it)
663 list< vector< double* > >::iterator f = this->_faces.begin(), fEnd = this->_faces.end();
664 for ( ; f != fEnd; ++f )
666 vector< double* >& polygon = *f;
667 for(vector<double*>::iterator it = polygon.begin() ; it != polygon.end() ; ++it)
673 this->_faces.clear();
677 //================================================================================
679 * \brief Destructor clears coordinates of faces
681 //================================================================================
683 UnitTetraIntersectionBary::~UnitTetraIntersectionBary()
685 clearPolygons(/*andFaces=*/true );