1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 // File : SMDS_VolumeTool.cxx
24 // Created : Tue Jul 13 12:22:13 2004
25 // Author : Edward AGAPOV (eap)
28 #pragma warning(disable:4786)
31 #include "SMDS_VolumeTool.hxx"
33 #include "SMDS_MeshElement.hxx"
34 #include "SMDS_MeshNode.hxx"
35 #include "SMDS_VtkVolume.hxx"
36 #include "SMDS_Mesh.hxx"
38 #include "utilities.h"
50 // ======================================================
51 // Node indices in faces depending on volume orientation
52 // making most faces normals external
53 // ======================================================
54 // For all elements, 0-th face is bottom based on the first nodes.
55 // For prismatic elements (tetra,hexa,prisms), 1-th face is a top one.
56 // For all elements, side faces follow order of bottom nodes
57 // ======================================================
65 // N0 +---|---+ N1 TETRAHEDRON
73 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
74 { 0, 1, 2, 0 }, // All faces have external normals
78 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
79 { 0, 2, 1, 0 }, // All faces have external normals
83 static int Tetra_nbN [] = { 3, 3, 3, 3 };
88 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
89 { 0, 1, 2, 3, 0 }, // All faces have external normals
95 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
96 { 0, 3, 2, 1, 0 }, // All faces but a bottom have external normals
101 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
112 // | / \ | PENTAHEDRON
118 static int Penta_F [5][5] = { // FORWARD
119 { 0, 1, 2, 0, 0 }, // All faces have external normals
120 { 3, 5, 4, 3, 3 }, // 0 is bottom, 1 is top face
124 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
130 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
137 // N4+----------+N7 |
138 // | | | | HEXAHEDRON
139 // | N1+------|---+N2
145 static int Hexa_F [6][5] = { // FORWARD
147 { 4, 7, 6, 5, 4 }, // all face normals are external
152 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
154 { 4, 5, 6, 7, 4 }, // all face normals are external
159 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
160 static int Hexa_oppF[] = { 1, 0, 4, 5, 2, 3 }; // oppopsite facet indices
179 static int HexPrism_F [8][7] = { // FORWARD
180 { 0, 1, 2, 3, 4, 5, 0 },
181 { 6,11,10, 9, 8, 7, 6 },
182 { 0, 6, 7, 1, 0, 0, 0 },
183 { 1, 7, 8, 2, 1, 1, 1 },
184 { 2, 8, 9, 3, 2, 2, 2 },
185 { 3, 9,10, 4, 3, 3, 3 },
186 { 4,10,11, 5, 4, 4, 4 },
187 { 5,11, 6, 0, 5, 5, 5 }};
188 static int HexPrism_RE [8][7] = { // REVERSED -> EXTERNAL
189 { 0, 5, 4, 3, 2, 1, 0 },
190 { 6,11,10, 9, 8, 7, 6 },
191 { 0, 6, 7, 1, 0, 0, 0 },
192 { 1, 7, 8, 2, 1, 1, 1 },
193 { 2, 8, 9, 3, 2, 2, 2 },
194 { 3, 9,10, 4, 3, 3, 3 },
195 { 4,10,11, 5, 4, 4, 4 },
196 { 5,11, 6, 0, 5, 5, 5 }};
197 static int HexPrism_nbN [] = { 6, 6, 4, 4, 4, 4, 4, 4 };
206 // N0 +---|---+ N1 TETRAHEDRON
214 static int QuadTetra_F [4][7] = { // FORWARD
215 { 0, 4, 1, 5, 2, 6, 0 }, // All faces have external normals
216 { 0, 7, 3, 8, 1, 4, 0 },
217 { 1, 8, 3, 9, 2, 5, 1 },
218 { 0, 6, 2, 9, 3, 7, 0 }};
219 static int QuadTetra_RE [4][7] = { // REVERSED -> FORWARD (EXTERNAL)
220 { 0, 6, 2, 5, 1, 4, 0 }, // All faces have external normals
221 { 0, 4, 1, 8, 3, 7, 0 },
222 { 1, 5, 2, 9, 3, 8, 1 },
223 { 0, 7, 3, 9, 2, 6, 0 }};
224 static int QuadTetra_nbN [] = { 6, 6, 6, 6 };
234 // | | 9 - middle point for (0,4) etc.
247 static int QuadPyram_F [5][9] = { // FORWARD
248 { 0, 5, 1, 6, 2, 7, 3, 8, 0 }, // All faces have external normals
249 { 0, 9, 4, 10,1, 5, 0, 4, 4 },
250 { 1, 10,4, 11,2, 6, 1, 4, 4 },
251 { 2, 11,4, 12,3, 7, 2, 4, 4 },
252 { 3, 12,4, 9, 0, 8, 3, 4, 4 }};
253 static int QuadPyram_RE [5][9] = { // REVERSED -> FORWARD (EXTERNAL)
254 { 0, 8, 3, 7, 2, 6, 1, 5, 0 }, // All faces but a bottom have external normals
255 { 0, 5, 1, 10,4, 9, 0, 4, 4 },
256 { 1, 6, 2, 11,4, 10,1, 4, 4 },
257 { 2, 7, 3, 12,4, 11,2, 4, 4 },
258 { 3, 8, 0, 9, 4, 12,3, 4, 4 }};
259 static int QuadPyram_nbN [] = { 8, 6, 6, 6, 6 };
283 static int QuadPenta_F [5][9] = { // FORWARD
284 { 0, 6, 1, 7, 2, 8, 0, 0, 0 },
285 { 3, 11,5, 10,4, 9, 3, 3, 3 },
286 { 0, 12,3, 9, 4, 13,1, 6, 0 },
287 { 1, 13,4, 10,5, 14,2, 7, 1 },
288 { 0, 8, 2, 14,5, 11,3, 12,0 }};
289 static int QuadPenta_RE [5][9] = { // REVERSED -> EXTERNAL
290 { 0, 8, 2, 7, 1, 6, 0, 0, 0 },
291 { 3, 9, 4, 10,5, 11,3, 3, 3 },
292 { 0, 6, 1, 13,4, 9, 3, 12,0 },
293 { 1, 7, 2, 14,5, 10,4, 13,1 },
294 { 0, 12,3, 11,5, 14,2, 8, 0 }};
295 static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 };
299 // N5+-----+-----+N6 +-----+-----+
301 // 12+ | 14+ | + | +25 + |
303 // N4+-----+-----+N7 | QUADRATIC +-----+-----+ | Central nodes
304 // | | 15 | | HEXAHEDRON | | | | of tri-quadratic
305 // | | | | | | | | HEXAHEDRON
306 // | 17+ | +18 | + 22+ | +
308 // | | | | | + | 26+ | + |
310 // 16+ | +19 | + | +24 + |
313 // | N1+-----+-|---+N2 | +-----+-|---+
315 // | +8 | +10 | + 20+ | +
317 // N0+-----+-----+N3 +-----+-----+
320 static int QuadHexa_F [6][9] = { // FORWARD
321 { 0, 8, 1, 9, 2, 10,3, 11,0 }, // all face normals are external,
322 { 4, 15,7, 14,6, 13,5, 12,4 },
323 { 0, 16,4, 12,5, 17,1, 8, 0 },
324 { 1, 17,5, 13,6, 18,2, 9, 1 },
325 { 3, 10,2, 18,6, 14,7, 19,3 },
326 { 0, 11,3, 19,7, 15,4, 16,0 }};
327 static int QuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL
328 { 0, 11,3, 10,2, 9, 1, 8, 0 }, // all face normals are external
329 { 4, 12,5, 13,6, 14,7, 15,4 },
330 { 0, 8, 1, 17,5, 12,4, 16,0 },
331 { 1, 9, 2, 18,6, 13,5, 17,1 },
332 { 3, 19,7, 14,6, 18,2, 10,3 },
333 { 0, 16,4, 15,7, 19,3, 11,0 }};
334 static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 };
336 static int TriQuadHexa_F [6][9] = { // FORWARD
337 { 0, 8, 1, 9, 2, 10,3, 11, 20 }, // all face normals are external
338 { 4, 15,7, 14,6, 13,5, 12, 25 },
339 { 0, 16,4, 12,5, 17,1, 8, 21 },
340 { 1, 17,5, 13,6, 18,2, 9, 22 },
341 { 3, 10,2, 18,6, 14,7, 19, 23 },
342 { 0, 11,3, 19,7, 15,4, 16, 24 }};
343 static int TriQuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL
344 { 0, 11,3, 10,2, 9, 1, 8, 20 }, // opposite faces are neighbouring,
345 { 4, 12,5, 13,6, 14,7, 15, 25 }, // all face normals are external
346 { 0, 8, 1, 17,5, 12,4, 16, 21 },
347 { 1, 9, 2, 18,6, 13,5, 17, 22 },
348 { 3, 19,7, 14,6, 18,2, 10, 23 },
349 { 0, 16,4, 15,7, 19,3, 11, 24 }};
350 static int TriQuadHexa_nbN [] = { 9, 9, 9, 9, 9, 9 };
353 // ========================================================
354 // to perform some calculations without linkage to CASCADE
355 // ========================================================
360 XYZ() { x = 0; y = 0; z = 0; }
361 XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
362 XYZ( const XYZ& other ) { x = other.x; y = other.y; z = other.z; }
363 XYZ( const SMDS_MeshNode* n ) { x = n->X(); y = n->Y(); z = n->Z(); }
364 inline XYZ operator-( const XYZ& other );
365 inline XYZ operator+( const XYZ& other );
366 inline XYZ Crossed( const XYZ& other );
367 inline double Dot( const XYZ& other );
368 inline double Magnitude();
369 inline double SquareMagnitude();
371 inline XYZ XYZ::operator-( const XYZ& Right ) {
372 return XYZ(x - Right.x, y - Right.y, z - Right.z);
374 inline XYZ XYZ::operator+( const XYZ& Right ) {
375 return XYZ(x + Right.x, y + Right.y, z + Right.z);
377 inline XYZ XYZ::Crossed( const XYZ& Right ) {
378 return XYZ (y * Right.z - z * Right.y,
379 z * Right.x - x * Right.z,
380 x * Right.y - y * Right.x);
382 inline double XYZ::Dot( const XYZ& Other ) {
383 return(x * Other.x + y * Other.y + z * Other.z);
385 inline double XYZ::Magnitude() {
386 return sqrt (x * x + y * y + z * z);
388 inline double XYZ::SquareMagnitude() {
389 return (x * x + y * y + z * z);
392 //================================================================================
394 * \brief Return linear type corresponding to a quadratic one
396 //================================================================================
398 SMDS_VolumeTool::VolumeType quadToLinear(SMDS_VolumeTool::VolumeType quadType)
400 SMDS_VolumeTool::VolumeType linType = SMDS_VolumeTool::VolumeType( int(quadType)-4 );
401 const int nbCornersByQuad = SMDS_VolumeTool::NbCornerNodes( quadType );
402 if ( SMDS_VolumeTool::NbCornerNodes( linType ) == nbCornersByQuad )
406 for ( ; iLin < SMDS_VolumeTool::NB_VOLUME_TYPES; ++iLin )
407 if ( SMDS_VolumeTool::NbCornerNodes( SMDS_VolumeTool::VolumeType( iLin )) == nbCornersByQuad)
408 return SMDS_VolumeTool::VolumeType( iLin );
410 return SMDS_VolumeTool::UNKNOWN;
415 //================================================================================
417 * \brief Saver/restorer of a SMDS_VolumeTool::myCurFace
419 //================================================================================
421 struct SMDS_VolumeTool::SaveFacet
423 SMDS_VolumeTool::Facet mySaved;
424 SMDS_VolumeTool::Facet& myToRestore;
425 SaveFacet( SMDS_VolumeTool::Facet& facet ): myToRestore( facet )
428 mySaved.myNodes.swap( facet.myNodes );
432 if ( myToRestore.myIndex != mySaved.myIndex )
433 myToRestore = mySaved;
434 myToRestore.myNodes.swap( mySaved.myNodes );
438 //=======================================================================
439 //function : SMDS_VolumeTool
441 //=======================================================================
443 SMDS_VolumeTool::SMDS_VolumeTool ()
448 //=======================================================================
449 //function : SMDS_VolumeTool
451 //=======================================================================
453 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume,
454 const bool ignoreCentralNodes)
456 Set( theVolume, ignoreCentralNodes );
459 //=======================================================================
460 //function : SMDS_VolumeTool
462 //=======================================================================
464 SMDS_VolumeTool::~SMDS_VolumeTool()
466 myCurFace.myNodeIndices = NULL;
469 //=======================================================================
470 //function : SetVolume
471 //purpose : Set volume to iterate on
472 //=======================================================================
474 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume,
475 const bool ignoreCentralNodes,
476 const std::vector<const SMDS_MeshNode*>* otherNodes)
481 myIgnoreCentralNodes = ignoreCentralNodes;
485 myVolumeNodes.clear();
486 myPolyIndices.clear();
487 myPolyQuantities.clear();
488 myPolyFacetOri.clear();
491 myExternalFaces = false;
493 myAllFacesNodeIndices_F = 0;
494 myAllFacesNodeIndices_RE = 0;
495 myAllFacesNbNodes = 0;
497 myCurFace.myIndex = -1;
498 myCurFace.myNodeIndices = NULL;
499 myCurFace.myNodes.clear();
502 if ( !theVolume || theVolume->GetType() != SMDSAbs_Volume )
505 myVolume = theVolume;
506 myNbFaces = theVolume->NbFaces();
507 if ( myVolume->IsPoly() )
509 myPolyedre = dynamic_cast<const SMDS_VtkVolume*>( myVolume );
510 myPolyFacetOri.resize( myNbFaces, 0 );
514 myVolumeNodes.resize( myVolume->NbNodes() );
517 if ( otherNodes->size() != myVolumeNodes.size() )
518 return ( myVolume = 0 );
519 for ( size_t i = 0; i < otherNodes->size(); ++i )
520 if ( ! ( myVolumeNodes[i] = (*otherNodes)[0] ))
521 return ( myVolume = 0 );
526 SMDS_ElemIteratorPtr nodeIt = myVolume->nodesIterator();
527 while ( nodeIt->more() )
528 myVolumeNodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
533 return ( myVolume = 0 );
537 // define volume orientation
539 if ( GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z ))
541 const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
542 int topNodeIndex = myVolume->NbCornerNodes() - 1;
543 while ( !IsLinked( 0, topNodeIndex, /*ignoreMediumNodes=*/true )) --topNodeIndex;
544 const SMDS_MeshNode* topNode = myVolumeNodes[ topNodeIndex ];
545 XYZ upDir (topNode->X() - botNode->X(),
546 topNode->Y() - botNode->Y(),
547 topNode->Z() - botNode->Z() );
548 myVolForward = ( botNormal.Dot( upDir ) < 0 );
551 myCurFace.myIndex = -1; // previous setFace(0) didn't take myVolForward into account
556 //=======================================================================
558 //purpose : Inverse volume
559 //=======================================================================
561 #define SWAP_NODES(nodes,i1,i2) \
563 const SMDS_MeshNode* tmp = nodes[ i1 ]; \
564 nodes[ i1 ] = nodes[ i2 ]; \
567 void SMDS_VolumeTool::Inverse ()
569 if ( !myVolume ) return;
571 if (myVolume->IsPoly()) {
572 MESSAGE("Warning: attempt to inverse polyhedral volume");
576 myVolForward = !myVolForward;
577 myCurFace.myIndex = -1;
579 // inverse top and bottom faces
580 switch ( myVolumeNodes.size() ) {
582 SWAP_NODES( myVolumeNodes, 1, 2 );
585 SWAP_NODES( myVolumeNodes, 1, 3 );
588 SWAP_NODES( myVolumeNodes, 1, 2 );
589 SWAP_NODES( myVolumeNodes, 4, 5 );
592 SWAP_NODES( myVolumeNodes, 1, 3 );
593 SWAP_NODES( myVolumeNodes, 5, 7 );
596 SWAP_NODES( myVolumeNodes, 1, 5 );
597 SWAP_NODES( myVolumeNodes, 2, 4 );
598 SWAP_NODES( myVolumeNodes, 7, 11 );
599 SWAP_NODES( myVolumeNodes, 8, 10 );
603 SWAP_NODES( myVolumeNodes, 1, 2 );
604 SWAP_NODES( myVolumeNodes, 4, 6 );
605 SWAP_NODES( myVolumeNodes, 8, 9 );
608 SWAP_NODES( myVolumeNodes, 1, 3 );
609 SWAP_NODES( myVolumeNodes, 5, 8 );
610 SWAP_NODES( myVolumeNodes, 6, 7 );
611 SWAP_NODES( myVolumeNodes, 10, 12 );
614 SWAP_NODES( myVolumeNodes, 1, 2 );
615 SWAP_NODES( myVolumeNodes, 4, 5 );
616 SWAP_NODES( myVolumeNodes, 6, 8 );
617 SWAP_NODES( myVolumeNodes, 9, 11 );
618 SWAP_NODES( myVolumeNodes, 13, 14 );
621 SWAP_NODES( myVolumeNodes, 1, 3 );
622 SWAP_NODES( myVolumeNodes, 5, 7 );
623 SWAP_NODES( myVolumeNodes, 8, 11 );
624 SWAP_NODES( myVolumeNodes, 9, 10 );
625 SWAP_NODES( myVolumeNodes, 12, 15 );
626 SWAP_NODES( myVolumeNodes, 13, 14 );
627 SWAP_NODES( myVolumeNodes, 17, 19 );
630 SWAP_NODES( myVolumeNodes, 1, 3 );
631 SWAP_NODES( myVolumeNodes, 5, 7 );
632 SWAP_NODES( myVolumeNodes, 8, 11 );
633 SWAP_NODES( myVolumeNodes, 9, 10 );
634 SWAP_NODES( myVolumeNodes, 12, 15 );
635 SWAP_NODES( myVolumeNodes, 13, 14 );
636 SWAP_NODES( myVolumeNodes, 17, 19 );
637 SWAP_NODES( myVolumeNodes, 21, 24 );
638 SWAP_NODES( myVolumeNodes, 22, 23 );
644 //=======================================================================
645 //function : GetVolumeType
647 //=======================================================================
649 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
654 switch( myVolumeNodes.size() ) {
655 case 4: return TETRA;
656 case 5: return PYRAM;
657 case 6: return PENTA;
659 case 12: return HEX_PRISM;
660 case 10: return QUAD_TETRA;
661 case 13: return QUAD_PYRAM;
662 case 15: return QUAD_PENTA;
663 case 20: return QUAD_HEXA;
664 case 27: return QUAD_HEXA;
671 //=======================================================================
672 //function : getTetraVolume
674 //=======================================================================
676 static double getTetraVolume(const SMDS_MeshNode* n1,
677 const SMDS_MeshNode* n2,
678 const SMDS_MeshNode* n3,
679 const SMDS_MeshNode* n4)
681 double p1[3], p2[3], p3[3], p4[3];
687 double Q1 = -(p1[ 0 ]-p2[ 0 ])*(p3[ 1 ]*p4[ 2 ]-p4[ 1 ]*p3[ 2 ]);
688 double Q2 = (p1[ 0 ]-p3[ 0 ])*(p2[ 1 ]*p4[ 2 ]-p4[ 1 ]*p2[ 2 ]);
689 double R1 = -(p1[ 0 ]-p4[ 0 ])*(p2[ 1 ]*p3[ 2 ]-p3[ 1 ]*p2[ 2 ]);
690 double R2 = -(p2[ 0 ]-p3[ 0 ])*(p1[ 1 ]*p4[ 2 ]-p4[ 1 ]*p1[ 2 ]);
691 double S1 = (p2[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p3[ 2 ]-p3[ 1 ]*p1[ 2 ]);
692 double S2 = -(p3[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p2[ 2 ]-p2[ 1 ]*p1[ 2 ]);
694 return (Q1+Q2+R1+R2+S1+S2)/6.0;
697 //=======================================================================
699 //purpose : Return element volume
700 //=======================================================================
702 double SMDS_VolumeTool::GetSize() const
708 if ( myVolume->IsPoly() )
713 // split a polyhedron into tetrahedrons
715 SaveFacet savedFacet( myCurFace );
716 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
717 for ( int f = 0; f < NbFaces(); ++f )
720 XYZ area (0,0,0), p1( myCurFace.myNodes[0] );
721 for ( int n = 0; n < myCurFace.myNbNodes; ++n )
723 XYZ p2( myCurFace.myNodes[ n+1 ]);
724 area = area + p1.Crossed( p2 );
733 const static int ind[] = {
734 0, 1, 3, 6, 11, 23, 31, 44, 58, 78 };
735 const static int vtab[][4] = { // decomposition into tetra in the order of enum VolumeType
768 // quadratic tetrahedron
793 // quadratic pentahedron
810 // quadratic hexahedron
835 int type = GetVolumeType();
837 int n2 = ind[type+1];
839 for (int i = n1; i < n2; i++) {
840 V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
841 myVolumeNodes[ vtab[i][1] ],
842 myVolumeNodes[ vtab[i][2] ],
843 myVolumeNodes[ vtab[i][3] ]);
849 //=======================================================================
850 //function : GetBaryCenter
852 //=======================================================================
854 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
860 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
861 X += myVolumeNodes[ i ]->X();
862 Y += myVolumeNodes[ i ]->Y();
863 Z += myVolumeNodes[ i ]->Z();
865 X /= myVolumeNodes.size();
866 Y /= myVolumeNodes.size();
867 Z /= myVolumeNodes.size();
872 //================================================================================
874 * \brief Classify a point
875 * \param tol - thickness of faces
877 //================================================================================
879 bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
881 // LIMITATION: for convex volumes only
883 for ( int iF = 0; iF < myNbFaces; ++iF )
886 if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z ))
888 if ( !IsFaceExternal( iF ))
889 faceNormal = XYZ() - faceNormal; // reverse
891 XYZ face2p( p - XYZ( myCurFace.myNodes[0] ));
892 if ( face2p.Dot( faceNormal ) > tol )
898 //=======================================================================
899 //function : SetExternalNormal
900 //purpose : Node order will be so that faces normals are external
901 //=======================================================================
903 void SMDS_VolumeTool::SetExternalNormal ()
905 myExternalFaces = true;
906 myCurFace.myIndex = -1;
909 //=======================================================================
910 //function : NbFaceNodes
911 //purpose : Return number of nodes in the array of face nodes
912 //=======================================================================
914 int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const
916 if ( !setFace( faceIndex ))
918 return myCurFace.myNbNodes;
921 //=======================================================================
922 //function : GetFaceNodes
923 //purpose : Return pointer to the array of face nodes.
924 // To comfort link iteration, the array
925 // length == NbFaceNodes( faceIndex ) + 1 and
926 // the last node == the first one.
927 //=======================================================================
929 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const
931 if ( !setFace( faceIndex ))
933 return &myCurFace.myNodes[0];
936 //=======================================================================
937 //function : GetFaceNodesIndices
938 //purpose : Return pointer to the array of face nodes indices
939 // To comfort link iteration, the array
940 // length == NbFaceNodes( faceIndex ) + 1 and
941 // the last node index == the first one.
942 //=======================================================================
944 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const
946 if ( !setFace( faceIndex ))
949 return myCurFace.myNodeIndices;
952 //=======================================================================
953 //function : GetFaceNodes
954 //purpose : Return a set of face nodes.
955 //=======================================================================
957 bool SMDS_VolumeTool::GetFaceNodes (int faceIndex,
958 set<const SMDS_MeshNode*>& theFaceNodes ) const
960 if ( !setFace( faceIndex ))
963 theFaceNodes.clear();
964 theFaceNodes.insert( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
971 struct NLink : public std::pair<int,int>
974 NLink(const SMDS_MeshNode* n1=0, const SMDS_MeshNode* n2=0, int ori=1 )
978 if (( myOri = ( n1->GetID() < n2->GetID() )))
981 second = n2->GetID();
987 second = n1->GetID();
993 myOri = first = second = 0;
996 //int Node1() const { return myOri == -1 ? second : first; }
998 //bool IsSameOri( const std::pair<int,int>& link ) const { return link.first == Node1(); }
1002 //=======================================================================
1003 //function : IsFaceExternal
1004 //purpose : Check normal orientation of a given face
1005 //=======================================================================
1007 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const
1009 if ( myExternalFaces || !myVolume )
1012 if ( !myPolyedre ) // all classical volumes have external facet normals
1015 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1017 if ( myPolyFacetOri[ faceIndex ])
1018 return myPolyFacetOri[ faceIndex ] > 0;
1020 int ori = 0; // -1-in, +1-out, 0-undef
1021 double minProj, maxProj;
1022 if ( projectNodesToNormal( faceIndex, minProj, maxProj ))
1024 // all nodes are on the same side of the facet
1025 ori = ( minProj < 0 ? +1 : -1 );
1026 me->myPolyFacetOri[ faceIndex ] = ori;
1028 if ( !myFwdLinks.empty() ) // concave polyhedron; collect oriented links
1029 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1031 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1032 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1037 SaveFacet savedFacet( myCurFace );
1039 // concave polyhedron
1041 if ( myFwdLinks.empty() ) // get links of the least ambiguously oriented facet
1043 for ( size_t i = 0; i < myPolyFacetOri.size() && !ori; ++i )
1044 ori = myPolyFacetOri[ i ];
1046 if ( !ori ) // none facet is oriented yet
1048 // find the least ambiguously oriented facet
1049 int faceMostConvex = -1;
1050 std::map< double, int > convexity2face;
1051 for ( size_t iF = 0; iF < myPolyFacetOri.size() && faceMostConvex < 0; ++iF )
1053 if ( projectNodesToNormal( iF, minProj, maxProj ))
1055 // all nodes are on the same side of the facet
1056 me->myPolyFacetOri[ iF ] = ( minProj < 0 ? +1 : -1 );
1057 faceMostConvex = iF;
1061 ori = ( -minProj < maxProj ? -1 : +1 );
1062 double convexity = std::min( -minProj, maxProj ) / std::max( -minProj, maxProj );
1063 convexity2face.insert( make_pair( convexity, iF * ori ));
1066 if ( faceMostConvex < 0 ) // none facet has nodes on the same side
1068 // use the least ambiguous facet
1069 faceMostConvex = convexity2face.begin()->second;
1070 ori = ( faceMostConvex < 0 ? -1 : +1 );
1071 faceMostConvex = std::abs( faceMostConvex );
1072 me->myPolyFacetOri[ faceMostConvex ] = ori;
1075 // collect links of the oriented facets in myFwdLinks
1076 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1078 ori = myPolyFacetOri[ iF ];
1079 if ( !ori ) continue;
1081 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1083 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1084 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1089 // compare orientation of links of the facet with myFwdLinks
1091 setFace( faceIndex );
1092 vector< NLink > links( myCurFace.myNbNodes ), links2;
1093 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1095 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1096 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1097 if ( l2o != myFwdLinks.end() )
1098 ori = link.myOri * l2o->second * -1;
1101 while ( !ori ) // the facet has no common links with already oriented facets
1103 // orient and collect links of other non-oriented facets
1104 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1106 if ( myPolyFacetOri[ iF ] ) continue; // already oriented
1110 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1112 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1113 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1114 if ( l2o != myFwdLinks.end() )
1115 ori = link.myOri * l2o->second * -1;
1116 links2.push_back( link );
1118 if ( ori ) // one more facet oriented
1120 me->myPolyFacetOri[ iF ] = ori;
1121 for ( size_t i = 0; i < links2.size(); ++i )
1122 me->myFwdLinks.insert( make_pair( links2[i], links2[i].myOri * ori ));
1127 return false; // error in algorithm: infinite loop
1129 // try to orient the facet again
1131 for ( size_t i = 0; i < links.size() && !ori; ++i )
1133 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( links[i] );
1134 if ( l2o != myFwdLinks.end() )
1135 ori = links[i].myOri * l2o->second * -1;
1137 me->myPolyFacetOri[ faceIndex ] = ori;
1143 //=======================================================================
1144 //function : projectNodesToNormal
1145 //purpose : compute min and max projections of all nodes to normal of a facet.
1146 //=======================================================================
1148 bool SMDS_VolumeTool::projectNodesToNormal( int faceIndex,
1150 double& maxProj ) const
1152 minProj = std::numeric_limits<double>::max();
1153 maxProj = std::numeric_limits<double>::min();
1156 if ( !GetFaceNormal( faceIndex, normal.x, normal.y, normal.z ))
1158 XYZ p0 ( myCurFace.myNodes[0] );
1159 for ( size_t i = 0; i < myVolumeNodes.size(); ++i )
1161 if ( std::find( myCurFace.myNodes.begin() + 1,
1162 myCurFace.myNodes.end(),
1163 myVolumeNodes[ i ] ) != myCurFace.myNodes.end() )
1164 continue; // node of the faceIndex-th facet
1166 double proj = normal.Dot( XYZ( myVolumeNodes[ i ]) - p0 );
1167 if ( proj < minProj ) minProj = proj;
1168 if ( proj > maxProj ) maxProj = proj;
1170 const double tol = 1e-7;
1173 bool diffSize = ( minProj * maxProj < 0 );
1176 // minProj = -minProj;
1178 // else if ( minProj < 0 )
1180 // minProj = -minProj;
1181 // maxProj = -maxProj;
1184 return !diffSize; // ? 0 : (minProj >= 0);
1187 //=======================================================================
1188 //function : GetFaceNormal
1189 //purpose : Return a normal to a face
1190 //=======================================================================
1192 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) const
1194 if ( !setFace( faceIndex ))
1197 const int iQuad = ( !myPolyedre && myCurFace.myNbNodes > 6 ) ? 2 : 1;
1198 XYZ p1 ( myCurFace.myNodes[0*iQuad] );
1199 XYZ p2 ( myCurFace.myNodes[1*iQuad] );
1200 XYZ p3 ( myCurFace.myNodes[2*iQuad] );
1201 XYZ aVec12( p2 - p1 );
1202 XYZ aVec13( p3 - p1 );
1203 XYZ cross = aVec12.Crossed( aVec13 );
1205 if ( myCurFace.myNbNodes >3*iQuad ) {
1206 XYZ p4 ( myCurFace.myNodes[3*iQuad] );
1207 XYZ aVec14( p4 - p1 );
1208 XYZ cross2 = aVec13.Crossed( aVec14 );
1209 cross = cross + cross2;
1212 double size = cross.Magnitude();
1213 if ( size <= numeric_limits<double>::min() )
1223 //================================================================================
1225 * \brief Return barycenter of a face
1227 //================================================================================
1229 bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const
1231 if ( !setFace( faceIndex ))
1235 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1237 X += myCurFace.myNodes[i]->X() / myCurFace.myNbNodes;
1238 Y += myCurFace.myNodes[i]->Y() / myCurFace.myNbNodes;
1239 Z += myCurFace.myNodes[i]->Z() / myCurFace.myNbNodes;
1244 //=======================================================================
1245 //function : GetFaceArea
1246 //purpose : Return face area
1247 //=======================================================================
1249 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
1252 if ( !setFace( faceIndex ))
1255 XYZ p1 ( myCurFace.myNodes[0] );
1256 XYZ p2 ( myCurFace.myNodes[1] );
1257 XYZ p3 ( myCurFace.myNodes[2] );
1258 XYZ aVec12( p2 - p1 );
1259 XYZ aVec13( p3 - p1 );
1260 area += aVec12.Crossed( aVec13 ).Magnitude();
1262 if (myVolume->IsPoly())
1264 for ( int i = 3; i < myCurFace.myNbNodes; ++i )
1266 XYZ pI ( myCurFace.myNodes[i] );
1267 XYZ aVecI( pI - p1 );
1268 area += aVec13.Crossed( aVecI ).Magnitude();
1274 if ( myCurFace.myNbNodes == 4 ) {
1275 XYZ p4 ( myCurFace.myNodes[3] );
1276 XYZ aVec14( p4 - p1 );
1277 area += aVec14.Crossed( aVec13 ).Magnitude();
1283 //================================================================================
1285 * \brief Return index of the node located at face center of a quadratic element like HEX27
1287 //================================================================================
1289 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
1291 if ( myAllFacesNbNodes && myVolumeNodes.size() == 27 ) // classic element with 27 nodes
1293 switch ( faceIndex ) {
1297 return faceIndex + 19;
1303 //=======================================================================
1304 //function : GetOppFaceIndex
1305 //purpose : Return index of the opposite face if it exists, else -1.
1306 //=======================================================================
1308 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1312 MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1316 const int nbHoriFaces = 2;
1318 if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1319 switch ( myVolumeNodes.size() ) {
1322 if ( faceIndex == 0 || faceIndex == 1 )
1323 ind = 1 - faceIndex;
1327 if ( faceIndex <= 1 ) // top or bottom
1328 ind = 1 - faceIndex;
1330 const int nbSideFaces = myAllFacesNbNodes[0];
1331 ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1336 ind = GetOppFaceIndexOfHex( faceIndex );
1344 //=======================================================================
1345 //function : GetOppFaceIndexOfHex
1346 //purpose : Return index of the opposite face of the hexahedron
1347 //=======================================================================
1349 int SMDS_VolumeTool::GetOppFaceIndexOfHex( int faceIndex )
1351 return Hexa_oppF[ faceIndex ];
1354 //=======================================================================
1355 //function : IsLinked
1356 //purpose : return true if theNode1 is linked with theNode2
1357 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1358 //=======================================================================
1360 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1361 const SMDS_MeshNode* theNode2,
1362 const bool theIgnoreMediumNodes) const
1367 if (myVolume->IsPoly()) {
1369 MESSAGE("Warning: bad volumic element");
1372 if ( !myAllFacesNbNodes ) {
1373 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1374 me->myPolyQuantities = myPolyedre->GetQuantities();
1375 myAllFacesNbNodes = &myPolyQuantities[0];
1377 int from, to = 0, d1 = 1, d2 = 2;
1378 if ( myPolyedre->IsQuadratic() ) {
1379 if ( theIgnoreMediumNodes ) {
1385 vector<const SMDS_MeshNode*>::const_iterator i;
1386 for (int iface = 0; iface < myNbFaces; iface++)
1389 to += myPolyQuantities[iface];
1390 i = std::find( myVolumeNodes.begin() + from, myVolumeNodes.begin() + to, theNode1 );
1391 if ( i != myVolumeNodes.end() )
1393 if (( theNode2 == *( i-d1 ) ||
1394 theNode2 == *( i+d1 )))
1397 (( theNode2 == *( i-d2 ) ||
1398 theNode2 == *( i+d2 ))))
1405 // find nodes indices
1406 int i1 = -1, i2 = -1, nbFound = 0;
1407 for ( size_t i = 0; i < myVolumeNodes.size() && nbFound < 2; i++ )
1409 if ( myVolumeNodes[ i ] == theNode1 )
1411 else if ( myVolumeNodes[ i ] == theNode2 )
1414 return IsLinked( i1, i2 );
1417 //=======================================================================
1418 //function : IsLinked
1419 //purpose : return true if the node with theNode1Index is linked
1420 // with the node with theNode2Index
1421 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1422 //=======================================================================
1424 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1425 const int theNode2Index,
1426 bool theIgnoreMediumNodes) const
1428 if ( myVolume->IsPoly() ) {
1429 return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1432 int minInd = min( theNode1Index, theNode2Index );
1433 int maxInd = max( theNode1Index, theNode2Index );
1435 if ( minInd < 0 || maxInd > (int)myVolumeNodes.size() - 1 || maxInd == minInd )
1438 VolumeType type = GetVolumeType();
1439 if ( myVolume->IsQuadratic() )
1441 int firstMediumInd = myVolume->NbCornerNodes();
1442 if ( minInd >= firstMediumInd )
1443 return false; // both nodes are medium - not linked
1444 if ( maxInd < firstMediumInd ) // both nodes are corners
1446 if ( theIgnoreMediumNodes )
1447 type = quadToLinear(type); // to check linkage of corner nodes only
1449 return false; // corner nodes are not linked directly in a quadratic cell
1457 switch ( maxInd - minInd ) {
1458 case 1: return minInd != 3;
1459 case 3: return minInd == 0 || minInd == 4;
1460 case 4: return true;
1467 switch ( maxInd - minInd ) {
1469 case 3: return true;
1474 switch ( maxInd - minInd ) {
1475 case 1: return minInd != 2;
1476 case 2: return minInd == 0 || minInd == 3;
1477 case 3: return true;
1484 case 0: if( maxInd==4 || maxInd==6 || maxInd==7 ) return true;
1485 case 1: if( maxInd==4 || maxInd==5 || maxInd==8 ) return true;
1486 case 2: if( maxInd==5 || maxInd==6 || maxInd==9 ) return true;
1487 case 3: if( maxInd==7 || maxInd==8 || maxInd==9 ) return true;
1495 case 0: if( maxInd==8 || maxInd==11 || maxInd==16 ) return true;
1496 case 1: if( maxInd==8 || maxInd==9 || maxInd==17 ) return true;
1497 case 2: if( maxInd==9 || maxInd==10 || maxInd==18 ) return true;
1498 case 3: if( maxInd==10 || maxInd==11 || maxInd==19 ) return true;
1499 case 4: if( maxInd==12 || maxInd==15 || maxInd==16 ) return true;
1500 case 5: if( maxInd==12 || maxInd==13 || maxInd==17 ) return true;
1501 case 6: if( maxInd==13 || maxInd==14 || maxInd==18 ) return true;
1502 case 7: if( maxInd==14 || maxInd==15 || maxInd==19 ) return true;
1510 case 0: if( maxInd==5 || maxInd==8 || maxInd==9 ) return true;
1511 case 1: if( maxInd==5 || maxInd==6 || maxInd==10 ) return true;
1512 case 2: if( maxInd==6 || maxInd==7 || maxInd==11 ) return true;
1513 case 3: if( maxInd==7 || maxInd==8 || maxInd==12 ) return true;
1514 case 4: if( maxInd==9 || maxInd==10 || maxInd==11 || maxInd==12 ) return true;
1522 case 0: if( maxInd==6 || maxInd==8 || maxInd==12 ) return true;
1523 case 1: if( maxInd==6 || maxInd==7 || maxInd==13 ) return true;
1524 case 2: if( maxInd==7 || maxInd==8 || maxInd==14 ) return true;
1525 case 3: if( maxInd==9 || maxInd==11 || maxInd==12 ) return true;
1526 case 4: if( maxInd==9 || maxInd==10 || maxInd==13 ) return true;
1527 case 5: if( maxInd==10 || maxInd==11 || maxInd==14 ) return true;
1534 const int diff = maxInd-minInd;
1535 if ( diff > 6 ) return false;// not linked top and bottom
1536 if ( diff == 6 ) return true; // linked top and bottom
1537 return diff == 1 || diff == 7;
1544 //=======================================================================
1545 //function : GetNodeIndex
1546 //purpose : Return an index of theNode
1547 //=======================================================================
1549 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1552 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
1553 if ( myVolumeNodes[ i ] == theNode )
1560 //================================================================================
1562 * \brief Fill vector with boundary faces existing in the mesh
1563 * \param faces - vector of found nodes
1564 * \retval int - nb of found faces
1566 //================================================================================
1568 int SMDS_VolumeTool::GetAllExistingFaces(vector<const SMDS_MeshElement*> & faces) const
1571 SaveFacet savedFacet( myCurFace );
1573 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1575 if ( const SMDS_MeshElement* face = SMDS_Mesh::FindFace( myCurFace.myNodes ))
1576 faces.push_back( face );
1579 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1580 const SMDS_MeshFace* face = 0;
1581 const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1582 switch ( NbFaceNodes( iF )) {
1584 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1586 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1588 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1589 nodes[3], nodes[4], nodes[5]); break;
1591 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1592 nodes[4], nodes[5], nodes[6], nodes[7]); break;
1595 faces.push_back( face );
1597 return faces.size();
1601 //================================================================================
1603 * \brief Fill vector with boundary edges existing in the mesh
1604 * \param edges - vector of found edges
1605 * \retval int - nb of found faces
1607 //================================================================================
1609 int SMDS_VolumeTool::GetAllExistingEdges(vector<const SMDS_MeshElement*> & edges) const
1612 edges.reserve( myVolumeNodes.size() * 2 );
1613 for ( size_t i = 0; i < myVolumeNodes.size()-1; ++i ) {
1614 for ( size_t j = i + 1; j < myVolumeNodes.size(); ++j ) {
1615 if ( IsLinked( i, j )) {
1616 const SMDS_MeshElement* edge =
1617 SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1619 edges.push_back( edge );
1623 return edges.size();
1626 //================================================================================
1628 * \brief Return minimal square distance between connected corner nodes
1630 //================================================================================
1632 double SMDS_VolumeTool::MinLinearSize2() const
1634 double minSize = 1e+100;
1635 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1637 SaveFacet savedFacet( myCurFace );
1639 // it seems that compute distance twice is faster than organization of a sole computing
1640 myCurFace.myIndex = -1;
1641 for ( int iF = 0; iF < myNbFaces; ++iF )
1644 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1646 XYZ n1( myCurFace.myNodes[ iN ]);
1647 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1648 minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
1655 //================================================================================
1657 * \brief Return maximal square distance between connected corner nodes
1659 //================================================================================
1661 double SMDS_VolumeTool::MaxLinearSize2() const
1663 double maxSize = -1e+100;
1664 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1666 SaveFacet savedFacet( myCurFace );
1668 // it seems that compute distance twice is faster than organization of a sole computing
1669 myCurFace.myIndex = -1;
1670 for ( int iF = 0; iF < myNbFaces; ++iF )
1673 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1675 XYZ n1( myCurFace.myNodes[ iN ]);
1676 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1677 maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude());
1684 //================================================================================
1686 * \brief fast check that only one volume is build on the face nodes
1687 * This check is valid for conformal meshes only
1689 //================================================================================
1691 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1693 const bool isFree = true;
1695 if ( !setFace( faceIndex ))
1698 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1700 const int di = myVolume->IsQuadratic() ? 2 : 1;
1701 const int nbN = ( myCurFace.myNbNodes/di <= 4 && !IsPoly()) ? 3 : myCurFace.myNbNodes/di; // nb nodes to check
1703 SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
1704 while ( eIt->more() )
1706 const SMDS_MeshElement* vol = eIt->next();
1707 if ( vol == myVolume )
1710 for ( iN = 1; iN < nbN; ++iN )
1711 if ( vol->GetNodeIndex( nodes[ iN*di ]) < 0 )
1713 if ( iN == nbN ) // nbN nodes are shared with vol
1715 // if ( vol->IsPoly() || vol->NbFaces() > 6 ) // vol is polyhed or hex prism
1717 // int nb = myCurFace.myNbNodes;
1718 // if ( myVolume->GetEntityType() != vol->GetEntityType() )
1719 // nb -= ( GetCenterNodeIndex(0) > 0 );
1720 // set<const SMDS_MeshNode*> faceNodes( nodes, nodes + nb );
1721 // if ( SMDS_VolumeTool( vol ).GetFaceIndex( faceNodes ) < 0 )
1724 if ( otherVol ) *otherVol = vol;
1728 if ( otherVol ) *otherVol = 0;
1732 //================================================================================
1734 * \brief Thorough check that only one volume is build on the face nodes
1736 //================================================================================
1738 bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1740 const bool isFree = true;
1742 if (!setFace( faceIndex ))
1745 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1746 const int nbFaceNodes = myCurFace.myNbNodes;
1748 // evaluate nb of face nodes shared by other volumes
1749 int maxNbShared = -1;
1750 typedef map< const SMDS_MeshElement*, int > TElemIntMap;
1751 TElemIntMap volNbShared;
1752 TElemIntMap::iterator vNbIt;
1753 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1754 const SMDS_MeshNode* n = nodes[ iNode ];
1755 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1756 while ( eIt->more() ) {
1757 const SMDS_MeshElement* elem = eIt->next();
1758 if ( elem != myVolume ) {
1759 vNbIt = volNbShared.insert( make_pair( elem, 0 )).first;
1761 if ( vNbIt->second > maxNbShared )
1762 maxNbShared = vNbIt->second;
1766 if ( maxNbShared < 3 )
1767 return isFree; // is free
1769 // find volumes laying on the opposite side of the face
1770 // and sharing all nodes
1771 XYZ intNormal; // internal normal
1772 GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
1773 if ( IsFaceExternal( faceIndex ))
1774 intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
1775 XYZ p0 ( nodes[0] ), baryCenter;
1776 for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); ) {
1777 const int& nbShared = (*vNbIt).second;
1778 if ( nbShared >= 3 ) {
1779 SMDS_VolumeTool volume( (*vNbIt).first );
1780 volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
1781 XYZ intNormal2( baryCenter - p0 );
1782 if ( intNormal.Dot( intNormal2 ) < 0 ) {
1784 if ( nbShared >= nbFaceNodes )
1786 // a volume shares the whole facet
1787 if ( otherVol ) *otherVol = vNbIt->first;
1794 // remove a volume from volNbShared map
1795 volNbShared.erase( vNbIt++ );
1798 // here volNbShared contains only volumes laying on the opposite side of
1799 // the face and sharing 3 or more but not all face nodes with myVolume
1800 if ( volNbShared.size() < 2 ) {
1801 return isFree; // is free
1804 // check if the whole area of a face is shared
1805 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1807 const SMDS_MeshNode* n = nodes[ iNode ];
1808 // check if n is shared by one of volumes of volNbShared
1809 bool isShared = false;
1810 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1811 while ( eIt->more() && !isShared )
1812 isShared = volNbShared.count( eIt->next() );
1816 if ( otherVol ) *otherVol = volNbShared.begin()->first;
1819 // if ( !myVolume->IsPoly() )
1821 // bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
1822 // for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1823 // SMDS_VolumeTool volume( (*vNbIt).first );
1824 // bool prevLinkShared = false;
1825 // int nbSharedLinks = 0;
1826 // for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1827 // bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
1828 // if ( linkShared )
1830 // if ( linkShared && prevLinkShared &&
1831 // volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
1832 // isShared[ iNode ] = true;
1833 // prevLinkShared = linkShared;
1835 // if ( nbSharedLinks == nbFaceNodes )
1836 // return !free; // is not free
1837 // if ( nbFaceNodes == 4 ) {
1838 // // check traingle parts 1 & 3
1839 // if ( isShared[1] && isShared[3] )
1840 // return !free; // is not free
1841 // // check triangle parts 0 & 2;
1842 // // 0 part could not be checked in the loop; check it here
1843 // if ( isShared[2] && prevLinkShared &&
1844 // volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
1845 // volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
1846 // return !free; // is not free
1853 //=======================================================================
1854 //function : GetFaceIndex
1855 //purpose : Return index of a face formed by theFaceNodes
1856 //=======================================================================
1858 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes,
1859 const int theFaceIndexHint ) const
1861 if ( theFaceIndexHint >= 0 )
1863 int nbNodes = NbFaceNodes( theFaceIndexHint );
1864 if ( nbNodes == (int) theFaceNodes.size() )
1866 const SMDS_MeshNode** nodes = GetFaceNodes( theFaceIndexHint );
1868 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1873 return theFaceIndexHint;
1876 for ( int iFace = 0; iFace < myNbFaces; iFace++ )
1878 if ( iFace == theFaceIndexHint )
1880 int nbNodes = NbFaceNodes( iFace );
1881 if ( nbNodes == (int) theFaceNodes.size() )
1883 const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
1885 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1896 //=======================================================================
1897 //function : GetFaceIndex
1898 //purpose : Return index of a face formed by theFaceNodes
1899 //=======================================================================
1901 /*int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
1903 for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1904 const int* nodes = GetFaceNodesIndices( iFace );
1905 int nbFaceNodes = NbFaceNodes( iFace );
1907 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1908 nodeSet.insert( nodes[ iNode ] );
1909 if ( theFaceNodesIndices == nodeSet )
1915 //=======================================================================
1916 //function : setFace
1918 //=======================================================================
1920 bool SMDS_VolumeTool::setFace( int faceIndex ) const
1925 if ( myCurFace.myIndex == faceIndex )
1928 myCurFace.myIndex = -1;
1930 if ( faceIndex < 0 || faceIndex >= NbFaces() )
1933 if (myVolume->IsPoly())
1936 MESSAGE("Warning: bad volumic element");
1941 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1942 if ( !myAllFacesNbNodes ) {
1943 me->myPolyQuantities = myPolyedre->GetQuantities();
1944 myAllFacesNbNodes = &myPolyQuantities[0];
1946 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
1947 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
1948 me->myPolyIndices.resize( myCurFace.myNbNodes + 1 );
1949 myCurFace.myNodeIndices = & me->myPolyIndices[0];
1950 int shift = std::accumulate( myAllFacesNbNodes, myAllFacesNbNodes+faceIndex, 0 );
1951 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
1953 myCurFace.myNodes [ iNode ] = myVolumeNodes[ shift + iNode ];
1954 myCurFace.myNodeIndices[ iNode ] = shift + iNode;
1956 myCurFace.myNodes [ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ]; // last = first
1957 myCurFace.myNodeIndices[ myCurFace.myNbNodes ] = myCurFace.myNodeIndices[ 0 ];
1959 // check orientation
1960 if (myExternalFaces)
1962 myCurFace.myIndex = faceIndex; // avoid infinite recursion in IsFaceExternal()
1963 myExternalFaces = false; // force normal computation by IsFaceExternal()
1964 if ( !IsFaceExternal( faceIndex ))
1965 std::reverse( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
1966 myExternalFaces = true;
1971 if ( !myAllFacesNodeIndices_F )
1973 // choose data for an element type
1974 switch ( myVolumeNodes.size() ) {
1976 myAllFacesNodeIndices_F = &Tetra_F [0][0];
1977 //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
1978 myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
1979 myAllFacesNbNodes = Tetra_nbN;
1980 myMaxFaceNbNodes = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
1983 myAllFacesNodeIndices_F = &Pyramid_F [0][0];
1984 //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
1985 myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
1986 myAllFacesNbNodes = Pyramid_nbN;
1987 myMaxFaceNbNodes = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
1990 myAllFacesNodeIndices_F = &Penta_F [0][0];
1991 //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
1992 myAllFacesNodeIndices_RE = &Penta_RE[0][0];
1993 myAllFacesNbNodes = Penta_nbN;
1994 myMaxFaceNbNodes = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
1997 myAllFacesNodeIndices_F = &Hexa_F [0][0];
1998 ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
1999 myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
2000 myAllFacesNbNodes = Hexa_nbN;
2001 myMaxFaceNbNodes = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
2004 myAllFacesNodeIndices_F = &QuadTetra_F [0][0];
2005 //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
2006 myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
2007 myAllFacesNbNodes = QuadTetra_nbN;
2008 myMaxFaceNbNodes = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
2011 myAllFacesNodeIndices_F = &QuadPyram_F [0][0];
2012 //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
2013 myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
2014 myAllFacesNbNodes = QuadPyram_nbN;
2015 myMaxFaceNbNodes = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
2018 myAllFacesNodeIndices_F = &QuadPenta_F [0][0];
2019 //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
2020 myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
2021 myAllFacesNbNodes = QuadPenta_nbN;
2022 myMaxFaceNbNodes = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
2026 myAllFacesNodeIndices_F = &QuadHexa_F [0][0];
2027 //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
2028 myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
2029 myAllFacesNbNodes = QuadHexa_nbN;
2030 myMaxFaceNbNodes = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
2031 if ( !myIgnoreCentralNodes && myVolumeNodes.size() == 27 )
2033 myAllFacesNodeIndices_F = &TriQuadHexa_F [0][0];
2034 //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
2035 myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
2036 myAllFacesNbNodes = TriQuadHexa_nbN;
2037 myMaxFaceNbNodes = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
2041 myAllFacesNodeIndices_F = &HexPrism_F [0][0];
2042 //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
2043 myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
2044 myAllFacesNbNodes = HexPrism_nbN;
2045 myMaxFaceNbNodes = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
2051 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2052 // if ( myExternalFaces )
2053 // myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2055 // myCurFace.myNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
2056 myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2059 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2060 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2061 myCurFace.myNodes[ iNode ] = myVolumeNodes[ myCurFace.myNodeIndices[ iNode ]];
2062 myCurFace.myNodes[ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ];
2065 myCurFace.myIndex = faceIndex;
2070 //=======================================================================
2071 //function : GetType
2072 //purpose : return VolumeType by nb of nodes in a volume
2073 //=======================================================================
2075 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
2077 switch ( nbNodes ) {
2078 case 4: return TETRA;
2079 case 5: return PYRAM;
2080 case 6: return PENTA;
2081 case 8: return HEXA;
2082 case 10: return QUAD_TETRA;
2083 case 13: return QUAD_PYRAM;
2084 case 15: return QUAD_PENTA;
2086 case 27: return QUAD_HEXA;
2087 case 12: return HEX_PRISM;
2088 default:return UNKNOWN;
2092 //=======================================================================
2093 //function : NbFaces
2094 //purpose : return nb of faces by volume type
2095 //=======================================================================
2097 int SMDS_VolumeTool::NbFaces( VolumeType type )
2101 case QUAD_TETRA: return 4;
2103 case QUAD_PYRAM: return 5;
2105 case QUAD_PENTA: return 5;
2107 case QUAD_HEXA : return 6;
2108 case HEX_PRISM : return 8;
2113 //================================================================================
2115 * \brief Useful to know nb of corner nodes of a quadratic volume
2116 * \param type - volume type
2117 * \retval int - nb of corner nodes
2119 //================================================================================
2121 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
2125 case QUAD_TETRA: return 4;
2127 case QUAD_PYRAM: return 5;
2129 case QUAD_PENTA: return 6;
2131 case QUAD_HEXA : return 8;
2132 case HEX_PRISM : return 12;
2139 //=======================================================================
2140 //function : GetFaceNodesIndices
2141 //purpose : Return the array of face nodes indices
2142 // To comfort link iteration, the array
2143 // length == NbFaceNodes( faceIndex ) + 1 and
2144 // the last node index == the first one.
2145 //=======================================================================
2147 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
2152 case TETRA: return Tetra_F[ faceIndex ];
2153 case PYRAM: return Pyramid_F[ faceIndex ];
2154 case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
2155 case HEXA: return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
2156 case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
2157 case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
2158 case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
2159 // what about SMDSEntity_TriQuad_Hexa?
2160 case QUAD_HEXA: return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
2161 case HEX_PRISM: return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
2167 //=======================================================================
2168 //function : NbFaceNodes
2169 //purpose : Return number of nodes in the array of face nodes
2170 //=======================================================================
2172 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
2176 case TETRA: return Tetra_nbN[ faceIndex ];
2177 case PYRAM: return Pyramid_nbN[ faceIndex ];
2178 case PENTA: return Penta_nbN[ faceIndex ];
2179 case HEXA: return Hexa_nbN[ faceIndex ];
2180 case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
2181 case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
2182 case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
2183 // what about SMDSEntity_TriQuad_Hexa?
2184 case QUAD_HEXA: return QuadHexa_nbN[ faceIndex ];
2185 case HEX_PRISM: return HexPrism_nbN[ faceIndex ];
2191 //=======================================================================
2192 //function : Element
2193 //purpose : return element
2194 //=======================================================================
2196 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
2198 return static_cast<const SMDS_MeshVolume*>( myVolume );
2201 //=======================================================================
2203 //purpose : return element ID
2204 //=======================================================================
2206 int SMDS_VolumeTool::ID() const
2208 return myVolume ? myVolume->GetID() : 0;