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 )
431 if ( myToRestore.myIndex != mySaved.myIndex )
432 myToRestore = mySaved;
436 //=======================================================================
437 //function : SMDS_VolumeTool
439 //=======================================================================
441 SMDS_VolumeTool::SMDS_VolumeTool ()
446 //=======================================================================
447 //function : SMDS_VolumeTool
449 //=======================================================================
451 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume,
452 const bool ignoreCentralNodes)
454 Set( theVolume, ignoreCentralNodes );
457 //=======================================================================
458 //function : SMDS_VolumeTool
460 //=======================================================================
462 SMDS_VolumeTool::~SMDS_VolumeTool()
464 myCurFace.myNodeIndices = NULL;
467 //=======================================================================
468 //function : SetVolume
469 //purpose : Set volume to iterate on
470 //=======================================================================
472 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume,
473 const bool ignoreCentralNodes)
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 = dynamic_cast<const SMDS_VtkVolume*>( myVolume );
507 myPolyFacetOri.resize( myNbFaces, 0 );
512 myVolumeNodes.resize( myVolume->NbNodes() );
513 SMDS_ElemIteratorPtr nodeIt = myVolume->nodesIterator();
514 while ( nodeIt->more() )
515 myVolumeNodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
519 return ( myVolume = 0 );
523 // define volume orientation
525 if ( GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z ))
527 const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
528 int topNodeIndex = myVolume->NbCornerNodes() - 1;
529 while ( !IsLinked( 0, topNodeIndex, /*ignoreMediumNodes=*/true )) --topNodeIndex;
530 const SMDS_MeshNode* topNode = myVolumeNodes[ topNodeIndex ];
531 XYZ upDir (topNode->X() - botNode->X(),
532 topNode->Y() - botNode->Y(),
533 topNode->Z() - botNode->Z() );
534 myVolForward = ( botNormal.Dot( upDir ) < 0 );
537 myCurFace.myIndex = -1; // previous setFace(0) didn't take myVolForward into account
542 //=======================================================================
544 //purpose : Inverse volume
545 //=======================================================================
547 #define SWAP_NODES(nodes,i1,i2) \
549 const SMDS_MeshNode* tmp = nodes[ i1 ]; \
550 nodes[ i1 ] = nodes[ i2 ]; \
553 void SMDS_VolumeTool::Inverse ()
555 if ( !myVolume ) return;
557 if (myVolume->IsPoly()) {
558 MESSAGE("Warning: attempt to inverse polyhedral volume");
562 myVolForward = !myVolForward;
563 myCurFace.myIndex = -1;
565 // inverse top and bottom faces
566 switch ( myVolumeNodes.size() ) {
568 SWAP_NODES( myVolumeNodes, 1, 2 );
571 SWAP_NODES( myVolumeNodes, 1, 3 );
574 SWAP_NODES( myVolumeNodes, 1, 2 );
575 SWAP_NODES( myVolumeNodes, 4, 5 );
578 SWAP_NODES( myVolumeNodes, 1, 3 );
579 SWAP_NODES( myVolumeNodes, 5, 7 );
582 SWAP_NODES( myVolumeNodes, 1, 5 );
583 SWAP_NODES( myVolumeNodes, 2, 4 );
584 SWAP_NODES( myVolumeNodes, 7, 11 );
585 SWAP_NODES( myVolumeNodes, 8, 10 );
589 SWAP_NODES( myVolumeNodes, 1, 2 );
590 SWAP_NODES( myVolumeNodes, 4, 6 );
591 SWAP_NODES( myVolumeNodes, 8, 9 );
594 SWAP_NODES( myVolumeNodes, 1, 3 );
595 SWAP_NODES( myVolumeNodes, 5, 8 );
596 SWAP_NODES( myVolumeNodes, 6, 7 );
597 SWAP_NODES( myVolumeNodes, 10, 12 );
600 SWAP_NODES( myVolumeNodes, 1, 2 );
601 SWAP_NODES( myVolumeNodes, 4, 5 );
602 SWAP_NODES( myVolumeNodes, 6, 8 );
603 SWAP_NODES( myVolumeNodes, 9, 11 );
604 SWAP_NODES( myVolumeNodes, 13, 14 );
607 SWAP_NODES( myVolumeNodes, 1, 3 );
608 SWAP_NODES( myVolumeNodes, 5, 7 );
609 SWAP_NODES( myVolumeNodes, 8, 11 );
610 SWAP_NODES( myVolumeNodes, 9, 10 );
611 SWAP_NODES( myVolumeNodes, 12, 15 );
612 SWAP_NODES( myVolumeNodes, 13, 14 );
613 SWAP_NODES( myVolumeNodes, 17, 19 );
616 SWAP_NODES( myVolumeNodes, 1, 3 );
617 SWAP_NODES( myVolumeNodes, 5, 7 );
618 SWAP_NODES( myVolumeNodes, 8, 11 );
619 SWAP_NODES( myVolumeNodes, 9, 10 );
620 SWAP_NODES( myVolumeNodes, 12, 15 );
621 SWAP_NODES( myVolumeNodes, 13, 14 );
622 SWAP_NODES( myVolumeNodes, 17, 19 );
623 SWAP_NODES( myVolumeNodes, 21, 24 );
624 SWAP_NODES( myVolumeNodes, 22, 23 );
630 //=======================================================================
631 //function : GetVolumeType
633 //=======================================================================
635 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
640 switch( myVolumeNodes.size() ) {
641 case 4: return TETRA;
642 case 5: return PYRAM;
643 case 6: return PENTA;
645 case 12: return HEX_PRISM;
646 case 10: return QUAD_TETRA;
647 case 13: return QUAD_PYRAM;
648 case 15: return QUAD_PENTA;
649 case 20: return QUAD_HEXA;
650 case 27: return QUAD_HEXA;
657 //=======================================================================
658 //function : getTetraVolume
660 //=======================================================================
662 static double getTetraVolume(const SMDS_MeshNode* n1,
663 const SMDS_MeshNode* n2,
664 const SMDS_MeshNode* n3,
665 const SMDS_MeshNode* n4)
667 double p1[3], p2[3], p3[3], p4[3];
673 double Q1 = -(p1[ 0 ]-p2[ 0 ])*(p3[ 1 ]*p4[ 2 ]-p4[ 1 ]*p3[ 2 ]);
674 double Q2 = (p1[ 0 ]-p3[ 0 ])*(p2[ 1 ]*p4[ 2 ]-p4[ 1 ]*p2[ 2 ]);
675 double R1 = -(p1[ 0 ]-p4[ 0 ])*(p2[ 1 ]*p3[ 2 ]-p3[ 1 ]*p2[ 2 ]);
676 double R2 = -(p2[ 0 ]-p3[ 0 ])*(p1[ 1 ]*p4[ 2 ]-p4[ 1 ]*p1[ 2 ]);
677 double S1 = (p2[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p3[ 2 ]-p3[ 1 ]*p1[ 2 ]);
678 double S2 = -(p3[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p2[ 2 ]-p2[ 1 ]*p1[ 2 ]);
680 return (Q1+Q2+R1+R2+S1+S2)/6.0;
683 //=======================================================================
685 //purpose : Return element volume
686 //=======================================================================
688 double SMDS_VolumeTool::GetSize() const
694 if ( myVolume->IsPoly() )
699 // split a polyhedron into tetrahedrons
701 SaveFacet savedFacet( myCurFace );
702 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
703 for ( int f = 0; f < NbFaces(); ++f )
706 XYZ area (0,0,0), p1( myCurFace.myNodes[0] );
707 for ( int n = 0; n < myCurFace.myNbNodes; ++n )
709 XYZ p2( myCurFace.myNodes[ n+1 ]);
710 area = area + p1.Crossed( p2 );
719 const static int ind[] = {
720 0, 1, 3, 6, 11, 23, 31, 44, 58, 78 };
721 const static int vtab[][4] = { // decomposition into tetra in the order of enum VolumeType
754 // quadratic tetrahedron
779 // quadratic pentahedron
796 // quadratic hexahedron
821 int type = GetVolumeType();
823 int n2 = ind[type+1];
825 for (int i = n1; i < n2; i++) {
826 V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
827 myVolumeNodes[ vtab[i][1] ],
828 myVolumeNodes[ vtab[i][2] ],
829 myVolumeNodes[ vtab[i][3] ]);
835 //=======================================================================
836 //function : GetBaryCenter
838 //=======================================================================
840 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
846 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
847 X += myVolumeNodes[ i ]->X();
848 Y += myVolumeNodes[ i ]->Y();
849 Z += myVolumeNodes[ i ]->Z();
851 X /= myVolumeNodes.size();
852 Y /= myVolumeNodes.size();
853 Z /= myVolumeNodes.size();
858 //================================================================================
860 * \brief Classify a point
861 * \param tol - thickness of faces
863 //================================================================================
865 bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
867 // LIMITATION: for convex volumes only
869 for ( int iF = 0; iF < myNbFaces; ++iF )
872 if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z ))
874 if ( !IsFaceExternal( iF ))
875 faceNormal = XYZ() - faceNormal; // reverse
877 XYZ face2p( p - XYZ( myCurFace.myNodes[0] ));
878 if ( face2p.Dot( faceNormal ) > tol )
884 //=======================================================================
885 //function : SetExternalNormal
886 //purpose : Node order will be so that faces normals are external
887 //=======================================================================
889 void SMDS_VolumeTool::SetExternalNormal ()
891 myExternalFaces = true;
892 myCurFace.myIndex = -1;
895 //=======================================================================
896 //function : NbFaceNodes
897 //purpose : Return number of nodes in the array of face nodes
898 //=======================================================================
900 int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const
902 if ( !setFace( faceIndex ))
904 return myCurFace.myNbNodes;
907 //=======================================================================
908 //function : GetFaceNodes
909 //purpose : Return pointer to the array of face nodes.
910 // To comfort link iteration, the array
911 // length == NbFaceNodes( faceIndex ) + 1 and
912 // the last node == the first one.
913 //=======================================================================
915 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const
917 if ( !setFace( faceIndex ))
919 return &myCurFace.myNodes[0];
922 //=======================================================================
923 //function : GetFaceNodesIndices
924 //purpose : Return pointer to the array of face nodes indices
925 // To comfort link iteration, the array
926 // length == NbFaceNodes( faceIndex ) + 1 and
927 // the last node index == the first one.
928 //=======================================================================
930 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const
932 if ( !setFace( faceIndex ))
935 return myCurFace.myNodeIndices;
938 //=======================================================================
939 //function : GetFaceNodes
940 //purpose : Return a set of face nodes.
941 //=======================================================================
943 bool SMDS_VolumeTool::GetFaceNodes (int faceIndex,
944 set<const SMDS_MeshNode*>& theFaceNodes ) const
946 if ( !setFace( faceIndex ))
949 theFaceNodes.clear();
950 theFaceNodes.insert( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
957 struct NLink : public std::pair<int,int>
960 NLink(const SMDS_MeshNode* n1=0, const SMDS_MeshNode* n2=0, int ori=1 )
964 if (( myOri = ( n1->GetID() < n2->GetID() )))
967 second = n2->GetID();
973 second = n1->GetID();
979 myOri = first = second = 0;
982 //int Node1() const { return myOri == -1 ? second : first; }
984 //bool IsSameOri( const std::pair<int,int>& link ) const { return link.first == Node1(); }
988 //=======================================================================
989 //function : IsFaceExternal
990 //purpose : Check normal orientation of a given face
991 //=======================================================================
993 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const
995 if ( myExternalFaces || !myVolume )
998 if ( !myPolyedre ) // all classical volumes have external facet normals
1001 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1003 if ( myPolyFacetOri[ faceIndex ])
1004 return myPolyFacetOri[ faceIndex ] > 0;
1006 int ori = 0; // -1-in, +1-out, 0-undef
1007 double minProj, maxProj;
1008 if ( projectNodesToNormal( faceIndex, minProj, maxProj ))
1010 // all nodes are on the same side of the facet
1011 ori = ( minProj < 0 ? +1 : -1 );
1012 me->myPolyFacetOri[ faceIndex ] = ori;
1014 if ( !myFwdLinks.empty() ) // concave polyhedron; collect oriented links
1015 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1017 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1018 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1023 SaveFacet savedFacet( myCurFace );
1025 // concave polyhedron
1027 if ( myFwdLinks.empty() ) // get links of the least ambiguously oriented facet
1029 for ( size_t i = 0; i < myPolyFacetOri.size() && !ori; ++i )
1030 ori = myPolyFacetOri[ i ];
1032 if ( !ori ) // none facet is oriented yet
1034 // find the least ambiguously oriented facet
1035 int faceMostConvex = -1;
1036 std::map< double, int > convexity2face;
1037 for ( size_t iF = 0; iF < myPolyFacetOri.size() && faceMostConvex < 0; ++iF )
1039 if ( projectNodesToNormal( iF, minProj, maxProj ))
1041 // all nodes are on the same side of the facet
1042 me->myPolyFacetOri[ iF ] = ( minProj < 0 ? +1 : -1 );
1043 faceMostConvex = iF;
1047 ori = ( -minProj < maxProj ? -1 : +1 );
1048 double convexity = std::min( -minProj, maxProj ) / std::max( -minProj, maxProj );
1049 convexity2face.insert( make_pair( convexity, iF * ori ));
1052 if ( faceMostConvex < 0 ) // none facet has nodes on the same side
1054 // use the least ambiguous facet
1055 faceMostConvex = convexity2face.begin()->second;
1056 ori = ( faceMostConvex < 0 ? -1 : +1 );
1057 faceMostConvex = std::abs( faceMostConvex );
1058 me->myPolyFacetOri[ faceMostConvex ] = ori;
1061 // collect links of the oriented facets in myFwdLinks
1062 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1064 ori = myPolyFacetOri[ iF ];
1065 if ( !ori ) continue;
1067 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1069 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1070 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1075 // compare orientation of links of the facet with myFwdLinks
1077 setFace( faceIndex );
1078 vector< NLink > links( myCurFace.myNbNodes ), links2;
1079 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1081 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1082 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1083 if ( l2o != myFwdLinks.end() )
1084 ori = link.myOri * l2o->second * -1;
1087 while ( !ori ) // the facet has no common links with already oriented facets
1089 // orient and collect links of other non-oriented facets
1090 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1092 if ( myPolyFacetOri[ iF ] ) continue; // already oriented
1096 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1098 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1099 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1100 if ( l2o != myFwdLinks.end() )
1101 ori = link.myOri * l2o->second * -1;
1102 links2.push_back( link );
1104 if ( ori ) // one more facet oriented
1106 me->myPolyFacetOri[ iF ] = ori;
1107 for ( size_t i = 0; i < links2.size(); ++i )
1108 me->myFwdLinks.insert( make_pair( links2[i], links2[i].myOri * ori ));
1113 return false; // error in algorithm: infinite loop
1115 // try to orient the facet again
1117 for ( size_t i = 0; i < links.size() && !ori; ++i )
1119 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( links[i] );
1120 if ( l2o != myFwdLinks.end() )
1121 ori = links[i].myOri * l2o->second * -1;
1123 me->myPolyFacetOri[ faceIndex ] = ori;
1129 //=======================================================================
1130 //function : projectNodesToNormal
1131 //purpose : compute min and max projections of all nodes to normal of a facet.
1132 //=======================================================================
1134 bool SMDS_VolumeTool::projectNodesToNormal( int faceIndex,
1136 double& maxProj ) const
1138 minProj = std::numeric_limits<double>::max();
1139 maxProj = std::numeric_limits<double>::min();
1142 if ( !GetFaceNormal( faceIndex, normal.x, normal.y, normal.z ))
1144 XYZ p0 ( myCurFace.myNodes[0] );
1145 for ( size_t i = 0; i < myVolumeNodes.size(); ++i )
1147 if ( std::find( myCurFace.myNodes.begin() + 1,
1148 myCurFace.myNodes.end(),
1149 myVolumeNodes[ i ] ) != myCurFace.myNodes.end() )
1150 continue; // node of the faceIndex-th facet
1152 double proj = normal.Dot( XYZ( myVolumeNodes[ i ]) - p0 );
1153 if ( proj < minProj ) minProj = proj;
1154 if ( proj > maxProj ) maxProj = proj;
1156 const double tol = 1e-7;
1159 bool diffSize = ( minProj * maxProj < 0 );
1162 // minProj = -minProj;
1164 // else if ( minProj < 0 )
1166 // minProj = -minProj;
1167 // maxProj = -maxProj;
1170 return !diffSize; // ? 0 : (minProj >= 0);
1173 //=======================================================================
1174 //function : GetFaceNormal
1175 //purpose : Return a normal to a face
1176 //=======================================================================
1178 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) const
1180 if ( !setFace( faceIndex ))
1183 const int iQuad = ( !myPolyedre && myCurFace.myNbNodes > 6 ) ? 2 : 1;
1184 XYZ p1 ( myCurFace.myNodes[0*iQuad] );
1185 XYZ p2 ( myCurFace.myNodes[1*iQuad] );
1186 XYZ p3 ( myCurFace.myNodes[2*iQuad] );
1187 XYZ aVec12( p2 - p1 );
1188 XYZ aVec13( p3 - p1 );
1189 XYZ cross = aVec12.Crossed( aVec13 );
1191 if ( myCurFace.myNbNodes >3*iQuad ) {
1192 XYZ p4 ( myCurFace.myNodes[3*iQuad] );
1193 XYZ aVec14( p4 - p1 );
1194 XYZ cross2 = aVec13.Crossed( aVec14 );
1195 cross = cross + cross2;
1198 double size = cross.Magnitude();
1199 if ( size <= numeric_limits<double>::min() )
1209 //================================================================================
1211 * \brief Return barycenter of a face
1213 //================================================================================
1215 bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const
1217 if ( !setFace( faceIndex ))
1221 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1223 X += myCurFace.myNodes[i]->X() / myCurFace.myNbNodes;
1224 Y += myCurFace.myNodes[i]->Y() / myCurFace.myNbNodes;
1225 Z += myCurFace.myNodes[i]->Z() / myCurFace.myNbNodes;
1230 //=======================================================================
1231 //function : GetFaceArea
1232 //purpose : Return face area
1233 //=======================================================================
1235 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
1238 if ( !setFace( faceIndex ))
1241 XYZ p1 ( myCurFace.myNodes[0] );
1242 XYZ p2 ( myCurFace.myNodes[1] );
1243 XYZ p3 ( myCurFace.myNodes[2] );
1244 XYZ aVec12( p2 - p1 );
1245 XYZ aVec13( p3 - p1 );
1246 area += aVec12.Crossed( aVec13 ).Magnitude();
1248 if (myVolume->IsPoly())
1250 for ( int i = 3; i < myCurFace.myNbNodes; ++i )
1252 XYZ pI ( myCurFace.myNodes[i] );
1253 XYZ aVecI( pI - p1 );
1254 area += aVec13.Crossed( aVecI ).Magnitude();
1260 if ( myCurFace.myNbNodes == 4 ) {
1261 XYZ p4 ( myCurFace.myNodes[3] );
1262 XYZ aVec14( p4 - p1 );
1263 area += aVec14.Crossed( aVec13 ).Magnitude();
1269 //================================================================================
1271 * \brief Return index of the node located at face center of a quadratic element like HEX27
1273 //================================================================================
1275 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
1277 if ( myAllFacesNbNodes && myVolumeNodes.size() == 27 ) // classic element with 27 nodes
1279 switch ( faceIndex ) {
1283 return faceIndex + 19;
1289 //=======================================================================
1290 //function : GetOppFaceIndex
1291 //purpose : Return index of the opposite face if it exists, else -1.
1292 //=======================================================================
1294 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1298 MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1302 const int nbHoriFaces = 2;
1304 if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1305 switch ( myVolumeNodes.size() ) {
1308 if ( faceIndex == 0 || faceIndex == 1 )
1309 ind = 1 - faceIndex;
1313 if ( faceIndex <= 1 ) // top or bottom
1314 ind = 1 - faceIndex;
1316 const int nbSideFaces = myAllFacesNbNodes[0];
1317 ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1322 ind = GetOppFaceIndexOfHex( faceIndex );
1330 //=======================================================================
1331 //function : GetOppFaceIndexOfHex
1332 //purpose : Return index of the opposite face of the hexahedron
1333 //=======================================================================
1335 int SMDS_VolumeTool::GetOppFaceIndexOfHex( int faceIndex )
1337 return Hexa_oppF[ faceIndex ];
1340 //=======================================================================
1341 //function : IsLinked
1342 //purpose : return true if theNode1 is linked with theNode2
1343 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1344 //=======================================================================
1346 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1347 const SMDS_MeshNode* theNode2,
1348 const bool theIgnoreMediumNodes) const
1353 if (myVolume->IsPoly()) {
1355 MESSAGE("Warning: bad volumic element");
1358 if ( !myAllFacesNbNodes ) {
1359 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1360 me->myPolyQuantities = myPolyedre->GetQuantities();
1361 myAllFacesNbNodes = &myPolyQuantities[0];
1363 int from, to = 0, d1 = 1, d2 = 2;
1364 if ( myPolyedre->IsQuadratic() ) {
1365 if ( theIgnoreMediumNodes ) {
1371 vector<const SMDS_MeshNode*>::const_iterator i;
1372 for (int iface = 0; iface < myNbFaces; iface++)
1375 to += myPolyQuantities[iface];
1376 i = std::find( myVolumeNodes.begin() + from, myVolumeNodes.begin() + to, theNode1 );
1377 if ( i != myVolumeNodes.end() )
1379 if (( theNode2 == *( i-d1 ) ||
1380 theNode2 == *( i+d1 )))
1383 (( theNode2 == *( i-d2 ) ||
1384 theNode2 == *( i+d2 ))))
1391 // find nodes indices
1392 int i1 = -1, i2 = -1, nbFound = 0;
1393 for ( size_t i = 0; i < myVolumeNodes.size() && nbFound < 2; i++ )
1395 if ( myVolumeNodes[ i ] == theNode1 )
1397 else if ( myVolumeNodes[ i ] == theNode2 )
1400 return IsLinked( i1, i2 );
1403 //=======================================================================
1404 //function : IsLinked
1405 //purpose : return true if the node with theNode1Index is linked
1406 // with the node with theNode2Index
1407 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1408 //=======================================================================
1410 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1411 const int theNode2Index,
1412 bool theIgnoreMediumNodes) const
1414 if ( myVolume->IsPoly() ) {
1415 return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1418 int minInd = min( theNode1Index, theNode2Index );
1419 int maxInd = max( theNode1Index, theNode2Index );
1421 if ( minInd < 0 || maxInd > (int)myVolumeNodes.size() - 1 || maxInd == minInd )
1424 VolumeType type = GetVolumeType();
1425 if ( myVolume->IsQuadratic() )
1427 int firstMediumInd = myVolume->NbCornerNodes();
1428 if ( minInd >= firstMediumInd )
1429 return false; // both nodes are medium - not linked
1430 if ( maxInd < firstMediumInd ) // both nodes are corners
1432 if ( theIgnoreMediumNodes )
1433 type = quadToLinear(type); // to check linkage of corner nodes only
1435 return false; // corner nodes are not linked directly in a quadratic cell
1443 switch ( maxInd - minInd ) {
1444 case 1: return minInd != 3;
1445 case 3: return minInd == 0 || minInd == 4;
1446 case 4: return true;
1453 switch ( maxInd - minInd ) {
1455 case 3: return true;
1460 switch ( maxInd - minInd ) {
1461 case 1: return minInd != 2;
1462 case 2: return minInd == 0 || minInd == 3;
1463 case 3: return true;
1470 case 0: if( maxInd==4 || maxInd==6 || maxInd==7 ) return true;
1471 case 1: if( maxInd==4 || maxInd==5 || maxInd==8 ) return true;
1472 case 2: if( maxInd==5 || maxInd==6 || maxInd==9 ) return true;
1473 case 3: if( maxInd==7 || maxInd==8 || maxInd==9 ) return true;
1481 case 0: if( maxInd==8 || maxInd==11 || maxInd==16 ) return true;
1482 case 1: if( maxInd==8 || maxInd==9 || maxInd==17 ) return true;
1483 case 2: if( maxInd==9 || maxInd==10 || maxInd==18 ) return true;
1484 case 3: if( maxInd==10 || maxInd==11 || maxInd==19 ) return true;
1485 case 4: if( maxInd==12 || maxInd==15 || maxInd==16 ) return true;
1486 case 5: if( maxInd==12 || maxInd==13 || maxInd==17 ) return true;
1487 case 6: if( maxInd==13 || maxInd==14 || maxInd==18 ) return true;
1488 case 7: if( maxInd==14 || maxInd==15 || maxInd==19 ) return true;
1496 case 0: if( maxInd==5 || maxInd==8 || maxInd==9 ) return true;
1497 case 1: if( maxInd==5 || maxInd==6 || maxInd==10 ) return true;
1498 case 2: if( maxInd==6 || maxInd==7 || maxInd==11 ) return true;
1499 case 3: if( maxInd==7 || maxInd==8 || maxInd==12 ) return true;
1500 case 4: if( maxInd==9 || maxInd==10 || maxInd==11 || maxInd==12 ) return true;
1508 case 0: if( maxInd==6 || maxInd==8 || maxInd==12 ) return true;
1509 case 1: if( maxInd==6 || maxInd==7 || maxInd==13 ) return true;
1510 case 2: if( maxInd==7 || maxInd==8 || maxInd==14 ) return true;
1511 case 3: if( maxInd==9 || maxInd==11 || maxInd==12 ) return true;
1512 case 4: if( maxInd==9 || maxInd==10 || maxInd==13 ) return true;
1513 case 5: if( maxInd==10 || maxInd==11 || maxInd==14 ) return true;
1520 const int diff = maxInd-minInd;
1521 if ( diff > 6 ) return false;// not linked top and bottom
1522 if ( diff == 6 ) return true; // linked top and bottom
1523 return diff == 1 || diff == 7;
1530 //=======================================================================
1531 //function : GetNodeIndex
1532 //purpose : Return an index of theNode
1533 //=======================================================================
1535 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1538 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
1539 if ( myVolumeNodes[ i ] == theNode )
1546 //================================================================================
1548 * \brief Fill vector with boundary faces existing in the mesh
1549 * \param faces - vector of found nodes
1550 * \retval int - nb of found faces
1552 //================================================================================
1554 int SMDS_VolumeTool::GetAllExistingFaces(vector<const SMDS_MeshElement*> & faces) const
1557 SaveFacet savedFacet( myCurFace );
1559 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1561 if ( const SMDS_MeshElement* face = SMDS_Mesh::FindFace( myCurFace.myNodes ))
1562 faces.push_back( face );
1565 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1566 const SMDS_MeshFace* face = 0;
1567 const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1568 switch ( NbFaceNodes( iF )) {
1570 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1572 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1574 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1575 nodes[3], nodes[4], nodes[5]); break;
1577 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1578 nodes[4], nodes[5], nodes[6], nodes[7]); break;
1581 faces.push_back( face );
1583 return faces.size();
1587 //================================================================================
1589 * \brief Fill vector with boundary edges existing in the mesh
1590 * \param edges - vector of found edges
1591 * \retval int - nb of found faces
1593 //================================================================================
1595 int SMDS_VolumeTool::GetAllExistingEdges(vector<const SMDS_MeshElement*> & edges) const
1598 edges.reserve( myVolumeNodes.size() * 2 );
1599 for ( size_t i = 0; i < myVolumeNodes.size()-1; ++i ) {
1600 for ( size_t j = i + 1; j < myVolumeNodes.size(); ++j ) {
1601 if ( IsLinked( i, j )) {
1602 const SMDS_MeshElement* edge =
1603 SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1605 edges.push_back( edge );
1609 return edges.size();
1612 //================================================================================
1614 * \brief Return minimal square distance between connected corner nodes
1616 //================================================================================
1618 double SMDS_VolumeTool::MinLinearSize2() const
1620 double minSize = 1e+100;
1621 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1623 SaveFacet savedFacet( myCurFace );
1625 // it seems that compute distance twice is faster than organization of a sole computing
1626 myCurFace.myIndex = -1;
1627 for ( int iF = 0; iF < myNbFaces; ++iF )
1630 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1632 XYZ n1( myCurFace.myNodes[ iN ]);
1633 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1634 minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
1641 //================================================================================
1643 * \brief Return maximal square distance between connected corner nodes
1645 //================================================================================
1647 double SMDS_VolumeTool::MaxLinearSize2() const
1649 double maxSize = -1e+100;
1650 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1652 SaveFacet savedFacet( myCurFace );
1654 // it seems that compute distance twice is faster than organization of a sole computing
1655 myCurFace.myIndex = -1;
1656 for ( int iF = 0; iF < myNbFaces; ++iF )
1659 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1661 XYZ n1( myCurFace.myNodes[ iN ]);
1662 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1663 maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude());
1670 //================================================================================
1672 * \brief fast check that only one volume is build on the face nodes
1673 * This check is valid for conformal meshes only
1675 //================================================================================
1677 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1679 const bool isFree = true;
1681 if ( !setFace( faceIndex ))
1684 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1686 const int di = myVolume->IsQuadratic() ? 2 : 1;
1687 const int nbN = ( myCurFace.myNbNodes/di <= 4 && !IsPoly()) ? 3 : myCurFace.myNbNodes/di; // nb nodes to check
1689 SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
1690 while ( eIt->more() )
1692 const SMDS_MeshElement* vol = eIt->next();
1693 if ( vol == myVolume )
1696 for ( iN = 1; iN < nbN; ++iN )
1697 if ( vol->GetNodeIndex( nodes[ iN*di ]) < 0 )
1699 if ( iN == nbN ) // nbN nodes are shared with vol
1701 // if ( vol->IsPoly() || vol->NbFaces() > 6 ) // vol is polyhed or hex prism
1703 // int nb = myCurFace.myNbNodes;
1704 // if ( myVolume->GetEntityType() != vol->GetEntityType() )
1705 // nb -= ( GetCenterNodeIndex(0) > 0 );
1706 // set<const SMDS_MeshNode*> faceNodes( nodes, nodes + nb );
1707 // if ( SMDS_VolumeTool( vol ).GetFaceIndex( faceNodes ) < 0 )
1710 if ( otherVol ) *otherVol = vol;
1714 if ( otherVol ) *otherVol = 0;
1718 //================================================================================
1720 * \brief Thorough check that only one volume is build on the face nodes
1722 //================================================================================
1724 bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1726 const bool isFree = true;
1728 if (!setFace( faceIndex ))
1731 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1732 const int nbFaceNodes = myCurFace.myNbNodes;
1734 // evaluate nb of face nodes shared by other volumes
1735 int maxNbShared = -1;
1736 typedef map< const SMDS_MeshElement*, int > TElemIntMap;
1737 TElemIntMap volNbShared;
1738 TElemIntMap::iterator vNbIt;
1739 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1740 const SMDS_MeshNode* n = nodes[ iNode ];
1741 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1742 while ( eIt->more() ) {
1743 const SMDS_MeshElement* elem = eIt->next();
1744 if ( elem != myVolume ) {
1745 vNbIt = volNbShared.insert( make_pair( elem, 0 )).first;
1747 if ( vNbIt->second > maxNbShared )
1748 maxNbShared = vNbIt->second;
1752 if ( maxNbShared < 3 )
1753 return isFree; // is free
1755 // find volumes laying on the opposite side of the face
1756 // and sharing all nodes
1757 XYZ intNormal; // internal normal
1758 GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
1759 if ( IsFaceExternal( faceIndex ))
1760 intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
1761 XYZ p0 ( nodes[0] ), baryCenter;
1762 for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); ) {
1763 const int& nbShared = (*vNbIt).second;
1764 if ( nbShared >= 3 ) {
1765 SMDS_VolumeTool volume( (*vNbIt).first );
1766 volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
1767 XYZ intNormal2( baryCenter - p0 );
1768 if ( intNormal.Dot( intNormal2 ) < 0 ) {
1770 if ( nbShared >= nbFaceNodes )
1772 // a volume shares the whole facet
1773 if ( otherVol ) *otherVol = vNbIt->first;
1780 // remove a volume from volNbShared map
1781 volNbShared.erase( vNbIt++ );
1784 // here volNbShared contains only volumes laying on the opposite side of
1785 // the face and sharing 3 or more but not all face nodes with myVolume
1786 if ( volNbShared.size() < 2 ) {
1787 return isFree; // is free
1790 // check if the whole area of a face is shared
1791 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1793 const SMDS_MeshNode* n = nodes[ iNode ];
1794 // check if n is shared by one of volumes of volNbShared
1795 bool isShared = false;
1796 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1797 while ( eIt->more() && !isShared )
1798 isShared = volNbShared.count( eIt->next() );
1802 if ( otherVol ) *otherVol = volNbShared.begin()->first;
1805 // if ( !myVolume->IsPoly() )
1807 // bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
1808 // for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1809 // SMDS_VolumeTool volume( (*vNbIt).first );
1810 // bool prevLinkShared = false;
1811 // int nbSharedLinks = 0;
1812 // for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1813 // bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
1814 // if ( linkShared )
1816 // if ( linkShared && prevLinkShared &&
1817 // volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
1818 // isShared[ iNode ] = true;
1819 // prevLinkShared = linkShared;
1821 // if ( nbSharedLinks == nbFaceNodes )
1822 // return !free; // is not free
1823 // if ( nbFaceNodes == 4 ) {
1824 // // check traingle parts 1 & 3
1825 // if ( isShared[1] && isShared[3] )
1826 // return !free; // is not free
1827 // // check triangle parts 0 & 2;
1828 // // 0 part could not be checked in the loop; check it here
1829 // if ( isShared[2] && prevLinkShared &&
1830 // volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
1831 // volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
1832 // return !free; // is not free
1839 //=======================================================================
1840 //function : GetFaceIndex
1841 //purpose : Return index of a face formed by theFaceNodes
1842 //=======================================================================
1844 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes,
1845 const int theFaceIndexHint ) const
1847 if ( theFaceIndexHint >= 0 )
1849 int nbNodes = NbFaceNodes( theFaceIndexHint );
1850 if ( nbNodes == (int) theFaceNodes.size() )
1852 const SMDS_MeshNode** nodes = GetFaceNodes( theFaceIndexHint );
1854 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1859 return theFaceIndexHint;
1862 for ( int iFace = 0; iFace < myNbFaces; iFace++ )
1864 if ( iFace == theFaceIndexHint )
1866 int nbNodes = NbFaceNodes( iFace );
1867 if ( nbNodes == (int) theFaceNodes.size() )
1869 const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
1871 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1882 //=======================================================================
1883 //function : GetFaceIndex
1884 //purpose : Return index of a face formed by theFaceNodes
1885 //=======================================================================
1887 /*int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
1889 for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1890 const int* nodes = GetFaceNodesIndices( iFace );
1891 int nbFaceNodes = NbFaceNodes( iFace );
1893 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1894 nodeSet.insert( nodes[ iNode ] );
1895 if ( theFaceNodesIndices == nodeSet )
1901 //=======================================================================
1902 //function : setFace
1904 //=======================================================================
1906 bool SMDS_VolumeTool::setFace( int faceIndex ) const
1911 if ( myCurFace.myIndex == faceIndex )
1914 myCurFace.myIndex = -1;
1916 if ( faceIndex < 0 || faceIndex >= NbFaces() )
1919 if (myVolume->IsPoly())
1922 MESSAGE("Warning: bad volumic element");
1927 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1928 if ( !myAllFacesNbNodes ) {
1929 me->myPolyQuantities = myPolyedre->GetQuantities();
1930 myAllFacesNbNodes = &myPolyQuantities[0];
1932 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
1933 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
1934 me->myPolyIndices.resize( myCurFace.myNbNodes + 1 );
1935 myCurFace.myNodeIndices = & me->myPolyIndices[0];
1936 int shift = std::accumulate( myAllFacesNbNodes, myAllFacesNbNodes+faceIndex, 0 );
1937 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
1939 myCurFace.myNodes [ iNode ] = myVolumeNodes[ shift + iNode ];
1940 myCurFace.myNodeIndices[ iNode ] = shift + iNode;
1942 myCurFace.myNodes [ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ]; // last = first
1943 myCurFace.myNodeIndices[ myCurFace.myNbNodes ] = myCurFace.myNodeIndices[ 0 ];
1945 // check orientation
1946 if (myExternalFaces)
1948 myCurFace.myIndex = faceIndex; // avoid infinite recursion in IsFaceExternal()
1949 myExternalFaces = false; // force normal computation by IsFaceExternal()
1950 if ( !IsFaceExternal( faceIndex ))
1951 std::reverse( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
1952 myExternalFaces = true;
1957 if ( !myAllFacesNodeIndices_F )
1959 // choose data for an element type
1960 switch ( myVolumeNodes.size() ) {
1962 myAllFacesNodeIndices_F = &Tetra_F [0][0];
1963 //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
1964 myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
1965 myAllFacesNbNodes = Tetra_nbN;
1966 myMaxFaceNbNodes = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
1969 myAllFacesNodeIndices_F = &Pyramid_F [0][0];
1970 //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
1971 myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
1972 myAllFacesNbNodes = Pyramid_nbN;
1973 myMaxFaceNbNodes = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
1976 myAllFacesNodeIndices_F = &Penta_F [0][0];
1977 //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
1978 myAllFacesNodeIndices_RE = &Penta_RE[0][0];
1979 myAllFacesNbNodes = Penta_nbN;
1980 myMaxFaceNbNodes = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
1983 myAllFacesNodeIndices_F = &Hexa_F [0][0];
1984 ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
1985 myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
1986 myAllFacesNbNodes = Hexa_nbN;
1987 myMaxFaceNbNodes = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
1990 myAllFacesNodeIndices_F = &QuadTetra_F [0][0];
1991 //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
1992 myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
1993 myAllFacesNbNodes = QuadTetra_nbN;
1994 myMaxFaceNbNodes = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
1997 myAllFacesNodeIndices_F = &QuadPyram_F [0][0];
1998 //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
1999 myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
2000 myAllFacesNbNodes = QuadPyram_nbN;
2001 myMaxFaceNbNodes = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
2004 myAllFacesNodeIndices_F = &QuadPenta_F [0][0];
2005 //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
2006 myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
2007 myAllFacesNbNodes = QuadPenta_nbN;
2008 myMaxFaceNbNodes = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
2012 myAllFacesNodeIndices_F = &QuadHexa_F [0][0];
2013 //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
2014 myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
2015 myAllFacesNbNodes = QuadHexa_nbN;
2016 myMaxFaceNbNodes = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
2017 if ( !myIgnoreCentralNodes && myVolumeNodes.size() == 27 )
2019 myAllFacesNodeIndices_F = &TriQuadHexa_F [0][0];
2020 //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
2021 myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
2022 myAllFacesNbNodes = TriQuadHexa_nbN;
2023 myMaxFaceNbNodes = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
2027 myAllFacesNodeIndices_F = &HexPrism_F [0][0];
2028 //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
2029 myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
2030 myAllFacesNbNodes = HexPrism_nbN;
2031 myMaxFaceNbNodes = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
2037 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2038 // if ( myExternalFaces )
2039 // myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2041 // myCurFace.myNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
2042 myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2045 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2046 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2047 myCurFace.myNodes[ iNode ] = myVolumeNodes[ myCurFace.myNodeIndices[ iNode ]];
2048 myCurFace.myNodes[ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ];
2051 myCurFace.myIndex = faceIndex;
2056 //=======================================================================
2057 //function : GetType
2058 //purpose : return VolumeType by nb of nodes in a volume
2059 //=======================================================================
2061 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
2063 switch ( nbNodes ) {
2064 case 4: return TETRA;
2065 case 5: return PYRAM;
2066 case 6: return PENTA;
2067 case 8: return HEXA;
2068 case 10: return QUAD_TETRA;
2069 case 13: return QUAD_PYRAM;
2070 case 15: return QUAD_PENTA;
2072 case 27: return QUAD_HEXA;
2073 case 12: return HEX_PRISM;
2074 default:return UNKNOWN;
2078 //=======================================================================
2079 //function : NbFaces
2080 //purpose : return nb of faces by volume type
2081 //=======================================================================
2083 int SMDS_VolumeTool::NbFaces( VolumeType type )
2087 case QUAD_TETRA: return 4;
2089 case QUAD_PYRAM: return 5;
2091 case QUAD_PENTA: return 5;
2093 case QUAD_HEXA : return 6;
2094 case HEX_PRISM : return 8;
2099 //================================================================================
2101 * \brief Useful to know nb of corner nodes of a quadratic volume
2102 * \param type - volume type
2103 * \retval int - nb of corner nodes
2105 //================================================================================
2107 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
2111 case QUAD_TETRA: return 4;
2113 case QUAD_PYRAM: return 5;
2115 case QUAD_PENTA: return 6;
2117 case QUAD_HEXA : return 8;
2118 case HEX_PRISM : return 12;
2125 //=======================================================================
2126 //function : GetFaceNodesIndices
2127 //purpose : Return the array of face nodes indices
2128 // To comfort link iteration, the array
2129 // length == NbFaceNodes( faceIndex ) + 1 and
2130 // the last node index == the first one.
2131 //=======================================================================
2133 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
2138 case TETRA: return Tetra_F[ faceIndex ];
2139 case PYRAM: return Pyramid_F[ faceIndex ];
2140 case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
2141 case HEXA: return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
2142 case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
2143 case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
2144 case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
2145 // what about SMDSEntity_TriQuad_Hexa?
2146 case QUAD_HEXA: return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
2147 case HEX_PRISM: return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
2153 //=======================================================================
2154 //function : NbFaceNodes
2155 //purpose : Return number of nodes in the array of face nodes
2156 //=======================================================================
2158 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
2162 case TETRA: return Tetra_nbN[ faceIndex ];
2163 case PYRAM: return Pyramid_nbN[ faceIndex ];
2164 case PENTA: return Penta_nbN[ faceIndex ];
2165 case HEXA: return Hexa_nbN[ faceIndex ];
2166 case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
2167 case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
2168 case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
2169 // what about SMDSEntity_TriQuad_Hexa?
2170 case QUAD_HEXA: return QuadHexa_nbN[ faceIndex ];
2171 case HEX_PRISM: return HexPrism_nbN[ faceIndex ];
2177 //=======================================================================
2178 //function : Element
2179 //purpose : return element
2180 //=======================================================================
2182 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
2184 return static_cast<const SMDS_MeshVolume*>( myVolume );
2187 //=======================================================================
2189 //purpose : return element ID
2190 //=======================================================================
2192 int SMDS_VolumeTool::ID() const
2194 return myVolume ? myVolume->GetID() : 0;