1 // Copyright (C) 2007-2019 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_Mesh.hxx"
37 #include <utilities.h>
47 // ======================================================
48 // Node indices in faces depending on volume orientation
49 // making most faces normals external
50 // ======================================================
51 // For all elements, 0-th face is bottom based on the first nodes.
52 // For prismatic elements (tetra,hexa,prisms), 1-th face is a top one.
53 // For all elements, side faces follow order of bottom nodes
54 // ======================================================
62 // N0 +---|---+ N1 TETRAHEDRON
70 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
71 { 0, 1, 2, 0 }, // All faces have external normals
75 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
76 { 0, 2, 1, 0 }, // All faces have external normals
80 static int Tetra_nbN [] = { 3, 3, 3, 3 };
85 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
86 { 0, 1, 2, 3, 0 }, // All faces have external normals
92 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
93 { 0, 3, 2, 1, 0 }, // All faces but a bottom have external normals
98 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
109 // | / \ | PENTAHEDRON
115 static int Penta_F [5][5] = { // FORWARD
116 { 0, 1, 2, 0, 0 }, // All faces have external normals
117 { 3, 5, 4, 3, 3 }, // 0 is bottom, 1 is top face
121 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
127 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
134 // N4+----------+N7 |
135 // | | | | HEXAHEDRON
136 // | N1+------|---+N2
142 static int Hexa_F [6][5] = { // FORWARD
144 { 4, 7, 6, 5, 4 }, // all face normals are external
149 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
151 { 4, 5, 6, 7, 4 }, // all face normals are external
156 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
157 static int Hexa_oppF[] = { 1, 0, 4, 5, 2, 3 }; // oppopsite facet indices
176 static int HexPrism_F [8][7] = { // FORWARD
177 { 0, 1, 2, 3, 4, 5, 0 },
178 { 6,11,10, 9, 8, 7, 6 },
179 { 0, 6, 7, 1, 0, 0, 0 },
180 { 1, 7, 8, 2, 1, 1, 1 },
181 { 2, 8, 9, 3, 2, 2, 2 },
182 { 3, 9,10, 4, 3, 3, 3 },
183 { 4,10,11, 5, 4, 4, 4 },
184 { 5,11, 6, 0, 5, 5, 5 }};
185 static int HexPrism_RE [8][7] = { // REVERSED -> EXTERNAL
186 { 0, 5, 4, 3, 2, 1, 0 },
187 { 6,11,10, 9, 8, 7, 6 },
188 { 0, 6, 7, 1, 0, 0, 0 },
189 { 1, 7, 8, 2, 1, 1, 1 },
190 { 2, 8, 9, 3, 2, 2, 2 },
191 { 3, 9,10, 4, 3, 3, 3 },
192 { 4,10,11, 5, 4, 4, 4 },
193 { 5,11, 6, 0, 5, 5, 5 }};
194 static int HexPrism_nbN [] = { 6, 6, 4, 4, 4, 4, 4, 4 };
203 // N0 +---|---+ N1 TETRAHEDRON
211 static int QuadTetra_F [4][7] = { // FORWARD
212 { 0, 4, 1, 5, 2, 6, 0 }, // All faces have external normals
213 { 0, 7, 3, 8, 1, 4, 0 },
214 { 1, 8, 3, 9, 2, 5, 1 },
215 { 0, 6, 2, 9, 3, 7, 0 }};
216 static int QuadTetra_RE [4][7] = { // REVERSED -> FORWARD (EXTERNAL)
217 { 0, 6, 2, 5, 1, 4, 0 }, // All faces have external normals
218 { 0, 4, 1, 8, 3, 7, 0 },
219 { 1, 5, 2, 9, 3, 8, 1 },
220 { 0, 7, 3, 9, 2, 6, 0 }};
221 static int QuadTetra_nbN [] = { 6, 6, 6, 6 };
231 // | | 9 - middle point for (0,4) etc.
244 static int QuadPyram_F [5][9] = { // FORWARD
245 { 0, 5, 1, 6, 2, 7, 3, 8, 0 }, // All faces have external normals
246 { 0, 9, 4, 10,1, 5, 0, 4, 4 },
247 { 1, 10,4, 11,2, 6, 1, 4, 4 },
248 { 2, 11,4, 12,3, 7, 2, 4, 4 },
249 { 3, 12,4, 9, 0, 8, 3, 4, 4 }};
250 static int QuadPyram_RE [5][9] = { // REVERSED -> FORWARD (EXTERNAL)
251 { 0, 8, 3, 7, 2, 6, 1, 5, 0 }, // All faces but a bottom have external normals
252 { 0, 5, 1, 10,4, 9, 0, 4, 4 },
253 { 1, 6, 2, 11,4, 10,1, 4, 4 },
254 { 2, 7, 3, 12,4, 11,2, 4, 4 },
255 { 3, 8, 0, 9, 4, 12,3, 4, 4 }};
256 static int QuadPyram_nbN [] = { 8, 6, 6, 6, 6 };
280 static int QuadPenta_F [5][9] = { // FORWARD
281 { 0, 6, 1, 7, 2, 8, 0, 0, 0 },
282 { 3, 11,5, 10,4, 9, 3, 3, 3 },
283 { 0, 12,3, 9, 4, 13,1, 6, 0 },
284 { 1, 13,4, 10,5, 14,2, 7, 1 },
285 { 0, 8, 2, 14,5, 11,3, 12,0 }};
286 static int QuadPenta_RE [5][9] = { // REVERSED -> EXTERNAL
287 { 0, 8, 2, 7, 1, 6, 0, 0, 0 },
288 { 3, 9, 4, 10,5, 11,3, 3, 3 },
289 { 0, 6, 1, 13,4, 9, 3, 12,0 },
290 { 1, 7, 2, 14,5, 10,4, 13,1 },
291 { 0, 12,3, 11,5, 14,2, 8, 0 }};
292 static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 };
296 // N5+-----+-----+N6 +-----+-----+
298 // 12+ | 14+ | + | +25 + |
300 // N4+-----+-----+N7 | QUADRATIC +-----+-----+ | Central nodes
301 // | | 15 | | HEXAHEDRON | | | | of tri-quadratic
302 // | | | | | | | | HEXAHEDRON
303 // | 17+ | +18 | + 22+ | +
305 // | | | | | + | 26+ | + |
307 // 16+ | +19 | + | +24 + |
310 // | N1+-----+-|---+N2 | +-----+-|---+
312 // | +8 | +10 | + 20+ | +
314 // N0+-----+-----+N3 +-----+-----+
317 static int QuadHexa_F [6][9] = { // FORWARD
318 { 0, 8, 1, 9, 2, 10,3, 11,0 }, // all face normals are external,
319 { 4, 15,7, 14,6, 13,5, 12,4 },
320 { 0, 16,4, 12,5, 17,1, 8, 0 },
321 { 1, 17,5, 13,6, 18,2, 9, 1 },
322 { 3, 10,2, 18,6, 14,7, 19,3 },
323 { 0, 11,3, 19,7, 15,4, 16,0 }};
324 static int QuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL
325 { 0, 11,3, 10,2, 9, 1, 8, 0 }, // all face normals are external
326 { 4, 12,5, 13,6, 14,7, 15,4 },
327 { 0, 8, 1, 17,5, 12,4, 16,0 },
328 { 1, 9, 2, 18,6, 13,5, 17,1 },
329 { 3, 19,7, 14,6, 18,2, 10,3 },
330 { 0, 16,4, 15,7, 19,3, 11,0 }};
331 static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 };
333 static int TriQuadHexa_F [6][9] = { // FORWARD
334 { 0, 8, 1, 9, 2, 10,3, 11, 20 }, // all face normals are external
335 { 4, 15,7, 14,6, 13,5, 12, 25 },
336 { 0, 16,4, 12,5, 17,1, 8, 21 },
337 { 1, 17,5, 13,6, 18,2, 9, 22 },
338 { 3, 10,2, 18,6, 14,7, 19, 23 },
339 { 0, 11,3, 19,7, 15,4, 16, 24 }};
340 static int TriQuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL
341 { 0, 11,3, 10,2, 9, 1, 8, 20 }, // opposite faces are neighbouring,
342 { 4, 12,5, 13,6, 14,7, 15, 25 }, // all face normals are external
343 { 0, 8, 1, 17,5, 12,4, 16, 21 },
344 { 1, 9, 2, 18,6, 13,5, 17, 22 },
345 { 3, 19,7, 14,6, 18,2, 10, 23 },
346 { 0, 16,4, 15,7, 19,3, 11, 24 }};
347 static int TriQuadHexa_nbN [] = { 9, 9, 9, 9, 9, 9 };
350 // ========================================================
351 // to perform some calculations without linkage to CASCADE
352 // ========================================================
357 XYZ() { x = 0; y = 0; z = 0; }
358 XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
359 XYZ( const XYZ& other ) { x = other.x; y = other.y; z = other.z; }
360 XYZ( const SMDS_MeshNode* n ) { x = n->X(); y = n->Y(); z = n->Z(); }
361 inline XYZ operator-( const XYZ& other );
362 inline XYZ operator+( const XYZ& other );
363 inline XYZ Crossed( const XYZ& other );
364 inline double Dot( const XYZ& other );
365 inline double Magnitude();
366 inline double SquareMagnitude();
368 inline XYZ XYZ::operator-( const XYZ& Right ) {
369 return XYZ(x - Right.x, y - Right.y, z - Right.z);
371 inline XYZ XYZ::operator+( const XYZ& Right ) {
372 return XYZ(x + Right.x, y + Right.y, z + Right.z);
374 inline XYZ XYZ::Crossed( const XYZ& Right ) {
375 return XYZ (y * Right.z - z * Right.y,
376 z * Right.x - x * Right.z,
377 x * Right.y - y * Right.x);
379 inline double XYZ::Dot( const XYZ& Other ) {
380 return(x * Other.x + y * Other.y + z * Other.z);
382 inline double XYZ::Magnitude() {
383 return sqrt (x * x + y * y + z * z);
385 inline double XYZ::SquareMagnitude() {
386 return (x * x + y * y + z * z);
389 //================================================================================
391 * \brief Return linear type corresponding to a quadratic one
393 //================================================================================
395 SMDS_VolumeTool::VolumeType quadToLinear(SMDS_VolumeTool::VolumeType quadType)
397 SMDS_VolumeTool::VolumeType linType = SMDS_VolumeTool::VolumeType( int(quadType)-4 );
398 const int nbCornersByQuad = SMDS_VolumeTool::NbCornerNodes( quadType );
399 if ( SMDS_VolumeTool::NbCornerNodes( linType ) == nbCornersByQuad )
403 for ( ; iLin < SMDS_VolumeTool::NB_VOLUME_TYPES; ++iLin )
404 if ( SMDS_VolumeTool::NbCornerNodes( SMDS_VolumeTool::VolumeType( iLin )) == nbCornersByQuad)
405 return SMDS_VolumeTool::VolumeType( iLin );
407 return SMDS_VolumeTool::UNKNOWN;
412 //================================================================================
414 * \brief Saver/restorer of a SMDS_VolumeTool::myCurFace
416 //================================================================================
418 struct SMDS_VolumeTool::SaveFacet
420 SMDS_VolumeTool::Facet mySaved;
421 SMDS_VolumeTool::Facet& myToRestore;
422 SaveFacet( SMDS_VolumeTool::Facet& facet ): myToRestore( facet )
425 mySaved.myNodes.swap( facet.myNodes );
429 if ( myToRestore.myIndex != mySaved.myIndex )
430 myToRestore = mySaved;
431 myToRestore.myNodes.swap( mySaved.myNodes );
435 //=======================================================================
436 //function : SMDS_VolumeTool
438 //=======================================================================
440 SMDS_VolumeTool::SMDS_VolumeTool ()
445 //=======================================================================
446 //function : SMDS_VolumeTool
448 //=======================================================================
450 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume,
451 const bool ignoreCentralNodes)
453 Set( theVolume, ignoreCentralNodes );
456 //=======================================================================
457 //function : SMDS_VolumeTool
459 //=======================================================================
461 SMDS_VolumeTool::~SMDS_VolumeTool()
463 myCurFace.myNodeIndices = NULL;
466 //=======================================================================
467 //function : SetVolume
468 //purpose : Set volume to iterate on
469 //=======================================================================
471 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume,
472 const bool ignoreCentralNodes,
473 const std::vector<const SMDS_MeshNode*>* otherNodes)
478 myIgnoreCentralNodes = ignoreCentralNodes;
482 myVolumeNodes.clear();
483 myPolyIndices.clear();
484 myPolyQuantities.clear();
485 myPolyFacetOri.clear();
488 myExternalFaces = false;
490 myAllFacesNodeIndices_F = 0;
491 myAllFacesNodeIndices_RE = 0;
492 myAllFacesNbNodes = 0;
494 myCurFace.myIndex = -1;
495 myCurFace.myNodeIndices = NULL;
496 myCurFace.myNodes.clear();
499 if ( !theVolume || theVolume->GetType() != SMDSAbs_Volume )
502 myVolume = theVolume;
503 myNbFaces = theVolume->NbFaces();
504 if ( myVolume->IsPoly() )
506 myPolyedre = SMDS_Mesh::DownCast<SMDS_MeshVolume>( myVolume );
507 myPolyFacetOri.resize( myNbFaces, 0 );
511 myVolumeNodes.resize( myVolume->NbNodes() );
514 if ( otherNodes->size() != myVolumeNodes.size() )
515 return ( myVolume = 0 );
516 for ( size_t i = 0; i < otherNodes->size(); ++i )
517 if ( ! ( myVolumeNodes[i] = (*otherNodes)[0] ))
518 return ( myVolume = 0 );
522 myVolumeNodes.assign( myVolume->begin_nodes(), myVolume->end_nodes() );
527 return ( myVolume = 0 );
531 // define volume orientation
533 if ( GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z ))
535 const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
536 int topNodeIndex = myVolume->NbCornerNodes() - 1;
537 while ( !IsLinked( 0, topNodeIndex, /*ignoreMediumNodes=*/true )) --topNodeIndex;
538 const SMDS_MeshNode* topNode = myVolumeNodes[ topNodeIndex ];
539 XYZ upDir (topNode->X() - botNode->X(),
540 topNode->Y() - botNode->Y(),
541 topNode->Z() - botNode->Z() );
542 myVolForward = ( botNormal.Dot( upDir ) < 0 );
545 myCurFace.myIndex = -1; // previous setFace(0) didn't take myVolForward into account
550 //=======================================================================
552 //purpose : Inverse volume
553 //=======================================================================
555 #define SWAP_NODES(nodes,i1,i2) \
557 const SMDS_MeshNode* tmp = nodes[ i1 ]; \
558 nodes[ i1 ] = nodes[ i2 ]; \
561 void SMDS_VolumeTool::Inverse ()
563 if ( !myVolume ) return;
565 if (myVolume->IsPoly()) {
566 MESSAGE("Warning: attempt to inverse polyhedral volume");
570 myVolForward = !myVolForward;
571 myCurFace.myIndex = -1;
573 // inverse top and bottom faces
574 switch ( myVolumeNodes.size() ) {
576 SWAP_NODES( myVolumeNodes, 1, 2 );
579 SWAP_NODES( myVolumeNodes, 1, 3 );
582 SWAP_NODES( myVolumeNodes, 1, 2 );
583 SWAP_NODES( myVolumeNodes, 4, 5 );
586 SWAP_NODES( myVolumeNodes, 1, 3 );
587 SWAP_NODES( myVolumeNodes, 5, 7 );
590 SWAP_NODES( myVolumeNodes, 1, 5 );
591 SWAP_NODES( myVolumeNodes, 2, 4 );
592 SWAP_NODES( myVolumeNodes, 7, 11 );
593 SWAP_NODES( myVolumeNodes, 8, 10 );
597 SWAP_NODES( myVolumeNodes, 1, 2 );
598 SWAP_NODES( myVolumeNodes, 4, 6 );
599 SWAP_NODES( myVolumeNodes, 8, 9 );
602 SWAP_NODES( myVolumeNodes, 1, 3 );
603 SWAP_NODES( myVolumeNodes, 5, 8 );
604 SWAP_NODES( myVolumeNodes, 6, 7 );
605 SWAP_NODES( myVolumeNodes, 10, 12 );
608 SWAP_NODES( myVolumeNodes, 1, 2 );
609 SWAP_NODES( myVolumeNodes, 4, 5 );
610 SWAP_NODES( myVolumeNodes, 6, 8 );
611 SWAP_NODES( myVolumeNodes, 9, 11 );
612 SWAP_NODES( myVolumeNodes, 13, 14 );
615 SWAP_NODES( myVolumeNodes, 1, 3 );
616 SWAP_NODES( myVolumeNodes, 5, 7 );
617 SWAP_NODES( myVolumeNodes, 8, 11 );
618 SWAP_NODES( myVolumeNodes, 9, 10 );
619 SWAP_NODES( myVolumeNodes, 12, 15 );
620 SWAP_NODES( myVolumeNodes, 13, 14 );
621 SWAP_NODES( myVolumeNodes, 17, 19 );
624 SWAP_NODES( myVolumeNodes, 1, 3 );
625 SWAP_NODES( myVolumeNodes, 5, 7 );
626 SWAP_NODES( myVolumeNodes, 8, 11 );
627 SWAP_NODES( myVolumeNodes, 9, 10 );
628 SWAP_NODES( myVolumeNodes, 12, 15 );
629 SWAP_NODES( myVolumeNodes, 13, 14 );
630 SWAP_NODES( myVolumeNodes, 17, 19 );
631 SWAP_NODES( myVolumeNodes, 21, 24 );
632 SWAP_NODES( myVolumeNodes, 22, 23 );
638 //=======================================================================
639 //function : GetVolumeType
641 //=======================================================================
643 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
648 switch( myVolumeNodes.size() ) {
649 case 4: return TETRA;
650 case 5: return PYRAM;
651 case 6: return PENTA;
653 case 12: return HEX_PRISM;
654 case 10: return QUAD_TETRA;
655 case 13: return QUAD_PYRAM;
656 case 15: return QUAD_PENTA;
657 case 20: return QUAD_HEXA;
658 case 27: return QUAD_HEXA;
665 //=======================================================================
666 //function : getTetraVolume
668 //=======================================================================
670 static double getTetraVolume(const SMDS_MeshNode* n1,
671 const SMDS_MeshNode* n2,
672 const SMDS_MeshNode* n3,
673 const SMDS_MeshNode* n4)
675 double p1[3], p2[3], p3[3], p4[3];
681 double Q1 = -(p1[ 0 ]-p2[ 0 ])*(p3[ 1 ]*p4[ 2 ]-p4[ 1 ]*p3[ 2 ]);
682 double Q2 = (p1[ 0 ]-p3[ 0 ])*(p2[ 1 ]*p4[ 2 ]-p4[ 1 ]*p2[ 2 ]);
683 double R1 = -(p1[ 0 ]-p4[ 0 ])*(p2[ 1 ]*p3[ 2 ]-p3[ 1 ]*p2[ 2 ]);
684 double R2 = -(p2[ 0 ]-p3[ 0 ])*(p1[ 1 ]*p4[ 2 ]-p4[ 1 ]*p1[ 2 ]);
685 double S1 = (p2[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p3[ 2 ]-p3[ 1 ]*p1[ 2 ]);
686 double S2 = -(p3[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p2[ 2 ]-p2[ 1 ]*p1[ 2 ]);
688 return (Q1+Q2+R1+R2+S1+S2)/6.0;
691 //=======================================================================
693 //purpose : Return element volume
694 //=======================================================================
696 double SMDS_VolumeTool::GetSize() const
702 if ( myVolume->IsPoly() )
707 // split a polyhedron into tetrahedrons
709 SaveFacet savedFacet( myCurFace );
710 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
711 for ( int f = 0; f < NbFaces(); ++f )
714 XYZ area (0,0,0), p1( myCurFace.myNodes[0] );
715 for ( int n = 0; n < myCurFace.myNbNodes; ++n )
717 XYZ p2( myCurFace.myNodes[ n+1 ]);
718 area = area + p1.Crossed( p2 );
727 const static int ind[] = {
728 0, 1, 3, 6, 11, 23, 31, 44, 58, 78 };
729 const static int vtab[][4] = { // decomposition into tetra in the order of enum VolumeType
762 // quadratic tetrahedron
787 // quadratic pentahedron
804 // quadratic hexahedron
829 int type = GetVolumeType();
831 int n2 = ind[type+1];
833 for (int i = n1; i < n2; i++) {
834 V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
835 myVolumeNodes[ vtab[i][1] ],
836 myVolumeNodes[ vtab[i][2] ],
837 myVolumeNodes[ vtab[i][3] ]);
843 //=======================================================================
844 //function : GetBaryCenter
846 //=======================================================================
848 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
854 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
855 X += myVolumeNodes[ i ]->X();
856 Y += myVolumeNodes[ i ]->Y();
857 Z += myVolumeNodes[ i ]->Z();
859 X /= myVolumeNodes.size();
860 Y /= myVolumeNodes.size();
861 Z /= myVolumeNodes.size();
866 //================================================================================
868 * \brief Classify a point
869 * \param tol - thickness of faces
871 //================================================================================
873 bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
875 // LIMITATION: for convex volumes only
877 for ( int iF = 0; iF < myNbFaces; ++iF )
880 if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z ))
882 if ( !IsFaceExternal( iF ))
883 faceNormal = XYZ() - faceNormal; // reverse
885 XYZ face2p( p - XYZ( myCurFace.myNodes[0] ));
886 if ( face2p.Dot( faceNormal ) > tol )
892 //=======================================================================
893 //function : SetExternalNormal
894 //purpose : Node order will be so that faces normals are external
895 //=======================================================================
897 void SMDS_VolumeTool::SetExternalNormal ()
899 myExternalFaces = true;
900 myCurFace.myIndex = -1;
903 //=======================================================================
904 //function : NbFaceNodes
905 //purpose : Return number of nodes in the array of face nodes
906 //=======================================================================
908 int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const
910 if ( !setFace( faceIndex ))
912 return myCurFace.myNbNodes;
915 //=======================================================================
916 //function : GetFaceNodes
917 //purpose : Return pointer to the array of face nodes.
918 // To comfort link iteration, the array
919 // length == NbFaceNodes( faceIndex ) + 1 and
920 // the last node == the first one.
921 //=======================================================================
923 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const
925 if ( !setFace( faceIndex ))
927 return &myCurFace.myNodes[0];
930 //=======================================================================
931 //function : GetFaceNodesIndices
932 //purpose : Return pointer to the array of face nodes indices
933 // To comfort link iteration, the array
934 // length == NbFaceNodes( faceIndex ) + 1 and
935 // the last node index == the first one.
936 //=======================================================================
938 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const
940 if ( !setFace( faceIndex ))
943 return myCurFace.myNodeIndices;
946 //=======================================================================
947 //function : GetFaceNodes
948 //purpose : Return a set of face nodes.
949 //=======================================================================
951 bool SMDS_VolumeTool::GetFaceNodes (int faceIndex,
952 std::set<const SMDS_MeshNode*>& theFaceNodes ) const
954 if ( !setFace( faceIndex ))
957 theFaceNodes.clear();
958 theFaceNodes.insert( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
965 struct NLink : public std::pair<int,int>
968 NLink(const SMDS_MeshNode* n1=0, const SMDS_MeshNode* n2=0, int ori=1 )
972 if (( myOri = ( n1->GetID() < n2->GetID() )))
975 second = n2->GetID();
981 second = n1->GetID();
987 myOri = first = second = 0;
990 //int Node1() const { return myOri == -1 ? second : first; }
992 //bool IsSameOri( const std::pair<int,int>& link ) const { return link.first == Node1(); }
996 //=======================================================================
997 //function : IsFaceExternal
998 //purpose : Check normal orientation of a given face
999 //=======================================================================
1001 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const
1003 if ( myExternalFaces || !myVolume )
1006 if ( !myPolyedre ) // all classical volumes have external facet normals
1009 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1011 if ( myPolyFacetOri[ faceIndex ])
1012 return myPolyFacetOri[ faceIndex ] > 0;
1014 int ori = 0; // -1-in, +1-out, 0-undef
1015 double minProj, maxProj;
1016 if ( projectNodesToNormal( faceIndex, minProj, maxProj ))
1018 // all nodes are on the same side of the facet
1019 ori = ( minProj < 0 ? +1 : -1 );
1020 me->myPolyFacetOri[ faceIndex ] = ori;
1022 if ( !myFwdLinks.empty() ) // concave polyhedron; collect oriented links
1023 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1025 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1026 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1031 SaveFacet savedFacet( myCurFace );
1033 // concave polyhedron
1035 if ( myFwdLinks.empty() ) // get links of the least ambiguously oriented facet
1037 for ( size_t i = 0; i < myPolyFacetOri.size() && !ori; ++i )
1038 ori = myPolyFacetOri[ i ];
1040 if ( !ori ) // none facet is oriented yet
1042 // find the least ambiguously oriented facet
1043 int faceMostConvex = -1;
1044 std::map< double, int > convexity2face;
1045 for ( size_t iF = 0; iF < myPolyFacetOri.size() && faceMostConvex < 0; ++iF )
1047 if ( projectNodesToNormal( iF, minProj, maxProj ))
1049 // all nodes are on the same side of the facet
1050 me->myPolyFacetOri[ iF ] = ( minProj < 0 ? +1 : -1 );
1051 faceMostConvex = iF;
1055 ori = ( -minProj < maxProj ? -1 : +1 );
1056 double convexity = std::min( -minProj, maxProj ) / std::max( -minProj, maxProj );
1057 convexity2face.insert( std::make_pair( convexity, iF * ori ));
1060 if ( faceMostConvex < 0 ) // none facet has nodes on the same side
1062 // use the least ambiguous facet
1063 faceMostConvex = convexity2face.begin()->second;
1064 ori = ( faceMostConvex < 0 ? -1 : +1 );
1065 faceMostConvex = std::abs( faceMostConvex );
1066 me->myPolyFacetOri[ faceMostConvex ] = ori;
1069 // collect links of the oriented facets in myFwdLinks
1070 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1072 ori = myPolyFacetOri[ iF ];
1073 if ( !ori ) continue;
1075 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1077 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1078 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1083 // compare orientation of links of the facet with myFwdLinks
1085 setFace( faceIndex );
1086 std::vector< NLink > links( myCurFace.myNbNodes ), links2;
1087 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1089 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1090 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1091 if ( l2o != myFwdLinks.end() )
1092 ori = link.myOri * l2o->second * -1;
1095 while ( !ori ) // the facet has no common links with already oriented facets
1097 // orient and collect links of other non-oriented facets
1098 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1100 if ( myPolyFacetOri[ iF ] ) continue; // already oriented
1104 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1106 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1107 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1108 if ( l2o != myFwdLinks.end() )
1109 ori = link.myOri * l2o->second * -1;
1110 links2.push_back( link );
1112 if ( ori ) // one more facet oriented
1114 me->myPolyFacetOri[ iF ] = ori;
1115 for ( size_t i = 0; i < links2.size(); ++i )
1116 me->myFwdLinks.insert( make_pair( links2[i], links2[i].myOri * ori ));
1121 return false; // error in algorithm: infinite loop
1123 // try to orient the facet again
1125 for ( size_t i = 0; i < links.size() && !ori; ++i )
1127 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( links[i] );
1128 if ( l2o != myFwdLinks.end() )
1129 ori = links[i].myOri * l2o->second * -1;
1131 me->myPolyFacetOri[ faceIndex ] = ori;
1137 //=======================================================================
1138 //function : projectNodesToNormal
1139 //purpose : compute min and max projections of all nodes to normal of a facet.
1140 //=======================================================================
1142 bool SMDS_VolumeTool::projectNodesToNormal( int faceIndex,
1144 double& maxProj ) const
1146 minProj = std::numeric_limits<double>::max();
1147 maxProj = std::numeric_limits<double>::min();
1150 if ( !GetFaceNormal( faceIndex, normal.x, normal.y, normal.z ))
1152 XYZ p0 ( myCurFace.myNodes[0] );
1153 for ( size_t i = 0; i < myVolumeNodes.size(); ++i )
1155 if ( std::find( myCurFace.myNodes.begin() + 1,
1156 myCurFace.myNodes.end(),
1157 myVolumeNodes[ i ] ) != myCurFace.myNodes.end() )
1158 continue; // node of the faceIndex-th facet
1160 double proj = normal.Dot( XYZ( myVolumeNodes[ i ]) - p0 );
1161 if ( proj < minProj ) minProj = proj;
1162 if ( proj > maxProj ) maxProj = proj;
1164 const double tol = 1e-7;
1167 bool diffSize = ( minProj * maxProj < 0 );
1170 // minProj = -minProj;
1172 // else if ( minProj < 0 )
1174 // minProj = -minProj;
1175 // maxProj = -maxProj;
1178 return !diffSize; // ? 0 : (minProj >= 0);
1181 //=======================================================================
1182 //function : GetFaceNormal
1183 //purpose : Return a normal to a face
1184 //=======================================================================
1186 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) const
1188 if ( !setFace( faceIndex ))
1191 const int iQuad = ( !myPolyedre && myCurFace.myNbNodes > 6 ) ? 2 : 1;
1192 XYZ p1 ( myCurFace.myNodes[0*iQuad] );
1193 XYZ p2 ( myCurFace.myNodes[1*iQuad] );
1194 XYZ p3 ( myCurFace.myNodes[2*iQuad] );
1195 XYZ aVec12( p2 - p1 );
1196 XYZ aVec13( p3 - p1 );
1197 XYZ cross = aVec12.Crossed( aVec13 );
1199 if ( myCurFace.myNbNodes >3*iQuad ) {
1200 XYZ p4 ( myCurFace.myNodes[3*iQuad] );
1201 XYZ aVec14( p4 - p1 );
1202 XYZ cross2 = aVec13.Crossed( aVec14 );
1203 cross = cross + cross2;
1206 double size = cross.Magnitude();
1207 if ( size <= std::numeric_limits<double>::min() )
1217 //================================================================================
1219 * \brief Return barycenter of a face
1221 //================================================================================
1223 bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const
1225 if ( !setFace( faceIndex ))
1229 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1231 X += myCurFace.myNodes[i]->X() / myCurFace.myNbNodes;
1232 Y += myCurFace.myNodes[i]->Y() / myCurFace.myNbNodes;
1233 Z += myCurFace.myNodes[i]->Z() / myCurFace.myNbNodes;
1238 //=======================================================================
1239 //function : GetFaceArea
1240 //purpose : Return face area
1241 //=======================================================================
1243 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
1246 if ( !setFace( faceIndex ))
1249 XYZ p1 ( myCurFace.myNodes[0] );
1250 XYZ p2 ( myCurFace.myNodes[1] );
1251 XYZ p3 ( myCurFace.myNodes[2] );
1252 XYZ aVec12( p2 - p1 );
1253 XYZ aVec13( p3 - p1 );
1254 area += aVec12.Crossed( aVec13 ).Magnitude();
1256 if (myVolume->IsPoly())
1258 for ( int i = 3; i < myCurFace.myNbNodes; ++i )
1260 XYZ pI ( myCurFace.myNodes[i] );
1261 XYZ aVecI( pI - p1 );
1262 area += aVec13.Crossed( aVecI ).Magnitude();
1268 if ( myCurFace.myNbNodes == 4 ) {
1269 XYZ p4 ( myCurFace.myNodes[3] );
1270 XYZ aVec14( p4 - p1 );
1271 area += aVec14.Crossed( aVec13 ).Magnitude();
1277 //================================================================================
1279 * \brief Return index of the node located at face center of a quadratic element like HEX27
1281 //================================================================================
1283 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
1285 if ( myAllFacesNbNodes && myVolumeNodes.size() == 27 ) // classic element with 27 nodes
1287 switch ( faceIndex ) {
1291 return faceIndex + 19;
1297 //=======================================================================
1298 //function : GetOppFaceIndex
1299 //purpose : Return index of the opposite face if it exists, else -1.
1300 //=======================================================================
1302 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1306 MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1310 const int nbHoriFaces = 2;
1312 if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1313 switch ( myVolumeNodes.size() ) {
1316 if ( faceIndex == 0 || faceIndex == 1 )
1317 ind = 1 - faceIndex;
1321 if ( faceIndex <= 1 ) // top or bottom
1322 ind = 1 - faceIndex;
1324 const int nbSideFaces = myAllFacesNbNodes[0];
1325 ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1330 ind = GetOppFaceIndexOfHex( faceIndex );
1338 //=======================================================================
1339 //function : GetOppFaceIndexOfHex
1340 //purpose : Return index of the opposite face of the hexahedron
1341 //=======================================================================
1343 int SMDS_VolumeTool::GetOppFaceIndexOfHex( int faceIndex )
1345 return Hexa_oppF[ faceIndex ];
1348 //=======================================================================
1349 //function : IsLinked
1350 //purpose : return true if theNode1 is linked with theNode2
1351 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1352 //=======================================================================
1354 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1355 const SMDS_MeshNode* theNode2,
1356 const bool theIgnoreMediumNodes) const
1361 if (myVolume->IsPoly()) {
1363 MESSAGE("Warning: bad volumic element");
1366 if ( !myAllFacesNbNodes ) {
1367 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1368 me->myPolyQuantities = myPolyedre->GetQuantities();
1369 myAllFacesNbNodes = &myPolyQuantities[0];
1371 int from, to = 0, d1 = 1, d2 = 2;
1372 if ( myPolyedre->IsQuadratic() ) {
1373 if ( theIgnoreMediumNodes ) {
1379 std::vector<const SMDS_MeshNode*>::const_iterator i;
1380 for (int iface = 0; iface < myNbFaces; iface++)
1383 to += myPolyQuantities[iface];
1384 i = std::find( myVolumeNodes.begin() + from, myVolumeNodes.begin() + to, theNode1 );
1385 if ( i != myVolumeNodes.end() )
1387 if (( theNode2 == *( i-d1 ) ||
1388 theNode2 == *( i+d1 )))
1391 (( theNode2 == *( i-d2 ) ||
1392 theNode2 == *( i+d2 ))))
1399 // find nodes indices
1400 int i1 = -1, i2 = -1, nbFound = 0;
1401 for ( size_t i = 0; i < myVolumeNodes.size() && nbFound < 2; i++ )
1403 if ( myVolumeNodes[ i ] == theNode1 )
1405 else if ( myVolumeNodes[ i ] == theNode2 )
1408 return IsLinked( i1, i2 );
1411 //=======================================================================
1412 //function : IsLinked
1413 //purpose : return true if the node with theNode1Index is linked
1414 // with the node with theNode2Index
1415 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1416 //=======================================================================
1418 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1419 const int theNode2Index,
1420 bool theIgnoreMediumNodes) const
1422 if ( myVolume->IsPoly() ) {
1423 return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1426 int minInd = std::min( theNode1Index, theNode2Index );
1427 int maxInd = std::max( theNode1Index, theNode2Index );
1429 if ( minInd < 0 || maxInd > (int)myVolumeNodes.size() - 1 || maxInd == minInd )
1432 VolumeType type = GetVolumeType();
1433 if ( myVolume->IsQuadratic() )
1435 int firstMediumInd = myVolume->NbCornerNodes();
1436 if ( minInd >= firstMediumInd )
1437 return false; // both nodes are medium - not linked
1438 if ( maxInd < firstMediumInd ) // both nodes are corners
1440 if ( theIgnoreMediumNodes )
1441 type = quadToLinear(type); // to check linkage of corner nodes only
1443 return false; // corner nodes are not linked directly in a quadratic cell
1451 switch ( maxInd - minInd ) {
1452 case 1: return minInd != 3;
1453 case 3: return minInd == 0 || minInd == 4;
1454 case 4: return true;
1461 switch ( maxInd - minInd ) {
1463 case 3: return true;
1468 switch ( maxInd - minInd ) {
1469 case 1: return minInd != 2;
1470 case 2: return minInd == 0 || minInd == 3;
1471 case 3: return true;
1478 case 0: if( maxInd==4 || maxInd==6 || maxInd==7 ) return true;
1479 case 1: if( maxInd==4 || maxInd==5 || maxInd==8 ) return true;
1480 case 2: if( maxInd==5 || maxInd==6 || maxInd==9 ) return true;
1481 case 3: if( maxInd==7 || maxInd==8 || maxInd==9 ) return true;
1489 case 0: if( maxInd==8 || maxInd==11 || maxInd==16 ) return true;
1490 case 1: if( maxInd==8 || maxInd==9 || maxInd==17 ) return true;
1491 case 2: if( maxInd==9 || maxInd==10 || maxInd==18 ) return true;
1492 case 3: if( maxInd==10 || maxInd==11 || maxInd==19 ) return true;
1493 case 4: if( maxInd==12 || maxInd==15 || maxInd==16 ) return true;
1494 case 5: if( maxInd==12 || maxInd==13 || maxInd==17 ) return true;
1495 case 6: if( maxInd==13 || maxInd==14 || maxInd==18 ) return true;
1496 case 7: if( maxInd==14 || maxInd==15 || maxInd==19 ) return true;
1504 case 0: if( maxInd==5 || maxInd==8 || maxInd==9 ) return true;
1505 case 1: if( maxInd==5 || maxInd==6 || maxInd==10 ) return true;
1506 case 2: if( maxInd==6 || maxInd==7 || maxInd==11 ) return true;
1507 case 3: if( maxInd==7 || maxInd==8 || maxInd==12 ) return true;
1508 case 4: if( maxInd==9 || maxInd==10 || maxInd==11 || maxInd==12 ) return true;
1516 case 0: if( maxInd==6 || maxInd==8 || maxInd==12 ) return true;
1517 case 1: if( maxInd==6 || maxInd==7 || maxInd==13 ) return true;
1518 case 2: if( maxInd==7 || maxInd==8 || maxInd==14 ) return true;
1519 case 3: if( maxInd==9 || maxInd==11 || maxInd==12 ) return true;
1520 case 4: if( maxInd==9 || maxInd==10 || maxInd==13 ) return true;
1521 case 5: if( maxInd==10 || maxInd==11 || maxInd==14 ) return true;
1528 const int diff = maxInd-minInd;
1529 if ( diff > 6 ) return false;// not linked top and bottom
1530 if ( diff == 6 ) return true; // linked top and bottom
1531 return diff == 1 || diff == 7;
1538 //=======================================================================
1539 //function : GetNodeIndex
1540 //purpose : Return an index of theNode
1541 //=======================================================================
1543 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1546 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
1547 if ( myVolumeNodes[ i ] == theNode )
1554 //================================================================================
1556 * \brief Fill vector with boundary faces existing in the mesh
1557 * \param faces - vector of found nodes
1558 * \retval int - nb of found faces
1560 //================================================================================
1562 int SMDS_VolumeTool::GetAllExistingFaces(std::vector<const SMDS_MeshElement*> & faces) const
1565 SaveFacet savedFacet( myCurFace );
1567 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1569 if ( const SMDS_MeshElement* face = SMDS_Mesh::FindFace( myCurFace.myNodes ))
1570 faces.push_back( face );
1573 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1574 const SMDS_MeshFace* face = 0;
1575 const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1576 switch ( NbFaceNodes( iF )) {
1578 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1580 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1582 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1583 nodes[3], nodes[4], nodes[5]); break;
1585 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1586 nodes[4], nodes[5], nodes[6], nodes[7]); break;
1589 faces.push_back( face );
1591 return faces.size();
1595 //================================================================================
1597 * \brief Fill vector with boundary edges existing in the mesh
1598 * \param edges - vector of found edges
1599 * \retval int - nb of found faces
1601 //================================================================================
1603 int SMDS_VolumeTool::GetAllExistingEdges(std::vector<const SMDS_MeshElement*> & edges) const
1606 edges.reserve( myVolumeNodes.size() * 2 );
1607 for ( size_t i = 0; i < myVolumeNodes.size()-1; ++i ) {
1608 for ( size_t j = i + 1; j < myVolumeNodes.size(); ++j ) {
1609 if ( IsLinked( i, j )) {
1610 const SMDS_MeshElement* edge =
1611 SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1613 edges.push_back( edge );
1617 return edges.size();
1620 //================================================================================
1622 * \brief Return minimal square distance between connected corner nodes
1624 //================================================================================
1626 double SMDS_VolumeTool::MinLinearSize2() const
1628 double minSize = 1e+100;
1629 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1631 SaveFacet savedFacet( myCurFace );
1633 // it seems that compute distance twice is faster than organization of a sole computing
1634 myCurFace.myIndex = -1;
1635 for ( int iF = 0; iF < myNbFaces; ++iF )
1638 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1640 XYZ n1( myCurFace.myNodes[ iN ]);
1641 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1642 minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
1649 //================================================================================
1651 * \brief Return maximal square distance between connected corner nodes
1653 //================================================================================
1655 double SMDS_VolumeTool::MaxLinearSize2() const
1657 double maxSize = -1e+100;
1658 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1660 SaveFacet savedFacet( myCurFace );
1662 // it seems that compute distance twice is faster than organization of a sole computing
1663 myCurFace.myIndex = -1;
1664 for ( int iF = 0; iF < myNbFaces; ++iF )
1667 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1669 XYZ n1( myCurFace.myNodes[ iN ]);
1670 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1671 maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude());
1678 //================================================================================
1680 * \brief Fast quickly check that only one volume is built on the face nodes
1681 * This check is valid for conformal meshes only
1683 //================================================================================
1685 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1687 const bool isFree = true;
1689 if ( !setFace( faceIndex ))
1692 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1694 const int di = myVolume->IsQuadratic() ? 2 : 1;
1695 const int nbN = ( myCurFace.myNbNodes/di <= 4 && !IsPoly()) ? 3 : myCurFace.myNbNodes/di; // nb nodes to check
1697 SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
1698 while ( eIt->more() )
1700 const SMDS_MeshElement* vol = eIt->next();
1701 if ( vol == myVolume )
1704 for ( iN = 1; iN < nbN; ++iN )
1705 if ( vol->GetNodeIndex( nodes[ iN*di ]) < 0 )
1707 if ( iN == nbN ) // nbN nodes are shared with vol
1709 // if ( vol->IsPoly() || vol->NbFaces() > 6 ) // vol is polyhed or hex prism
1711 // int nb = myCurFace.myNbNodes;
1712 // if ( myVolume->GetEntityType() != vol->GetEntityType() )
1713 // nb -= ( GetCenterNodeIndex(0) > 0 );
1714 // std::set<const SMDS_MeshNode*> faceNodes( nodes, nodes + nb );
1715 // if ( SMDS_VolumeTool( vol ).GetFaceIndex( faceNodes ) < 0 )
1718 if ( otherVol ) *otherVol = vol;
1722 if ( otherVol ) *otherVol = 0;
1726 //================================================================================
1728 * \brief Thorough check that only one volume is built on the face nodes
1730 //================================================================================
1732 bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1734 const bool isFree = true;
1736 if (!setFace( faceIndex ))
1739 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1740 const int nbFaceNodes = myCurFace.myNbNodes;
1742 // evaluate nb of face nodes shared by other volumes
1743 int maxNbShared = -1;
1744 typedef std::map< const SMDS_MeshElement*, int > TElemIntMap;
1745 TElemIntMap volNbShared;
1746 TElemIntMap::iterator vNbIt;
1747 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1748 const SMDS_MeshNode* n = nodes[ iNode ];
1749 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1750 while ( eIt->more() ) {
1751 const SMDS_MeshElement* elem = eIt->next();
1752 if ( elem != myVolume ) {
1753 vNbIt = volNbShared.insert( std::make_pair( elem, 0 )).first;
1755 if ( vNbIt->second > maxNbShared )
1756 maxNbShared = vNbIt->second;
1760 if ( maxNbShared < 3 )
1761 return isFree; // is free
1763 // find volumes laying on the opposite side of the face
1764 // and sharing all nodes
1765 XYZ intNormal; // internal normal
1766 GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
1767 if ( IsFaceExternal( faceIndex ))
1768 intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
1769 XYZ p0 ( nodes[0] ), baryCenter;
1770 for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); ) {
1771 const int& nbShared = (*vNbIt).second;
1772 if ( nbShared >= 3 ) {
1773 SMDS_VolumeTool volume( (*vNbIt).first );
1774 volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
1775 XYZ intNormal2( baryCenter - p0 );
1776 if ( intNormal.Dot( intNormal2 ) < 0 ) {
1778 if ( nbShared >= nbFaceNodes )
1780 // a volume shares the whole facet
1781 if ( otherVol ) *otherVol = vNbIt->first;
1788 // remove a volume from volNbShared map
1789 volNbShared.erase( vNbIt++ );
1792 // here volNbShared contains only volumes laying on the opposite side of
1793 // the face and sharing 3 or more but not all face nodes with myVolume
1794 if ( volNbShared.size() < 2 ) {
1795 return isFree; // is free
1798 // check if the whole area of a face is shared
1799 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1801 const SMDS_MeshNode* n = nodes[ iNode ];
1802 // check if n is shared by one of volumes of volNbShared
1803 bool isShared = false;
1804 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1805 while ( eIt->more() && !isShared )
1806 isShared = volNbShared.count( eIt->next() );
1810 if ( otherVol ) *otherVol = volNbShared.begin()->first;
1813 // if ( !myVolume->IsPoly() )
1815 // bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
1816 // for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1817 // SMDS_VolumeTool volume( (*vNbIt).first );
1818 // bool prevLinkShared = false;
1819 // int nbSharedLinks = 0;
1820 // for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1821 // bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
1822 // if ( linkShared )
1824 // if ( linkShared && prevLinkShared &&
1825 // volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
1826 // isShared[ iNode ] = true;
1827 // prevLinkShared = linkShared;
1829 // if ( nbSharedLinks == nbFaceNodes )
1830 // return !free; // is not free
1831 // if ( nbFaceNodes == 4 ) {
1832 // // check traingle parts 1 & 3
1833 // if ( isShared[1] && isShared[3] )
1834 // return !free; // is not free
1835 // // check triangle parts 0 & 2;
1836 // // 0 part could not be checked in the loop; check it here
1837 // if ( isShared[2] && prevLinkShared &&
1838 // volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
1839 // volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
1840 // return !free; // is not free
1847 //=======================================================================
1848 //function : GetFaceIndex
1849 //purpose : Return index of a face formed by theFaceNodes
1850 //=======================================================================
1852 int SMDS_VolumeTool::GetFaceIndex( const std::set<const SMDS_MeshNode*>& theFaceNodes,
1853 const int theFaceIndexHint ) const
1855 if ( theFaceIndexHint >= 0 )
1857 int nbNodes = NbFaceNodes( theFaceIndexHint );
1858 if ( nbNodes == (int) theFaceNodes.size() )
1860 const SMDS_MeshNode** nodes = GetFaceNodes( theFaceIndexHint );
1862 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1867 return theFaceIndexHint;
1870 for ( int iFace = 0; iFace < myNbFaces; iFace++ )
1872 if ( iFace == theFaceIndexHint )
1874 int nbNodes = NbFaceNodes( iFace );
1875 if ( nbNodes == (int) theFaceNodes.size() )
1877 const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
1879 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1890 //=======================================================================
1891 //function : GetFaceIndex
1892 //purpose : Return index of a face formed by theFaceNodes
1893 //=======================================================================
1895 /*int SMDS_VolumeTool::GetFaceIndex( const std::set<int>& theFaceNodesIndices )
1897 for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1898 const int* nodes = GetFaceNodesIndices( iFace );
1899 int nbFaceNodes = NbFaceNodes( iFace );
1900 std::set<int> nodeSet;
1901 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1902 nodeSet.insert( nodes[ iNode ] );
1903 if ( theFaceNodesIndices == nodeSet )
1909 //=======================================================================
1910 //function : setFace
1912 //=======================================================================
1914 bool SMDS_VolumeTool::setFace( int faceIndex ) const
1919 if ( myCurFace.myIndex == faceIndex )
1922 myCurFace.myIndex = -1;
1924 if ( faceIndex < 0 || faceIndex >= NbFaces() )
1927 if (myVolume->IsPoly())
1929 if ( !myPolyedre ) {
1930 MESSAGE("Warning: bad volumic element");
1935 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1936 if ( !myAllFacesNbNodes ) {
1937 me->myPolyQuantities = myPolyedre->GetQuantities();
1938 myAllFacesNbNodes = &myPolyQuantities[0];
1940 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
1941 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
1942 me->myPolyIndices.resize( myCurFace.myNbNodes + 1 );
1943 myCurFace.myNodeIndices = & me->myPolyIndices[0];
1944 int shift = std::accumulate( myAllFacesNbNodes, myAllFacesNbNodes+faceIndex, 0 );
1945 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
1947 myCurFace.myNodes [ iNode ] = myVolumeNodes[ shift + iNode ];
1948 myCurFace.myNodeIndices[ iNode ] = shift + iNode;
1950 myCurFace.myNodes [ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ]; // last = first
1951 myCurFace.myNodeIndices[ myCurFace.myNbNodes ] = myCurFace.myNodeIndices[ 0 ];
1953 // check orientation
1954 if (myExternalFaces)
1956 myCurFace.myIndex = faceIndex; // avoid infinite recursion in IsFaceExternal()
1957 myExternalFaces = false; // force normal computation by IsFaceExternal()
1958 if ( !IsFaceExternal( faceIndex ))
1959 std::reverse( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
1960 myExternalFaces = true;
1965 if ( !myAllFacesNodeIndices_F )
1967 // choose data for an element type
1968 switch ( myVolumeNodes.size() ) {
1970 myAllFacesNodeIndices_F = &Tetra_F [0][0];
1971 //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
1972 myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
1973 myAllFacesNbNodes = Tetra_nbN;
1974 myMaxFaceNbNodes = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
1977 myAllFacesNodeIndices_F = &Pyramid_F [0][0];
1978 //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
1979 myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
1980 myAllFacesNbNodes = Pyramid_nbN;
1981 myMaxFaceNbNodes = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
1984 myAllFacesNodeIndices_F = &Penta_F [0][0];
1985 //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
1986 myAllFacesNodeIndices_RE = &Penta_RE[0][0];
1987 myAllFacesNbNodes = Penta_nbN;
1988 myMaxFaceNbNodes = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
1991 myAllFacesNodeIndices_F = &Hexa_F [0][0];
1992 ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
1993 myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
1994 myAllFacesNbNodes = Hexa_nbN;
1995 myMaxFaceNbNodes = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
1998 myAllFacesNodeIndices_F = &QuadTetra_F [0][0];
1999 //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
2000 myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
2001 myAllFacesNbNodes = QuadTetra_nbN;
2002 myMaxFaceNbNodes = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
2005 myAllFacesNodeIndices_F = &QuadPyram_F [0][0];
2006 //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
2007 myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
2008 myAllFacesNbNodes = QuadPyram_nbN;
2009 myMaxFaceNbNodes = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
2012 myAllFacesNodeIndices_F = &QuadPenta_F [0][0];
2013 //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
2014 myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
2015 myAllFacesNbNodes = QuadPenta_nbN;
2016 myMaxFaceNbNodes = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
2020 myAllFacesNodeIndices_F = &QuadHexa_F [0][0];
2021 //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
2022 myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
2023 myAllFacesNbNodes = QuadHexa_nbN;
2024 myMaxFaceNbNodes = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
2025 if ( !myIgnoreCentralNodes && myVolumeNodes.size() == 27 )
2027 myAllFacesNodeIndices_F = &TriQuadHexa_F [0][0];
2028 //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
2029 myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
2030 myAllFacesNbNodes = TriQuadHexa_nbN;
2031 myMaxFaceNbNodes = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
2035 myAllFacesNodeIndices_F = &HexPrism_F [0][0];
2036 //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
2037 myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
2038 myAllFacesNbNodes = HexPrism_nbN;
2039 myMaxFaceNbNodes = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
2045 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2046 // if ( myExternalFaces )
2047 // myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2049 // myCurFace.myNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
2050 myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2053 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2054 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2055 myCurFace.myNodes[ iNode ] = myVolumeNodes[ myCurFace.myNodeIndices[ iNode ]];
2056 myCurFace.myNodes[ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ];
2059 myCurFace.myIndex = faceIndex;
2064 //=======================================================================
2065 //function : GetType
2066 //purpose : return VolumeType by nb of nodes in a volume
2067 //=======================================================================
2069 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
2071 switch ( nbNodes ) {
2072 case 4: return TETRA;
2073 case 5: return PYRAM;
2074 case 6: return PENTA;
2075 case 8: return HEXA;
2076 case 10: return QUAD_TETRA;
2077 case 13: return QUAD_PYRAM;
2078 case 15: return QUAD_PENTA;
2080 case 27: return QUAD_HEXA;
2081 case 12: return HEX_PRISM;
2082 default:return UNKNOWN;
2086 //=======================================================================
2087 //function : NbFaces
2088 //purpose : return nb of faces by volume type
2089 //=======================================================================
2091 int SMDS_VolumeTool::NbFaces( VolumeType type )
2095 case QUAD_TETRA: return 4;
2097 case QUAD_PYRAM: return 5;
2099 case QUAD_PENTA: return 5;
2101 case QUAD_HEXA : return 6;
2102 case HEX_PRISM : return 8;
2107 //================================================================================
2109 * \brief Useful to know nb of corner nodes of a quadratic volume
2110 * \param type - volume type
2111 * \retval int - nb of corner nodes
2113 //================================================================================
2115 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
2119 case QUAD_TETRA: return 4;
2121 case QUAD_PYRAM: return 5;
2123 case QUAD_PENTA: return 6;
2125 case QUAD_HEXA : return 8;
2126 case HEX_PRISM : return 12;
2133 //=======================================================================
2134 //function : GetFaceNodesIndices
2135 //purpose : Return the array of face nodes indices
2136 // To comfort link iteration, the array
2137 // length == NbFaceNodes( faceIndex ) + 1 and
2138 // the last node index == the first one.
2139 //=======================================================================
2141 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
2146 case TETRA: return Tetra_F[ faceIndex ];
2147 case PYRAM: return Pyramid_F[ faceIndex ];
2148 case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
2149 case HEXA: return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
2150 case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
2151 case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
2152 case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
2153 // what about SMDSEntity_TriQuad_Hexa?
2154 case QUAD_HEXA: return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
2155 case HEX_PRISM: return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
2161 //=======================================================================
2162 //function : NbFaceNodes
2163 //purpose : Return number of nodes in the array of face nodes
2164 //=======================================================================
2166 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
2170 case TETRA: return Tetra_nbN[ faceIndex ];
2171 case PYRAM: return Pyramid_nbN[ faceIndex ];
2172 case PENTA: return Penta_nbN[ faceIndex ];
2173 case HEXA: return Hexa_nbN[ faceIndex ];
2174 case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
2175 case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
2176 case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
2177 // what about SMDSEntity_TriQuad_Hexa?
2178 case QUAD_HEXA: return QuadHexa_nbN[ faceIndex ];
2179 case HEX_PRISM: return HexPrism_nbN[ faceIndex ];
2185 //=======================================================================
2186 //function : Element
2187 //purpose : return element
2188 //=======================================================================
2190 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
2192 return static_cast<const SMDS_MeshVolume*>( myVolume );
2195 //=======================================================================
2197 //purpose : return element ID
2198 //=======================================================================
2200 int SMDS_VolumeTool::ID() const
2202 return myVolume ? myVolume->GetID() : 0;