1 // Copyright (C) 2007-2020 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>
48 // ======================================================
49 // Node indices in faces depending on volume orientation
50 // making most faces normals external
51 // ======================================================
52 // For all elements, 0-th face is bottom based on the first nodes.
53 // For prismatic elements (tetra,hexa,prisms), 1-th face is a top one.
54 // For all elements, side faces follow order of bottom nodes
55 // ======================================================
63 // N0 +---|---+ N1 TETRAHEDRON
71 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
72 { 0, 1, 2, 0 }, // All faces have external normals
76 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
77 { 0, 2, 1, 0 }, // All faces have external normals
81 static int Tetra_nbN [] = { 3, 3, 3, 3 };
86 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
87 { 0, 1, 2, 3, 0 }, // All faces have external normals
93 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
94 { 0, 3, 2, 1, 0 }, // All faces but a bottom have external normals
99 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
110 // | / \ | PENTAHEDRON
116 static int Penta_F [5][5] = { // FORWARD
117 { 0, 1, 2, 0, 0 }, // All faces have external normals
118 { 3, 5, 4, 3, 3 }, // 0 is bottom, 1 is top face
122 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
128 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
135 // N4+----------+N7 |
136 // | | | | HEXAHEDRON
137 // | N1+------|---+N2
143 static int Hexa_F [6][5] = { // FORWARD
145 { 4, 7, 6, 5, 4 }, // all face normals are external
150 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
152 { 4, 5, 6, 7, 4 }, // all face normals are external
157 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
158 static int Hexa_oppF[] = { 1, 0, 4, 5, 2, 3 }; // oppopsite facet indices
177 static int HexPrism_F [8][7] = { // FORWARD
178 { 0, 1, 2, 3, 4, 5, 0 },
179 { 6,11,10, 9, 8, 7, 6 },
180 { 0, 6, 7, 1, 0, 0, 0 },
181 { 1, 7, 8, 2, 1, 1, 1 },
182 { 2, 8, 9, 3, 2, 2, 2 },
183 { 3, 9,10, 4, 3, 3, 3 },
184 { 4,10,11, 5, 4, 4, 4 },
185 { 5,11, 6, 0, 5, 5, 5 }};
186 static int HexPrism_RE [8][7] = { // REVERSED -> EXTERNAL
187 { 0, 5, 4, 3, 2, 1, 0 },
188 { 6,11,10, 9, 8, 7, 6 },
189 { 0, 6, 7, 1, 0, 0, 0 },
190 { 1, 7, 8, 2, 1, 1, 1 },
191 { 2, 8, 9, 3, 2, 2, 2 },
192 { 3, 9,10, 4, 3, 3, 3 },
193 { 4,10,11, 5, 4, 4, 4 },
194 { 5,11, 6, 0, 5, 5, 5 }};
195 static int HexPrism_nbN [] = { 6, 6, 4, 4, 4, 4, 4, 4 };
204 // N0 +---|---+ N1 TETRAHEDRON
212 static int QuadTetra_F [4][7] = { // FORWARD
213 { 0, 4, 1, 5, 2, 6, 0 }, // All faces have external normals
214 { 0, 7, 3, 8, 1, 4, 0 },
215 { 1, 8, 3, 9, 2, 5, 1 },
216 { 0, 6, 2, 9, 3, 7, 0 }};
217 static int QuadTetra_RE [4][7] = { // REVERSED -> FORWARD (EXTERNAL)
218 { 0, 6, 2, 5, 1, 4, 0 }, // All faces have external normals
219 { 0, 4, 1, 8, 3, 7, 0 },
220 { 1, 5, 2, 9, 3, 8, 1 },
221 { 0, 7, 3, 9, 2, 6, 0 }};
222 static int QuadTetra_nbN [] = { 6, 6, 6, 6 };
232 // | | 9 - middle point for (0,4) etc.
245 static int QuadPyram_F [5][9] = { // FORWARD
246 { 0, 5, 1, 6, 2, 7, 3, 8, 0 }, // All faces have external normals
247 { 0, 9, 4, 10,1, 5, 0, 4, 4 },
248 { 1, 10,4, 11,2, 6, 1, 4, 4 },
249 { 2, 11,4, 12,3, 7, 2, 4, 4 },
250 { 3, 12,4, 9, 0, 8, 3, 4, 4 }};
251 static int QuadPyram_RE [5][9] = { // REVERSED -> FORWARD (EXTERNAL)
252 { 0, 8, 3, 7, 2, 6, 1, 5, 0 }, // All faces but a bottom have external normals
253 { 0, 5, 1, 10,4, 9, 0, 4, 4 },
254 { 1, 6, 2, 11,4, 10,1, 4, 4 },
255 { 2, 7, 3, 12,4, 11,2, 4, 4 },
256 { 3, 8, 0, 9, 4, 12,3, 4, 4 }};
257 static int QuadPyram_nbN [] = { 8, 6, 6, 6, 6 };
281 static int QuadPenta_F [5][9] = { // FORWARD
282 { 0, 6, 1, 7, 2, 8, 0, 0, 0 },
283 { 3, 11,5, 10,4, 9, 3, 3, 3 },
284 { 0, 12,3, 9, 4, 13,1, 6, 0 },
285 { 1, 13,4, 10,5, 14,2, 7, 1 },
286 { 0, 8, 2, 14,5, 11,3, 12,0 }};
287 static int QuadPenta_RE [5][9] = { // REVERSED -> EXTERNAL
288 { 0, 8, 2, 7, 1, 6, 0, 0, 0 },
289 { 3, 9, 4, 10,5, 11,3, 3, 3 },
290 { 0, 6, 1, 13,4, 9, 3, 12,0 },
291 { 1, 7, 2, 14,5, 10,4, 13,1 },
292 { 0, 12,3, 11,5, 14,2, 8, 0 }};
293 static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 };
297 // N5+-----+-----+N6 +-----+-----+
299 // 12+ | 14+ | + | +25 + |
301 // N4+-----+-----+N7 | QUADRATIC +-----+-----+ | Central nodes
302 // | | 15 | | HEXAHEDRON | | | | of tri-quadratic
303 // | | | | | | | | HEXAHEDRON
304 // | 17+ | +18 | + 22+ | +
306 // | | | | | + | 26+ | + |
308 // 16+ | +19 | + | +24 + |
311 // | N1+-----+-|---+N2 | +-----+-|---+
313 // | +8 | +10 | + 20+ | +
315 // N0+-----+-----+N3 +-----+-----+
318 static int QuadHexa_F [6][9] = { // FORWARD
319 { 0, 8, 1, 9, 2, 10,3, 11,0 }, // all face normals are external,
320 { 4, 15,7, 14,6, 13,5, 12,4 },
321 { 0, 16,4, 12,5, 17,1, 8, 0 },
322 { 1, 17,5, 13,6, 18,2, 9, 1 },
323 { 3, 10,2, 18,6, 14,7, 19,3 },
324 { 0, 11,3, 19,7, 15,4, 16,0 }};
325 static int QuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL
326 { 0, 11,3, 10,2, 9, 1, 8, 0 }, // all face normals are external
327 { 4, 12,5, 13,6, 14,7, 15,4 },
328 { 0, 8, 1, 17,5, 12,4, 16,0 },
329 { 1, 9, 2, 18,6, 13,5, 17,1 },
330 { 3, 19,7, 14,6, 18,2, 10,3 },
331 { 0, 16,4, 15,7, 19,3, 11,0 }};
332 static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 };
334 static int TriQuadHexa_F [6][9] = { // FORWARD
335 { 0, 8, 1, 9, 2, 10,3, 11, 20 }, // all face normals are external
336 { 4, 15,7, 14,6, 13,5, 12, 25 },
337 { 0, 16,4, 12,5, 17,1, 8, 21 },
338 { 1, 17,5, 13,6, 18,2, 9, 22 },
339 { 3, 10,2, 18,6, 14,7, 19, 23 },
340 { 0, 11,3, 19,7, 15,4, 16, 24 }};
341 static int TriQuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL
342 { 0, 11,3, 10,2, 9, 1, 8, 20 }, // opposite faces are neighbouring,
343 { 4, 12,5, 13,6, 14,7, 15, 25 }, // all face normals are external
344 { 0, 8, 1, 17,5, 12,4, 16, 21 },
345 { 1, 9, 2, 18,6, 13,5, 17, 22 },
346 { 3, 19,7, 14,6, 18,2, 10, 23 },
347 { 0, 16,4, 15,7, 19,3, 11, 24 }};
348 static int TriQuadHexa_nbN [] = { 9, 9, 9, 9, 9, 9 };
351 // ========================================================
352 // to perform some calculations without linkage to CASCADE
353 // ========================================================
358 XYZ() { x = 0; y = 0; z = 0; }
359 XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
360 XYZ( const XYZ& other ) { x = other.x; y = other.y; z = other.z; }
361 XYZ( const SMDS_MeshNode* n ) { x = n->X(); y = n->Y(); z = n->Z(); }
362 double* data() { return &x; }
363 inline XYZ operator-( const XYZ& other );
364 inline XYZ operator+( const XYZ& other );
365 inline XYZ Crossed( const XYZ& other );
366 inline double Dot( const XYZ& other );
367 inline double Magnitude();
368 inline double SquareMagnitude();
370 inline XYZ XYZ::operator-( const XYZ& Right ) {
371 return XYZ(x - Right.x, y - Right.y, z - Right.z);
373 inline XYZ XYZ::operator+( const XYZ& Right ) {
374 return XYZ(x + Right.x, y + Right.y, z + Right.z);
376 inline XYZ XYZ::Crossed( const XYZ& Right ) {
377 return XYZ (y * Right.z - z * Right.y,
378 z * Right.x - x * Right.z,
379 x * Right.y - y * Right.x);
381 inline double XYZ::Dot( const XYZ& Other ) {
382 return(x * Other.x + y * Other.y + z * Other.z);
384 inline double XYZ::Magnitude() {
385 return sqrt (x * x + y * y + z * z);
387 inline double XYZ::SquareMagnitude() {
388 return (x * x + y * y + z * z);
391 //================================================================================
393 * \brief Return linear type corresponding to a quadratic one
395 //================================================================================
397 SMDS_VolumeTool::VolumeType quadToLinear(SMDS_VolumeTool::VolumeType quadType)
399 SMDS_VolumeTool::VolumeType linType = SMDS_VolumeTool::VolumeType( int(quadType)-4 );
400 const int nbCornersByQuad = SMDS_VolumeTool::NbCornerNodes( quadType );
401 if ( SMDS_VolumeTool::NbCornerNodes( linType ) == nbCornersByQuad )
405 for ( ; iLin < SMDS_VolumeTool::NB_VOLUME_TYPES; ++iLin )
406 if ( SMDS_VolumeTool::NbCornerNodes( SMDS_VolumeTool::VolumeType( iLin )) == nbCornersByQuad)
407 return SMDS_VolumeTool::VolumeType( iLin );
409 return SMDS_VolumeTool::UNKNOWN;
414 //================================================================================
416 * \brief Saver/restorer of a SMDS_VolumeTool::myCurFace
418 //================================================================================
420 struct SMDS_VolumeTool::SaveFacet
422 SMDS_VolumeTool::Facet mySaved;
423 SMDS_VolumeTool::Facet& myToRestore;
424 SaveFacet( SMDS_VolumeTool::Facet& facet ): myToRestore( facet )
427 mySaved.myNodes.swap( facet.myNodes );
431 if ( myToRestore.myIndex != mySaved.myIndex )
432 myToRestore = mySaved;
433 myToRestore.myNodes.swap( mySaved.myNodes );
437 //=======================================================================
438 //function : SMDS_VolumeTool
440 //=======================================================================
442 SMDS_VolumeTool::SMDS_VolumeTool ()
447 //=======================================================================
448 //function : SMDS_VolumeTool
450 //=======================================================================
452 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume,
453 const bool ignoreCentralNodes)
455 Set( theVolume, ignoreCentralNodes );
458 //=======================================================================
459 //function : SMDS_VolumeTool
461 //=======================================================================
463 SMDS_VolumeTool::~SMDS_VolumeTool()
465 myCurFace.myNodeIndices = NULL;
468 //=======================================================================
469 //function : SetVolume
470 //purpose : Set volume to iterate on
471 //=======================================================================
473 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume,
474 const bool ignoreCentralNodes,
475 const std::vector<const SMDS_MeshNode*>* otherNodes)
480 myIgnoreCentralNodes = ignoreCentralNodes;
484 myVolumeNodes.clear();
485 myPolyIndices.clear();
486 myPolyQuantities.clear();
487 myPolyFacetOri.clear();
490 myExternalFaces = false;
492 myAllFacesNodeIndices_F = 0;
493 myAllFacesNodeIndices_RE = 0;
494 myAllFacesNbNodes = 0;
496 myCurFace.myIndex = -1;
497 myCurFace.myNodeIndices = NULL;
498 myCurFace.myNodes.clear();
501 if ( !theVolume || theVolume->GetType() != SMDSAbs_Volume )
504 myVolume = theVolume;
505 myNbFaces = theVolume->NbFaces();
506 if ( myVolume->IsPoly() )
508 myPolyedre = SMDS_Mesh::DownCast<SMDS_MeshVolume>( myVolume );
509 myPolyFacetOri.resize( myNbFaces, 0 );
513 myVolumeNodes.resize( myVolume->NbNodes() );
516 if ( otherNodes->size() != myVolumeNodes.size() )
517 return ( myVolume = 0 );
518 for ( size_t i = 0; i < otherNodes->size(); ++i )
519 if ( ! ( myVolumeNodes[i] = (*otherNodes)[0] ))
520 return ( myVolume = 0 );
524 myVolumeNodes.assign( myVolume->begin_nodes(), myVolume->end_nodes() );
529 return ( myVolume = 0 );
533 // define volume orientation
535 if ( GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z ))
537 const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
538 int topNodeIndex = myVolume->NbCornerNodes() - 1;
539 while ( !IsLinked( 0, topNodeIndex, /*ignoreMediumNodes=*/true )) --topNodeIndex;
540 const SMDS_MeshNode* topNode = myVolumeNodes[ topNodeIndex ];
541 XYZ upDir (topNode->X() - botNode->X(),
542 topNode->Y() - botNode->Y(),
543 topNode->Z() - botNode->Z() );
544 myVolForward = ( botNormal.Dot( upDir ) < 0 );
547 myCurFace.myIndex = -1; // previous setFace(0) didn't take myVolForward into account
552 //=======================================================================
554 //purpose : Inverse volume
555 //=======================================================================
557 #define SWAP_NODES(nodes,i1,i2) \
559 const SMDS_MeshNode* tmp = nodes[ i1 ]; \
560 nodes[ i1 ] = nodes[ i2 ]; \
563 void SMDS_VolumeTool::Inverse ()
565 if ( !myVolume ) return;
567 if (myVolume->IsPoly()) {
568 MESSAGE("Warning: attempt to inverse polyhedral volume");
572 myVolForward = !myVolForward;
573 myCurFace.myIndex = -1;
575 // inverse top and bottom faces
576 switch ( myVolumeNodes.size() ) {
578 SWAP_NODES( myVolumeNodes, 1, 2 );
581 SWAP_NODES( myVolumeNodes, 1, 3 );
584 SWAP_NODES( myVolumeNodes, 1, 2 );
585 SWAP_NODES( myVolumeNodes, 4, 5 );
588 SWAP_NODES( myVolumeNodes, 1, 3 );
589 SWAP_NODES( myVolumeNodes, 5, 7 );
592 SWAP_NODES( myVolumeNodes, 1, 5 );
593 SWAP_NODES( myVolumeNodes, 2, 4 );
594 SWAP_NODES( myVolumeNodes, 7, 11 );
595 SWAP_NODES( myVolumeNodes, 8, 10 );
599 SWAP_NODES( myVolumeNodes, 1, 2 );
600 SWAP_NODES( myVolumeNodes, 4, 6 );
601 SWAP_NODES( myVolumeNodes, 8, 9 );
604 SWAP_NODES( myVolumeNodes, 1, 3 );
605 SWAP_NODES( myVolumeNodes, 5, 8 );
606 SWAP_NODES( myVolumeNodes, 6, 7 );
607 SWAP_NODES( myVolumeNodes, 10, 12 );
610 SWAP_NODES( myVolumeNodes, 1, 2 );
611 SWAP_NODES( myVolumeNodes, 4, 5 );
612 SWAP_NODES( myVolumeNodes, 6, 8 );
613 SWAP_NODES( myVolumeNodes, 9, 11 );
614 SWAP_NODES( myVolumeNodes, 13, 14 );
617 SWAP_NODES( myVolumeNodes, 1, 3 );
618 SWAP_NODES( myVolumeNodes, 5, 7 );
619 SWAP_NODES( myVolumeNodes, 8, 11 );
620 SWAP_NODES( myVolumeNodes, 9, 10 );
621 SWAP_NODES( myVolumeNodes, 12, 15 );
622 SWAP_NODES( myVolumeNodes, 13, 14 );
623 SWAP_NODES( myVolumeNodes, 17, 19 );
626 SWAP_NODES( myVolumeNodes, 1, 3 );
627 SWAP_NODES( myVolumeNodes, 5, 7 );
628 SWAP_NODES( myVolumeNodes, 8, 11 );
629 SWAP_NODES( myVolumeNodes, 9, 10 );
630 SWAP_NODES( myVolumeNodes, 12, 15 );
631 SWAP_NODES( myVolumeNodes, 13, 14 );
632 SWAP_NODES( myVolumeNodes, 17, 19 );
633 SWAP_NODES( myVolumeNodes, 21, 24 );
634 SWAP_NODES( myVolumeNodes, 22, 23 );
640 //=======================================================================
641 //function : GetVolumeType
643 //=======================================================================
645 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
650 switch( myVolumeNodes.size() ) {
651 case 4: return TETRA;
652 case 5: return PYRAM;
653 case 6: return PENTA;
655 case 12: return HEX_PRISM;
656 case 10: return QUAD_TETRA;
657 case 13: return QUAD_PYRAM;
658 case 15: return QUAD_PENTA;
659 case 20: return QUAD_HEXA;
660 case 27: return QUAD_HEXA;
667 //=======================================================================
668 //function : getTetraVolume
670 //=======================================================================
672 static double getTetraVolume(const SMDS_MeshNode* n1,
673 const SMDS_MeshNode* n2,
674 const SMDS_MeshNode* n3,
675 const SMDS_MeshNode* n4)
677 double p1[3], p2[3], p3[3], p4[3];
683 double Q1 = -(p1[ 0 ]-p2[ 0 ])*(p3[ 1 ]*p4[ 2 ]-p4[ 1 ]*p3[ 2 ]);
684 double Q2 = (p1[ 0 ]-p3[ 0 ])*(p2[ 1 ]*p4[ 2 ]-p4[ 1 ]*p2[ 2 ]);
685 double R1 = -(p1[ 0 ]-p4[ 0 ])*(p2[ 1 ]*p3[ 2 ]-p3[ 1 ]*p2[ 2 ]);
686 double R2 = -(p2[ 0 ]-p3[ 0 ])*(p1[ 1 ]*p4[ 2 ]-p4[ 1 ]*p1[ 2 ]);
687 double S1 = (p2[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p3[ 2 ]-p3[ 1 ]*p1[ 2 ]);
688 double S2 = -(p3[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p2[ 2 ]-p2[ 1 ]*p1[ 2 ]);
690 return (Q1+Q2+R1+R2+S1+S2)/6.0;
693 //=======================================================================
695 //purpose : Return element volume
696 //=======================================================================
698 double SMDS_VolumeTool::GetSize() const
704 if ( myVolume->IsPoly() )
709 SaveFacet savedFacet( myCurFace );
711 // split a polyhedron into tetrahedrons
714 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
715 for ( int f = 0; f < NbFaces(); ++f )
718 XYZ area (0,0,0), p1( myCurFace.myNodes[0] );
719 for ( int n = 0; n < myCurFace.myNbNodes; ++n )
721 XYZ p2( myCurFace.myNodes[ n+1 ]);
722 area = area + p1.Crossed( p2 );
726 oriOk = oriOk && IsFaceExternal( f );
729 if ( !oriOk && V > 0 )
734 const static int ind[] = {
735 0, 1, 3, 6, 11, 23, 31, 44, 58, 78 };
736 const static int vtab[][4] = { // decomposition into tetra in the order of enum VolumeType
769 // quadratic tetrahedron
794 // quadratic pentahedron
811 // quadratic hexahedron
836 int type = GetVolumeType();
838 int n2 = ind[type+1];
840 for (int i = n1; i < n2; i++) {
841 V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
842 myVolumeNodes[ vtab[i][1] ],
843 myVolumeNodes[ vtab[i][2] ],
844 myVolumeNodes[ vtab[i][3] ]);
850 //=======================================================================
851 //function : GetBaryCenter
853 //=======================================================================
855 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
861 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
862 X += myVolumeNodes[ i ]->X();
863 Y += myVolumeNodes[ i ]->Y();
864 Z += myVolumeNodes[ i ]->Z();
866 X /= myVolumeNodes.size();
867 Y /= myVolumeNodes.size();
868 Z /= myVolumeNodes.size();
873 //================================================================================
875 * \brief Classify a point
876 * \param tol - thickness of faces
878 //================================================================================
880 bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
882 // LIMITATION: for convex volumes only
884 for ( int iF = 0; iF < myNbFaces; ++iF )
887 if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z ))
889 if ( !IsFaceExternal( iF ))
890 faceNormal = XYZ() - faceNormal; // reverse
892 XYZ face2p( p - XYZ( myCurFace.myNodes[0] ));
893 if ( face2p.Dot( faceNormal ) > tol )
899 //=======================================================================
900 //function : SetExternalNormal
901 //purpose : Node order will be so that faces normals are external
902 //=======================================================================
904 void SMDS_VolumeTool::SetExternalNormal ()
906 myExternalFaces = true;
907 myCurFace.myIndex = -1;
910 //=======================================================================
911 //function : NbFaceNodes
912 //purpose : Return number of nodes in the array of face nodes
913 //=======================================================================
915 int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const
917 if ( !setFace( faceIndex ))
919 return myCurFace.myNbNodes;
922 //=======================================================================
923 //function : GetFaceNodes
924 //purpose : Return pointer to the array of face nodes.
925 // To comfort link iteration, the array
926 // length == NbFaceNodes( faceIndex ) + 1 and
927 // the last node == the first one.
928 //=======================================================================
930 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const
932 if ( !setFace( faceIndex ))
934 return &myCurFace.myNodes[0];
937 //=======================================================================
938 //function : GetFaceNodesIndices
939 //purpose : Return pointer to the array of face nodes indices
940 // To comfort link iteration, the array
941 // length == NbFaceNodes( faceIndex ) + 1 and
942 // the last node index == the first one.
943 //=======================================================================
945 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const
947 if ( !setFace( faceIndex ))
950 return myCurFace.myNodeIndices;
953 //=======================================================================
954 //function : GetFaceNodes
955 //purpose : Return a set of face nodes.
956 //=======================================================================
958 bool SMDS_VolumeTool::GetFaceNodes (int faceIndex,
959 std::set<const SMDS_MeshNode*>& theFaceNodes ) const
961 if ( !setFace( faceIndex ))
964 theFaceNodes.clear();
965 theFaceNodes.insert( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
972 struct NLink : public std::pair<int,int>
975 NLink(const SMDS_MeshNode* n1=0, const SMDS_MeshNode* n2=0, int ori=1 )
979 if (( myOri = ( n1->GetID() < n2->GetID() )))
982 second = n2->GetID();
988 second = n1->GetID();
994 myOri = first = second = 0;
997 //int Node1() const { return myOri == -1 ? second : first; }
999 //bool IsSameOri( const std::pair<int,int>& link ) const { return link.first == Node1(); }
1003 //=======================================================================
1004 //function : IsFaceExternal
1005 //purpose : Check normal orientation of a given face
1006 //=======================================================================
1008 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const
1010 if ( myExternalFaces || !myVolume )
1013 if ( !myPolyedre ) // all classical volumes have external facet normals
1016 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1018 if ( myPolyFacetOri[ faceIndex ])
1019 return myPolyFacetOri[ faceIndex ] > 0;
1021 int ori = 0; // -1-in, +1-out, 0-undef
1022 double minProj, maxProj;
1023 if ( projectNodesToNormal( faceIndex, minProj, maxProj ))
1025 // all nodes are on the same side of the facet
1026 ori = ( minProj < 0 ? +1 : -1 );
1027 me->myPolyFacetOri[ faceIndex ] = ori;
1029 if ( !myFwdLinks.empty() ) // concave polyhedron; collect oriented links
1030 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1032 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1033 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1038 SaveFacet savedFacet( myCurFace );
1040 // concave polyhedron
1042 if ( myFwdLinks.empty() ) // get links of the least ambiguously oriented facet
1044 for ( size_t i = 0; i < myPolyFacetOri.size() && !ori; ++i )
1045 ori = myPolyFacetOri[ i ];
1047 if ( !ori ) // none facet is oriented yet
1049 // find the least ambiguously oriented facet
1050 int faceMostConvex = -1;
1051 std::map< double, int > convexity2face;
1052 for ( size_t iF = 0; iF < myPolyFacetOri.size() && faceMostConvex < 0; ++iF )
1054 if ( projectNodesToNormal( iF, minProj, maxProj ))
1056 // all nodes are on the same side of the facet
1057 me->myPolyFacetOri[ iF ] = ( minProj < 0 ? +1 : -1 );
1058 faceMostConvex = iF;
1062 ori = ( -minProj < maxProj ? -1 : +1 );
1063 double convexity = std::min( -minProj, maxProj ) / std::max( -minProj, maxProj );
1064 convexity2face.insert( std::make_pair( convexity, iF * ori ));
1067 if ( faceMostConvex < 0 ) // none facet has nodes on the same side
1069 // use the least ambiguous facet
1070 faceMostConvex = convexity2face.begin()->second;
1071 ori = ( faceMostConvex < 0 ? -1 : +1 );
1072 faceMostConvex = std::abs( faceMostConvex );
1073 me->myPolyFacetOri[ faceMostConvex ] = ori;
1076 // collect links of the oriented facets in myFwdLinks
1077 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1079 ori = myPolyFacetOri[ iF ];
1080 if ( !ori ) continue;
1082 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1084 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1085 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1090 // compare orientation of links of the facet with myFwdLinks
1092 setFace( faceIndex );
1093 std::vector< NLink > links( myCurFace.myNbNodes ), links2;
1094 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1096 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1097 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1098 if ( l2o != myFwdLinks.end() )
1099 ori = link.myOri * l2o->second * -1;
1102 while ( !ori ) // the facet has no common links with already oriented facets
1104 // orient and collect links of other non-oriented facets
1105 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1107 if ( myPolyFacetOri[ iF ] ) continue; // already oriented
1111 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1113 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1114 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1115 if ( l2o != myFwdLinks.end() )
1116 ori = link.myOri * l2o->second * -1;
1117 links2.push_back( link );
1119 if ( ori ) // one more facet oriented
1121 me->myPolyFacetOri[ iF ] = ori;
1122 for ( size_t i = 0; i < links2.size(); ++i )
1123 me->myFwdLinks.insert( make_pair( links2[i], links2[i].myOri * ori ));
1128 return false; // error in algorithm: infinite loop
1130 // try to orient the facet again
1132 for ( size_t i = 0; i < links.size() && !ori; ++i )
1134 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( links[i] );
1135 if ( l2o != myFwdLinks.end() )
1136 ori = links[i].myOri * l2o->second * -1;
1138 me->myPolyFacetOri[ faceIndex ] = ori;
1144 //=======================================================================
1145 //function : projectNodesToNormal
1146 //purpose : compute min and max projections of all nodes to normal of a facet.
1147 //=======================================================================
1149 bool SMDS_VolumeTool::projectNodesToNormal( int faceIndex,
1152 double* normalXYZ ) const
1154 minProj = std::numeric_limits<double>::max();
1155 maxProj = std::numeric_limits<double>::min();
1158 if ( !GetFaceNormal( faceIndex, normal.x, normal.y, normal.z ))
1161 memcpy( normalXYZ, normal.data(), 3*sizeof(double));
1163 XYZ p0 ( myCurFace.myNodes[0] );
1164 for ( size_t i = 0; i < myVolumeNodes.size(); ++i )
1166 if ( std::find( myCurFace.myNodes.begin() + 1,
1167 myCurFace.myNodes.end(),
1168 myVolumeNodes[ i ] ) != myCurFace.myNodes.end() )
1169 continue; // node of the faceIndex-th facet
1171 double proj = normal.Dot( XYZ( myVolumeNodes[ i ]) - p0 );
1172 if ( proj < minProj ) minProj = proj;
1173 if ( proj > maxProj ) maxProj = proj;
1175 const double tol = 1e-7;
1178 bool diffSize = ( minProj * maxProj < 0 );
1181 // minProj = -minProj;
1183 // else if ( minProj < 0 )
1185 // minProj = -minProj;
1186 // maxProj = -maxProj;
1189 return !diffSize; // ? 0 : (minProj >= 0);
1192 //=======================================================================
1193 //function : GetFaceNormal
1194 //purpose : Return a normal to a face
1195 //=======================================================================
1197 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) const
1199 if ( !setFace( faceIndex ))
1202 const int iQuad = ( !myPolyedre && myCurFace.myNbNodes > 6 ) ? 2 : 1;
1203 XYZ p1 ( myCurFace.myNodes[0*iQuad] );
1204 XYZ p2 ( myCurFace.myNodes[1*iQuad] );
1205 XYZ p3 ( myCurFace.myNodes[2*iQuad] );
1206 XYZ aVec12( p2 - p1 );
1207 XYZ aVec13( p3 - p1 );
1208 XYZ cross = aVec12.Crossed( aVec13 );
1210 for ( int i = 3*iQuad; i < myCurFace.myNbNodes; i += iQuad )
1212 XYZ p4 ( myCurFace.myNodes[i] );
1213 XYZ aVec14( p4 - p1 );
1214 XYZ cross2 = aVec13.Crossed( aVec14 );
1215 cross = cross + cross2;
1219 double size = cross.Magnitude();
1220 if ( size <= std::numeric_limits<double>::min() )
1230 //================================================================================
1232 * \brief Return barycenter of a face
1234 //================================================================================
1236 bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const
1238 if ( !setFace( faceIndex ))
1242 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1244 X += myCurFace.myNodes[i]->X() / myCurFace.myNbNodes;
1245 Y += myCurFace.myNodes[i]->Y() / myCurFace.myNbNodes;
1246 Z += myCurFace.myNodes[i]->Z() / myCurFace.myNbNodes;
1251 //=======================================================================
1252 //function : GetFaceArea
1253 //purpose : Return face area
1254 //=======================================================================
1256 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
1259 if ( !setFace( faceIndex ))
1262 XYZ p1 ( myCurFace.myNodes[0] );
1263 XYZ p2 ( myCurFace.myNodes[1] );
1264 XYZ p3 ( myCurFace.myNodes[2] );
1265 XYZ aVec12( p2 - p1 );
1266 XYZ aVec13( p3 - p1 );
1267 area += aVec12.Crossed( aVec13 ).Magnitude();
1269 if (myVolume->IsPoly())
1271 for ( int i = 3; i < myCurFace.myNbNodes; ++i )
1273 XYZ pI ( myCurFace.myNodes[i] );
1274 XYZ aVecI( pI - p1 );
1275 area += aVec13.Crossed( aVecI ).Magnitude();
1281 if ( myCurFace.myNbNodes == 4 ) {
1282 XYZ p4 ( myCurFace.myNodes[3] );
1283 XYZ aVec14( p4 - p1 );
1284 area += aVec14.Crossed( aVec13 ).Magnitude();
1290 //================================================================================
1292 * \brief Return index of the node located at face center of a quadratic element like HEX27
1294 //================================================================================
1296 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
1298 if ( myAllFacesNbNodes && myVolumeNodes.size() == 27 ) // classic element with 27 nodes
1300 switch ( faceIndex ) {
1304 return faceIndex + 19;
1310 //=======================================================================
1311 //function : GetOppFaceIndex
1312 //purpose : Return index of the opposite face if it exists, else -1.
1313 //=======================================================================
1315 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1319 MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1323 const int nbHoriFaces = 2;
1325 if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1326 switch ( myVolumeNodes.size() ) {
1329 if ( faceIndex == 0 || faceIndex == 1 )
1330 ind = 1 - faceIndex;
1334 if ( faceIndex <= 1 ) // top or bottom
1335 ind = 1 - faceIndex;
1337 const int nbSideFaces = myAllFacesNbNodes[0];
1338 ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1343 ind = GetOppFaceIndexOfHex( faceIndex );
1351 //=======================================================================
1352 //function : GetOppFaceIndexOfHex
1353 //purpose : Return index of the opposite face of the hexahedron
1354 //=======================================================================
1356 int SMDS_VolumeTool::GetOppFaceIndexOfHex( int faceIndex )
1358 return Hexa_oppF[ faceIndex ];
1361 //=======================================================================
1362 //function : IsLinked
1363 //purpose : return true if theNode1 is linked with theNode2
1364 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1365 //=======================================================================
1367 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1368 const SMDS_MeshNode* theNode2,
1369 const bool theIgnoreMediumNodes) const
1374 if (myVolume->IsPoly()) {
1376 MESSAGE("Warning: bad volumic element");
1379 if ( !myAllFacesNbNodes ) {
1380 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1381 me->myPolyQuantities = myPolyedre->GetQuantities();
1382 myAllFacesNbNodes = &myPolyQuantities[0];
1384 int from, to = 0, d1 = 1, d2 = 2;
1385 if ( myPolyedre->IsQuadratic() ) {
1386 if ( theIgnoreMediumNodes ) {
1392 std::vector<const SMDS_MeshNode*>::const_iterator i;
1393 for (int iface = 0; iface < myNbFaces; iface++)
1396 to += myPolyQuantities[iface];
1397 i = std::find( myVolumeNodes.begin() + from, myVolumeNodes.begin() + to, theNode1 );
1398 if ( i != myVolumeNodes.end() )
1400 if (( theNode2 == *( i-d1 ) ||
1401 theNode2 == *( i+d1 )))
1404 (( theNode2 == *( i-d2 ) ||
1405 theNode2 == *( i+d2 ))))
1412 // find nodes indices
1413 int i1 = -1, i2 = -1, nbFound = 0;
1414 for ( size_t i = 0; i < myVolumeNodes.size() && nbFound < 2; i++ )
1416 if ( myVolumeNodes[ i ] == theNode1 )
1418 else if ( myVolumeNodes[ i ] == theNode2 )
1421 return IsLinked( i1, i2 );
1424 //=======================================================================
1425 //function : IsLinked
1426 //purpose : return true if the node with theNode1Index is linked
1427 // with the node with theNode2Index
1428 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1429 //=======================================================================
1431 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1432 const int theNode2Index,
1433 bool theIgnoreMediumNodes) const
1435 if ( myVolume->IsPoly() ) {
1436 return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1439 int minInd = std::min( theNode1Index, theNode2Index );
1440 int maxInd = std::max( theNode1Index, theNode2Index );
1442 if ( minInd < 0 || maxInd > (int)myVolumeNodes.size() - 1 || maxInd == minInd )
1445 VolumeType type = GetVolumeType();
1446 if ( myVolume->IsQuadratic() )
1448 int firstMediumInd = myVolume->NbCornerNodes();
1449 if ( minInd >= firstMediumInd )
1450 return false; // both nodes are medium - not linked
1451 if ( maxInd < firstMediumInd ) // both nodes are corners
1453 if ( theIgnoreMediumNodes )
1454 type = quadToLinear(type); // to check linkage of corner nodes only
1456 return false; // corner nodes are not linked directly in a quadratic cell
1464 switch ( maxInd - minInd ) {
1465 case 1: return minInd != 3;
1466 case 3: return minInd == 0 || minInd == 4;
1467 case 4: return true;
1474 switch ( maxInd - minInd ) {
1476 case 3: return true;
1481 switch ( maxInd - minInd ) {
1482 case 1: return minInd != 2;
1483 case 2: return minInd == 0 || minInd == 3;
1484 case 3: return true;
1491 case 0: return ( maxInd==4 || maxInd==6 || maxInd==7 );
1492 case 1: return ( maxInd==4 || maxInd==5 || maxInd==8 );
1493 case 2: return ( maxInd==5 || maxInd==6 || maxInd==9 );
1494 case 3: return ( maxInd==7 || maxInd==8 || maxInd==9 );
1502 case 0: return ( maxInd==8 || maxInd==11 || maxInd==16 );
1503 case 1: return ( maxInd==8 || maxInd==9 || maxInd==17 );
1504 case 2: return ( maxInd==9 || maxInd==10 || maxInd==18 );
1505 case 3: return ( maxInd==10 || maxInd==11 || maxInd==19 );
1506 case 4: return ( maxInd==12 || maxInd==15 || maxInd==16 );
1507 case 5: return ( maxInd==12 || maxInd==13 || maxInd==17 );
1508 case 6: return ( maxInd==13 || maxInd==14 || maxInd==18 );
1509 case 7: return ( maxInd==14 || maxInd==15 || maxInd==19 );
1517 case 0: return ( maxInd==5 || maxInd==8 || maxInd==9 );
1518 case 1: return ( maxInd==5 || maxInd==6 || maxInd==10 );
1519 case 2: return ( maxInd==6 || maxInd==7 || maxInd==11 );
1520 case 3: return ( maxInd==7 || maxInd==8 || maxInd==12 );
1521 case 4: return ( maxInd==9 || maxInd==10 || maxInd==11 || maxInd==12 );
1529 case 0: return ( maxInd==6 || maxInd==8 || maxInd==12 );
1530 case 1: return ( maxInd==6 || maxInd==7 || maxInd==13 );
1531 case 2: return ( maxInd==7 || maxInd==8 || maxInd==14 );
1532 case 3: return ( maxInd==9 || maxInd==11 || maxInd==12 );
1533 case 4: return ( maxInd==9 || maxInd==10 || maxInd==13 );
1534 case 5: return ( maxInd==10 || maxInd==11 || maxInd==14 );
1541 const int diff = maxInd-minInd;
1542 if ( diff > 6 ) return false;// not linked top and bottom
1543 if ( diff == 6 ) return true; // linked top and bottom
1544 return diff == 1 || diff == 7;
1551 //=======================================================================
1552 //function : GetNodeIndex
1553 //purpose : Return an index of theNode
1554 //=======================================================================
1556 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1559 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
1560 if ( myVolumeNodes[ i ] == theNode )
1567 //================================================================================
1569 * \brief Fill vector with boundary faces existing in the mesh
1570 * \param faces - vector of found nodes
1571 * \retval int - nb of found faces
1573 //================================================================================
1575 int SMDS_VolumeTool::GetAllExistingFaces(std::vector<const SMDS_MeshElement*> & faces) const
1578 SaveFacet savedFacet( myCurFace );
1580 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1582 if ( const SMDS_MeshElement* face = SMDS_Mesh::FindFace( myCurFace.myNodes ))
1583 faces.push_back( face );
1586 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1587 const SMDS_MeshFace* face = 0;
1588 const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1589 switch ( NbFaceNodes( iF )) {
1591 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1593 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1595 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1596 nodes[3], nodes[4], nodes[5]); break;
1598 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1599 nodes[4], nodes[5], nodes[6], nodes[7]); break;
1602 faces.push_back( face );
1604 return faces.size();
1608 //================================================================================
1610 * \brief Fill vector with boundary edges existing in the mesh
1611 * \param edges - vector of found edges
1612 * \retval int - nb of found faces
1614 //================================================================================
1616 int SMDS_VolumeTool::GetAllExistingEdges(std::vector<const SMDS_MeshElement*> & edges) const
1619 edges.reserve( myVolumeNodes.size() * 2 );
1620 for ( size_t i = 0; i < myVolumeNodes.size()-1; ++i ) {
1621 for ( size_t j = i + 1; j < myVolumeNodes.size(); ++j ) {
1622 if ( IsLinked( i, j )) {
1623 const SMDS_MeshElement* edge =
1624 SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1626 edges.push_back( edge );
1630 return edges.size();
1633 //================================================================================
1635 * \brief Return minimal square distance between connected corner nodes
1637 //================================================================================
1639 double SMDS_VolumeTool::MinLinearSize2() const
1641 double minSize = 1e+100;
1642 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1644 SaveFacet savedFacet( myCurFace );
1646 // it seems that compute distance twice is faster than organization of a sole computing
1647 myCurFace.myIndex = -1;
1648 for ( int iF = 0; iF < myNbFaces; ++iF )
1651 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1653 XYZ n1( myCurFace.myNodes[ iN ]);
1654 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1655 minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
1662 //================================================================================
1664 * \brief Return maximal square distance between connected corner nodes
1666 //================================================================================
1668 double SMDS_VolumeTool::MaxLinearSize2() const
1670 double maxSize = -1e+100;
1671 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1673 SaveFacet savedFacet( myCurFace );
1675 // it seems that compute distance twice is faster than organization of a sole computing
1676 myCurFace.myIndex = -1;
1677 for ( int iF = 0; iF < myNbFaces; ++iF )
1680 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1682 XYZ n1( myCurFace.myNodes[ iN ]);
1683 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1684 maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude());
1691 //================================================================================
1693 * \brief Fast quickly check that only one volume is built on the face nodes
1694 * This check is valid for conformal meshes only
1696 //================================================================================
1698 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1700 const bool isFree = true;
1702 if ( !setFace( faceIndex ))
1705 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1707 const int di = myVolume->IsQuadratic() ? 2 : 1;
1708 const int nbN = ( myCurFace.myNbNodes/di <= 4 && !IsPoly()) ? 3 : myCurFace.myNbNodes/di; // nb nodes to check
1710 SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
1711 while ( eIt->more() )
1713 const SMDS_MeshElement* vol = eIt->next();
1714 if ( vol == myVolume )
1717 for ( iN = 1; iN < nbN; ++iN )
1718 if ( vol->GetNodeIndex( nodes[ iN*di ]) < 0 )
1720 if ( iN == nbN ) // nbN nodes are shared with vol
1722 // if ( vol->IsPoly() || vol->NbFaces() > 6 ) // vol is polyhed or hex prism
1724 // int nb = myCurFace.myNbNodes;
1725 // if ( myVolume->GetEntityType() != vol->GetEntityType() )
1726 // nb -= ( GetCenterNodeIndex(0) > 0 );
1727 // std::set<const SMDS_MeshNode*> faceNodes( nodes, nodes + nb );
1728 // if ( SMDS_VolumeTool( vol ).GetFaceIndex( faceNodes ) < 0 )
1731 if ( otherVol ) *otherVol = vol;
1735 if ( otherVol ) *otherVol = 0;
1739 //================================================================================
1741 * \brief Thorough check that only one volume is built on the face nodes
1743 //================================================================================
1745 bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1747 const bool isFree = true;
1749 if (!setFace( faceIndex ))
1752 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1753 const int nbFaceNodes = myCurFace.myNbNodes;
1755 // evaluate nb of face nodes shared by other volumes
1756 int maxNbShared = -1;
1757 typedef std::map< const SMDS_MeshElement*, int > TElemIntMap;
1758 TElemIntMap volNbShared;
1759 TElemIntMap::iterator vNbIt;
1760 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1761 const SMDS_MeshNode* n = nodes[ iNode ];
1762 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1763 while ( eIt->more() ) {
1764 const SMDS_MeshElement* elem = eIt->next();
1765 if ( elem != myVolume ) {
1766 vNbIt = volNbShared.insert( std::make_pair( elem, 0 )).first;
1768 if ( vNbIt->second > maxNbShared )
1769 maxNbShared = vNbIt->second;
1773 if ( maxNbShared < 3 )
1774 return isFree; // is free
1776 // find volumes laying on the opposite side of the face
1777 // and sharing all nodes
1778 XYZ intNormal; // internal normal
1779 GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
1780 if ( IsFaceExternal( faceIndex ))
1781 intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
1782 XYZ p0 ( nodes[0] ), baryCenter;
1783 for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); ) {
1784 const int& nbShared = (*vNbIt).second;
1785 if ( nbShared >= 3 ) {
1786 SMDS_VolumeTool volume( (*vNbIt).first );
1787 volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
1788 XYZ intNormal2( baryCenter - p0 );
1789 if ( intNormal.Dot( intNormal2 ) < 0 ) {
1791 if ( nbShared >= nbFaceNodes )
1793 // a volume shares the whole facet
1794 if ( otherVol ) *otherVol = vNbIt->first;
1801 // remove a volume from volNbShared map
1802 volNbShared.erase( vNbIt++ );
1805 // here volNbShared contains only volumes laying on the opposite side of
1806 // the face and sharing 3 or more but not all face nodes with myVolume
1807 if ( volNbShared.size() < 2 ) {
1808 return isFree; // is free
1811 // check if the whole area of a face is shared
1812 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1814 const SMDS_MeshNode* n = nodes[ iNode ];
1815 // check if n is shared by one of volumes of volNbShared
1816 bool isShared = false;
1817 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1818 while ( eIt->more() && !isShared )
1819 isShared = volNbShared.count( eIt->next() );
1823 if ( otherVol ) *otherVol = volNbShared.begin()->first;
1826 // if ( !myVolume->IsPoly() )
1828 // bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
1829 // for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1830 // SMDS_VolumeTool volume( (*vNbIt).first );
1831 // bool prevLinkShared = false;
1832 // int nbSharedLinks = 0;
1833 // for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1834 // bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
1835 // if ( linkShared )
1837 // if ( linkShared && prevLinkShared &&
1838 // volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
1839 // isShared[ iNode ] = true;
1840 // prevLinkShared = linkShared;
1842 // if ( nbSharedLinks == nbFaceNodes )
1843 // return !free; // is not free
1844 // if ( nbFaceNodes == 4 ) {
1845 // // check traingle parts 1 & 3
1846 // if ( isShared[1] && isShared[3] )
1847 // return !free; // is not free
1848 // // check triangle parts 0 & 2;
1849 // // 0 part could not be checked in the loop; check it here
1850 // if ( isShared[2] && prevLinkShared &&
1851 // volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
1852 // volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
1853 // return !free; // is not free
1860 //=======================================================================
1861 //function : GetFaceIndex
1862 //purpose : Return index of a face formed by theFaceNodes
1863 //=======================================================================
1865 int SMDS_VolumeTool::GetFaceIndex( const std::set<const SMDS_MeshNode*>& theFaceNodes,
1866 const int theFaceIndexHint ) const
1868 if ( theFaceIndexHint >= 0 )
1870 int nbNodes = NbFaceNodes( theFaceIndexHint );
1871 if ( nbNodes == (int) theFaceNodes.size() )
1873 const SMDS_MeshNode** nodes = GetFaceNodes( theFaceIndexHint );
1875 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1880 return theFaceIndexHint;
1883 for ( int iFace = 0; iFace < myNbFaces; iFace++ )
1885 if ( iFace == theFaceIndexHint )
1887 int nbNodes = NbFaceNodes( iFace );
1888 if ( nbNodes == (int) theFaceNodes.size() )
1890 const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
1892 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1903 //=======================================================================
1904 //function : GetFaceIndex
1905 //purpose : Return index of a face formed by theFaceNodes
1906 //=======================================================================
1908 /*int SMDS_VolumeTool::GetFaceIndex( const std::set<int>& theFaceNodesIndices )
1910 for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1911 const int* nodes = GetFaceNodesIndices( iFace );
1912 int nbFaceNodes = NbFaceNodes( iFace );
1913 std::set<int> nodeSet;
1914 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1915 nodeSet.insert( nodes[ iNode ] );
1916 if ( theFaceNodesIndices == nodeSet )
1922 //=======================================================================
1923 //function : setFace
1925 //=======================================================================
1927 bool SMDS_VolumeTool::setFace( int faceIndex ) const
1932 if ( myCurFace.myIndex == faceIndex )
1935 myCurFace.myIndex = -1;
1937 if ( faceIndex < 0 || faceIndex >= NbFaces() )
1940 if (myVolume->IsPoly())
1942 if ( !myPolyedre ) {
1943 MESSAGE("Warning: bad volumic element");
1948 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1949 if ( !myAllFacesNbNodes ) {
1950 me->myPolyQuantities = myPolyedre->GetQuantities();
1951 myAllFacesNbNodes = &myPolyQuantities[0];
1953 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
1954 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
1955 me->myPolyIndices.resize( myCurFace.myNbNodes + 1 );
1956 myCurFace.myNodeIndices = & me->myPolyIndices[0];
1957 int shift = std::accumulate( myAllFacesNbNodes, myAllFacesNbNodes+faceIndex, 0 );
1958 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
1960 myCurFace.myNodes [ iNode ] = myVolumeNodes[ shift + iNode ];
1961 myCurFace.myNodeIndices[ iNode ] = shift + iNode;
1963 myCurFace.myNodes [ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ]; // last = first
1964 myCurFace.myNodeIndices[ myCurFace.myNbNodes ] = myCurFace.myNodeIndices[ 0 ];
1966 // check orientation
1967 if (myExternalFaces)
1969 myCurFace.myIndex = faceIndex; // avoid infinite recursion in IsFaceExternal()
1970 myExternalFaces = false; // force normal computation by IsFaceExternal()
1971 if ( !IsFaceExternal( faceIndex ))
1972 std::reverse( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
1973 myExternalFaces = true;
1978 if ( !myAllFacesNodeIndices_F )
1980 // choose data for an element type
1981 switch ( myVolumeNodes.size() ) {
1983 myAllFacesNodeIndices_F = &Tetra_F [0][0];
1984 //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
1985 myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
1986 myAllFacesNbNodes = Tetra_nbN;
1987 myMaxFaceNbNodes = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
1990 myAllFacesNodeIndices_F = &Pyramid_F [0][0];
1991 //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
1992 myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
1993 myAllFacesNbNodes = Pyramid_nbN;
1994 myMaxFaceNbNodes = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
1997 myAllFacesNodeIndices_F = &Penta_F [0][0];
1998 //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
1999 myAllFacesNodeIndices_RE = &Penta_RE[0][0];
2000 myAllFacesNbNodes = Penta_nbN;
2001 myMaxFaceNbNodes = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
2004 myAllFacesNodeIndices_F = &Hexa_F [0][0];
2005 ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
2006 myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
2007 myAllFacesNbNodes = Hexa_nbN;
2008 myMaxFaceNbNodes = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
2011 myAllFacesNodeIndices_F = &QuadTetra_F [0][0];
2012 //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
2013 myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
2014 myAllFacesNbNodes = QuadTetra_nbN;
2015 myMaxFaceNbNodes = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
2018 myAllFacesNodeIndices_F = &QuadPyram_F [0][0];
2019 //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
2020 myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
2021 myAllFacesNbNodes = QuadPyram_nbN;
2022 myMaxFaceNbNodes = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
2025 myAllFacesNodeIndices_F = &QuadPenta_F [0][0];
2026 //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
2027 myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
2028 myAllFacesNbNodes = QuadPenta_nbN;
2029 myMaxFaceNbNodes = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
2033 myAllFacesNodeIndices_F = &QuadHexa_F [0][0];
2034 //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
2035 myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
2036 myAllFacesNbNodes = QuadHexa_nbN;
2037 myMaxFaceNbNodes = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
2038 if ( !myIgnoreCentralNodes && myVolumeNodes.size() == 27 )
2040 myAllFacesNodeIndices_F = &TriQuadHexa_F [0][0];
2041 //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
2042 myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
2043 myAllFacesNbNodes = TriQuadHexa_nbN;
2044 myMaxFaceNbNodes = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
2048 myAllFacesNodeIndices_F = &HexPrism_F [0][0];
2049 //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
2050 myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
2051 myAllFacesNbNodes = HexPrism_nbN;
2052 myMaxFaceNbNodes = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
2058 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2059 // if ( myExternalFaces )
2060 // myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2062 // myCurFace.myNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
2063 myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2066 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2067 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2068 myCurFace.myNodes[ iNode ] = myVolumeNodes[ myCurFace.myNodeIndices[ iNode ]];
2069 myCurFace.myNodes[ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ];
2072 myCurFace.myIndex = faceIndex;
2077 //=======================================================================
2078 //function : GetType
2079 //purpose : return VolumeType by nb of nodes in a volume
2080 //=======================================================================
2082 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
2084 switch ( nbNodes ) {
2085 case 4: return TETRA;
2086 case 5: return PYRAM;
2087 case 6: return PENTA;
2088 case 8: return HEXA;
2089 case 10: return QUAD_TETRA;
2090 case 13: return QUAD_PYRAM;
2091 case 15: return QUAD_PENTA;
2093 case 27: return QUAD_HEXA;
2094 case 12: return HEX_PRISM;
2095 default:return UNKNOWN;
2099 //=======================================================================
2100 //function : NbFaces
2101 //purpose : return nb of faces by volume type
2102 //=======================================================================
2104 int SMDS_VolumeTool::NbFaces( VolumeType type )
2108 case QUAD_TETRA: return 4;
2110 case QUAD_PYRAM: return 5;
2112 case QUAD_PENTA: return 5;
2114 case QUAD_HEXA : return 6;
2115 case HEX_PRISM : return 8;
2120 //================================================================================
2122 * \brief Useful to know nb of corner nodes of a quadratic volume
2123 * \param type - volume type
2124 * \retval int - nb of corner nodes
2126 //================================================================================
2128 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
2132 case QUAD_TETRA: return 4;
2134 case QUAD_PYRAM: return 5;
2136 case QUAD_PENTA: return 6;
2138 case QUAD_HEXA : return 8;
2139 case HEX_PRISM : return 12;
2146 //=======================================================================
2147 //function : GetFaceNodesIndices
2148 //purpose : Return the array of face nodes indices
2149 // To comfort link iteration, the array
2150 // length == NbFaceNodes( faceIndex ) + 1 and
2151 // the last node index == the first one.
2152 //=======================================================================
2154 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
2159 case TETRA: return Tetra_F[ faceIndex ];
2160 case PYRAM: return Pyramid_F[ faceIndex ];
2161 case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
2162 case HEXA: return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
2163 case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
2164 case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
2165 case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
2166 // what about SMDSEntity_TriQuad_Hexa?
2167 case QUAD_HEXA: return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
2168 case HEX_PRISM: return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
2174 //=======================================================================
2175 //function : NbFaceNodes
2176 //purpose : Return number of nodes in the array of face nodes
2177 //=======================================================================
2179 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
2183 case TETRA: return Tetra_nbN[ faceIndex ];
2184 case PYRAM: return Pyramid_nbN[ faceIndex ];
2185 case PENTA: return Penta_nbN[ faceIndex ];
2186 case HEXA: return Hexa_nbN[ faceIndex ];
2187 case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
2188 case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
2189 case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
2190 // what about SMDSEntity_TriQuad_Hexa?
2191 case QUAD_HEXA: return QuadHexa_nbN[ faceIndex ];
2192 case HEX_PRISM: return HexPrism_nbN[ faceIndex ];
2198 //=======================================================================
2199 //function : Element
2200 //purpose : return element
2201 //=======================================================================
2203 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
2205 return static_cast<const SMDS_MeshVolume*>( myVolume );
2208 //=======================================================================
2210 //purpose : return element ID
2211 //=======================================================================
2213 int SMDS_VolumeTool::ID() const
2215 return myVolume ? myVolume->GetID() : 0;