1 // Copyright (C) 2007-2019 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 // File : SMDS_VolumeTool.cxx
24 // Created : Tue Jul 13 12:22:13 2004
25 // Author : Edward AGAPOV (eap)
28 #pragma warning(disable:4786)
31 #include "SMDS_VolumeTool.hxx"
33 #include "SMDS_MeshElement.hxx"
34 #include "SMDS_MeshNode.hxx"
35 #include "SMDS_Mesh.hxx"
37 #include <utilities.h>
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 // split a polyhedron into tetrahedrons
711 SaveFacet savedFacet( myCurFace );
712 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
713 XYZ origin( 1 + 1e-6, 22 + 2e-6, 333 + 3e-6 ); // for invalid poly: avoid lying on a facet plane
714 for ( int f = 0; f < NbFaces(); ++f )
717 XYZ area (0,0,0), p1( myCurFace.myNodes[0] );
718 for ( int n = 0; n < myCurFace.myNbNodes; ++n )
720 XYZ p2( myCurFace.myNodes[ n+1 ]);
721 area = area + p1.Crossed( p2 );
724 V += ( p1 - origin ).Dot( area );
730 const static int ind[] = {
731 0, 1, 3, 6, 11, 23, 31, 44, 58, 78 };
732 const static int vtab[][4] = { // decomposition into tetra in the order of enum VolumeType
765 // quadratic tetrahedron
790 // quadratic pentahedron
807 // quadratic hexahedron
832 int type = GetVolumeType();
834 int n2 = ind[type+1];
836 for (int i = n1; i < n2; i++) {
837 V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
838 myVolumeNodes[ vtab[i][1] ],
839 myVolumeNodes[ vtab[i][2] ],
840 myVolumeNodes[ vtab[i][3] ]);
846 //=======================================================================
847 //function : GetBaryCenter
849 //=======================================================================
851 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
857 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
858 X += myVolumeNodes[ i ]->X();
859 Y += myVolumeNodes[ i ]->Y();
860 Z += myVolumeNodes[ i ]->Z();
862 X /= myVolumeNodes.size();
863 Y /= myVolumeNodes.size();
864 Z /= myVolumeNodes.size();
869 //================================================================================
871 * \brief Classify a point
872 * \param tol - thickness of faces
874 //================================================================================
876 bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
878 // LIMITATION: for convex volumes only
880 for ( int iF = 0; iF < myNbFaces; ++iF )
883 if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z ))
885 if ( !IsFaceExternal( iF ))
886 faceNormal = XYZ() - faceNormal; // reverse
888 XYZ face2p( p - XYZ( myCurFace.myNodes[0] ));
889 if ( face2p.Dot( faceNormal ) > tol )
895 //=======================================================================
896 //function : SetExternalNormal
897 //purpose : Node order will be so that faces normals are external
898 //=======================================================================
900 void SMDS_VolumeTool::SetExternalNormal ()
902 myExternalFaces = true;
903 myCurFace.myIndex = -1;
906 //=======================================================================
907 //function : NbFaceNodes
908 //purpose : Return number of nodes in the array of face nodes
909 //=======================================================================
911 int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const
913 if ( !setFace( faceIndex ))
915 return myCurFace.myNbNodes;
918 //=======================================================================
919 //function : GetFaceNodes
920 //purpose : Return pointer to the array of face nodes.
921 // To comfort link iteration, the array
922 // length == NbFaceNodes( faceIndex ) + 1 and
923 // the last node == the first one.
924 //=======================================================================
926 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const
928 if ( !setFace( faceIndex ))
930 return &myCurFace.myNodes[0];
933 //=======================================================================
934 //function : GetFaceNodesIndices
935 //purpose : Return pointer to the array of face nodes indices
936 // To comfort link iteration, the array
937 // length == NbFaceNodes( faceIndex ) + 1 and
938 // the last node index == the first one.
939 //=======================================================================
941 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const
943 if ( !setFace( faceIndex ))
946 return myCurFace.myNodeIndices;
949 //=======================================================================
950 //function : GetFaceNodes
951 //purpose : Return a set of face nodes.
952 //=======================================================================
954 bool SMDS_VolumeTool::GetFaceNodes (int faceIndex,
955 std::set<const SMDS_MeshNode*>& theFaceNodes ) const
957 if ( !setFace( faceIndex ))
960 theFaceNodes.clear();
961 theFaceNodes.insert( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
968 struct NLink : public std::pair<int,int>
971 NLink(const SMDS_MeshNode* n1=0, const SMDS_MeshNode* n2=0, int ori=1 )
975 if (( myOri = ( n1->GetID() < n2->GetID() )))
978 second = n2->GetID();
984 second = n1->GetID();
990 myOri = first = second = 0;
993 //int Node1() const { return myOri == -1 ? second : first; }
995 //bool IsSameOri( const std::pair<int,int>& link ) const { return link.first == Node1(); }
999 //=======================================================================
1000 //function : IsFaceExternal
1001 //purpose : Check normal orientation of a given face
1002 //=======================================================================
1004 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const
1006 if ( myExternalFaces || !myVolume )
1009 if ( !myPolyedre ) // all classical volumes have external facet normals
1012 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1014 if ( myPolyFacetOri[ faceIndex ])
1015 return myPolyFacetOri[ faceIndex ] > 0;
1017 int ori = 0; // -1-in, +1-out, 0-undef
1018 double minProj, maxProj;
1019 if ( projectNodesToNormal( faceIndex, minProj, maxProj ))
1021 // all nodes are on the same side of the facet
1022 ori = ( minProj < 0 ? +1 : -1 );
1023 me->myPolyFacetOri[ faceIndex ] = ori;
1025 if ( !myFwdLinks.empty() ) // concave polyhedron; collect oriented links
1026 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1028 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1029 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1034 SaveFacet savedFacet( myCurFace );
1036 // concave polyhedron
1038 if ( myFwdLinks.empty() ) // get links of the least ambiguously oriented facet
1040 for ( size_t i = 0; i < myPolyFacetOri.size() && !ori; ++i )
1041 ori = myPolyFacetOri[ i ];
1043 if ( !ori ) // none facet is oriented yet
1045 // find the least ambiguously oriented facet
1046 int faceMostConvex = -1;
1047 std::map< double, int > convexity2face;
1048 for ( size_t iF = 0; iF < myPolyFacetOri.size() && faceMostConvex < 0; ++iF )
1050 if ( projectNodesToNormal( iF, minProj, maxProj ))
1052 // all nodes are on the same side of the facet
1053 me->myPolyFacetOri[ iF ] = ( minProj < 0 ? +1 : -1 );
1054 faceMostConvex = iF;
1058 ori = ( -minProj < maxProj ? -1 : +1 );
1059 double convexity = std::min( -minProj, maxProj ) / std::max( -minProj, maxProj );
1060 convexity2face.insert( std::make_pair( convexity, iF * ori ));
1063 if ( faceMostConvex < 0 ) // none facet has nodes on the same side
1065 // use the least ambiguous facet
1066 faceMostConvex = convexity2face.begin()->second;
1067 ori = ( faceMostConvex < 0 ? -1 : +1 );
1068 faceMostConvex = std::abs( faceMostConvex );
1069 me->myPolyFacetOri[ faceMostConvex ] = ori;
1072 // collect links of the oriented facets in myFwdLinks
1073 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1075 ori = myPolyFacetOri[ iF ];
1076 if ( !ori ) continue;
1078 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1080 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1081 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1086 // compare orientation of links of the facet with myFwdLinks
1088 setFace( faceIndex );
1089 std::vector< NLink > links( myCurFace.myNbNodes ), links2;
1090 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1092 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1093 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1094 if ( l2o != myFwdLinks.end() )
1095 ori = link.myOri * l2o->second * -1;
1098 while ( !ori ) // the facet has no common links with already oriented facets
1100 // orient and collect links of other non-oriented facets
1101 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1103 if ( myPolyFacetOri[ iF ] ) continue; // already oriented
1107 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1109 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1110 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1111 if ( l2o != myFwdLinks.end() )
1112 ori = link.myOri * l2o->second * -1;
1113 links2.push_back( link );
1115 if ( ori ) // one more facet oriented
1117 me->myPolyFacetOri[ iF ] = ori;
1118 for ( size_t i = 0; i < links2.size(); ++i )
1119 me->myFwdLinks.insert( make_pair( links2[i], links2[i].myOri * ori ));
1124 return false; // error in algorithm: infinite loop
1126 // try to orient the facet again
1128 for ( size_t i = 0; i < links.size() && !ori; ++i )
1130 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( links[i] );
1131 if ( l2o != myFwdLinks.end() )
1132 ori = links[i].myOri * l2o->second * -1;
1134 me->myPolyFacetOri[ faceIndex ] = ori;
1140 //=======================================================================
1141 //function : projectNodesToNormal
1142 //purpose : compute min and max projections of all nodes to normal of a facet.
1143 //=======================================================================
1145 bool SMDS_VolumeTool::projectNodesToNormal( int faceIndex,
1148 double* normalXYZ ) const
1150 minProj = std::numeric_limits<double>::max();
1151 maxProj = std::numeric_limits<double>::min();
1154 if ( !GetFaceNormal( faceIndex, normal.x, normal.y, normal.z ))
1157 memcpy( normalXYZ, normal.data(), 3*sizeof(double));
1159 XYZ p0 ( myCurFace.myNodes[0] );
1160 for ( size_t i = 0; i < myVolumeNodes.size(); ++i )
1162 if ( std::find( myCurFace.myNodes.begin() + 1,
1163 myCurFace.myNodes.end(),
1164 myVolumeNodes[ i ] ) != myCurFace.myNodes.end() )
1165 continue; // node of the faceIndex-th facet
1167 double proj = normal.Dot( XYZ( myVolumeNodes[ i ]) - p0 );
1168 if ( proj < minProj ) minProj = proj;
1169 if ( proj > maxProj ) maxProj = proj;
1171 const double tol = 1e-7;
1174 bool diffSize = ( minProj * maxProj < 0 );
1177 // minProj = -minProj;
1179 // else if ( minProj < 0 )
1181 // minProj = -minProj;
1182 // maxProj = -maxProj;
1185 return !diffSize; // ? 0 : (minProj >= 0);
1188 //=======================================================================
1189 //function : GetFaceNormal
1190 //purpose : Return a normal to a face
1191 //=======================================================================
1193 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) const
1195 if ( !setFace( faceIndex ))
1198 const int iQuad = ( !myPolyedre && myCurFace.myNbNodes > 6 ) ? 2 : 1;
1199 XYZ p1 ( myCurFace.myNodes[0*iQuad] );
1200 XYZ p2 ( myCurFace.myNodes[1*iQuad] );
1201 XYZ p3 ( myCurFace.myNodes[2*iQuad] );
1202 XYZ aVec12( p2 - p1 );
1203 XYZ aVec13( p3 - p1 );
1204 XYZ cross = aVec12.Crossed( aVec13 );
1206 if ( myCurFace.myNbNodes >3*iQuad ) {
1207 XYZ p4 ( myCurFace.myNodes[3*iQuad] );
1208 XYZ aVec14( p4 - p1 );
1209 XYZ cross2 = aVec13.Crossed( aVec14 );
1210 cross = cross + cross2;
1213 double size = cross.Magnitude();
1214 if ( size <= std::numeric_limits<double>::min() )
1224 //================================================================================
1226 * \brief Return barycenter of a face
1228 //================================================================================
1230 bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const
1232 if ( !setFace( faceIndex ))
1236 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1238 X += myCurFace.myNodes[i]->X() / myCurFace.myNbNodes;
1239 Y += myCurFace.myNodes[i]->Y() / myCurFace.myNbNodes;
1240 Z += myCurFace.myNodes[i]->Z() / myCurFace.myNbNodes;
1245 //=======================================================================
1246 //function : GetFaceArea
1247 //purpose : Return face area
1248 //=======================================================================
1250 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
1253 if ( !setFace( faceIndex ))
1256 XYZ p1 ( myCurFace.myNodes[0] );
1257 XYZ p2 ( myCurFace.myNodes[1] );
1258 XYZ p3 ( myCurFace.myNodes[2] );
1259 XYZ aVec12( p2 - p1 );
1260 XYZ aVec13( p3 - p1 );
1261 area += aVec12.Crossed( aVec13 ).Magnitude();
1263 if (myVolume->IsPoly())
1265 for ( int i = 3; i < myCurFace.myNbNodes; ++i )
1267 XYZ pI ( myCurFace.myNodes[i] );
1268 XYZ aVecI( pI - p1 );
1269 area += aVec13.Crossed( aVecI ).Magnitude();
1275 if ( myCurFace.myNbNodes == 4 ) {
1276 XYZ p4 ( myCurFace.myNodes[3] );
1277 XYZ aVec14( p4 - p1 );
1278 area += aVec14.Crossed( aVec13 ).Magnitude();
1284 //================================================================================
1286 * \brief Return index of the node located at face center of a quadratic element like HEX27
1288 //================================================================================
1290 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
1292 if ( myAllFacesNbNodes && myVolumeNodes.size() == 27 ) // classic element with 27 nodes
1294 switch ( faceIndex ) {
1298 return faceIndex + 19;
1304 //=======================================================================
1305 //function : GetOppFaceIndex
1306 //purpose : Return index of the opposite face if it exists, else -1.
1307 //=======================================================================
1309 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1313 MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1317 const int nbHoriFaces = 2;
1319 if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1320 switch ( myVolumeNodes.size() ) {
1323 if ( faceIndex == 0 || faceIndex == 1 )
1324 ind = 1 - faceIndex;
1328 if ( faceIndex <= 1 ) // top or bottom
1329 ind = 1 - faceIndex;
1331 const int nbSideFaces = myAllFacesNbNodes[0];
1332 ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1337 ind = GetOppFaceIndexOfHex( faceIndex );
1345 //=======================================================================
1346 //function : GetOppFaceIndexOfHex
1347 //purpose : Return index of the opposite face of the hexahedron
1348 //=======================================================================
1350 int SMDS_VolumeTool::GetOppFaceIndexOfHex( int faceIndex )
1352 return Hexa_oppF[ faceIndex ];
1355 //=======================================================================
1356 //function : IsLinked
1357 //purpose : return true if theNode1 is linked with theNode2
1358 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1359 //=======================================================================
1361 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1362 const SMDS_MeshNode* theNode2,
1363 const bool theIgnoreMediumNodes) const
1368 if (myVolume->IsPoly()) {
1370 MESSAGE("Warning: bad volumic element");
1373 if ( !myAllFacesNbNodes ) {
1374 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1375 me->myPolyQuantities = myPolyedre->GetQuantities();
1376 myAllFacesNbNodes = &myPolyQuantities[0];
1378 int from, to = 0, d1 = 1, d2 = 2;
1379 if ( myPolyedre->IsQuadratic() ) {
1380 if ( theIgnoreMediumNodes ) {
1386 std::vector<const SMDS_MeshNode*>::const_iterator i;
1387 for (int iface = 0; iface < myNbFaces; iface++)
1390 to += myPolyQuantities[iface];
1391 i = std::find( myVolumeNodes.begin() + from, myVolumeNodes.begin() + to, theNode1 );
1392 if ( i != myVolumeNodes.end() )
1394 if (( theNode2 == *( i-d1 ) ||
1395 theNode2 == *( i+d1 )))
1398 (( theNode2 == *( i-d2 ) ||
1399 theNode2 == *( i+d2 ))))
1406 // find nodes indices
1407 int i1 = -1, i2 = -1, nbFound = 0;
1408 for ( size_t i = 0; i < myVolumeNodes.size() && nbFound < 2; i++ )
1410 if ( myVolumeNodes[ i ] == theNode1 )
1412 else if ( myVolumeNodes[ i ] == theNode2 )
1415 return IsLinked( i1, i2 );
1418 //=======================================================================
1419 //function : IsLinked
1420 //purpose : return true if the node with theNode1Index is linked
1421 // with the node with theNode2Index
1422 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1423 //=======================================================================
1425 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1426 const int theNode2Index,
1427 bool theIgnoreMediumNodes) const
1429 if ( myVolume->IsPoly() ) {
1430 return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1433 int minInd = std::min( theNode1Index, theNode2Index );
1434 int maxInd = std::max( theNode1Index, theNode2Index );
1436 if ( minInd < 0 || maxInd > (int)myVolumeNodes.size() - 1 || maxInd == minInd )
1439 VolumeType type = GetVolumeType();
1440 if ( myVolume->IsQuadratic() )
1442 int firstMediumInd = myVolume->NbCornerNodes();
1443 if ( minInd >= firstMediumInd )
1444 return false; // both nodes are medium - not linked
1445 if ( maxInd < firstMediumInd ) // both nodes are corners
1447 if ( theIgnoreMediumNodes )
1448 type = quadToLinear(type); // to check linkage of corner nodes only
1450 return false; // corner nodes are not linked directly in a quadratic cell
1458 switch ( maxInd - minInd ) {
1459 case 1: return minInd != 3;
1460 case 3: return minInd == 0 || minInd == 4;
1461 case 4: return true;
1468 switch ( maxInd - minInd ) {
1470 case 3: return true;
1475 switch ( maxInd - minInd ) {
1476 case 1: return minInd != 2;
1477 case 2: return minInd == 0 || minInd == 3;
1478 case 3: return true;
1485 case 0: if( maxInd==4 || maxInd==6 || maxInd==7 ) return true;
1486 case 1: if( maxInd==4 || maxInd==5 || maxInd==8 ) return true;
1487 case 2: if( maxInd==5 || maxInd==6 || maxInd==9 ) return true;
1488 case 3: if( maxInd==7 || maxInd==8 || maxInd==9 ) return true;
1496 case 0: if( maxInd==8 || maxInd==11 || maxInd==16 ) return true;
1497 case 1: if( maxInd==8 || maxInd==9 || maxInd==17 ) return true;
1498 case 2: if( maxInd==9 || maxInd==10 || maxInd==18 ) return true;
1499 case 3: if( maxInd==10 || maxInd==11 || maxInd==19 ) return true;
1500 case 4: if( maxInd==12 || maxInd==15 || maxInd==16 ) return true;
1501 case 5: if( maxInd==12 || maxInd==13 || maxInd==17 ) return true;
1502 case 6: if( maxInd==13 || maxInd==14 || maxInd==18 ) return true;
1503 case 7: if( maxInd==14 || maxInd==15 || maxInd==19 ) return true;
1511 case 0: if( maxInd==5 || maxInd==8 || maxInd==9 ) return true;
1512 case 1: if( maxInd==5 || maxInd==6 || maxInd==10 ) return true;
1513 case 2: if( maxInd==6 || maxInd==7 || maxInd==11 ) return true;
1514 case 3: if( maxInd==7 || maxInd==8 || maxInd==12 ) return true;
1515 case 4: if( maxInd==9 || maxInd==10 || maxInd==11 || maxInd==12 ) return true;
1523 case 0: if( maxInd==6 || maxInd==8 || maxInd==12 ) return true;
1524 case 1: if( maxInd==6 || maxInd==7 || maxInd==13 ) return true;
1525 case 2: if( maxInd==7 || maxInd==8 || maxInd==14 ) return true;
1526 case 3: if( maxInd==9 || maxInd==11 || maxInd==12 ) return true;
1527 case 4: if( maxInd==9 || maxInd==10 || maxInd==13 ) return true;
1528 case 5: if( maxInd==10 || maxInd==11 || maxInd==14 ) return true;
1535 const int diff = maxInd-minInd;
1536 if ( diff > 6 ) return false;// not linked top and bottom
1537 if ( diff == 6 ) return true; // linked top and bottom
1538 return diff == 1 || diff == 7;
1545 //=======================================================================
1546 //function : GetNodeIndex
1547 //purpose : Return an index of theNode
1548 //=======================================================================
1550 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1553 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
1554 if ( myVolumeNodes[ i ] == theNode )
1561 //================================================================================
1563 * \brief Fill vector with boundary faces existing in the mesh
1564 * \param faces - vector of found nodes
1565 * \retval int - nb of found faces
1567 //================================================================================
1569 int SMDS_VolumeTool::GetAllExistingFaces(std::vector<const SMDS_MeshElement*> & faces) const
1572 SaveFacet savedFacet( myCurFace );
1574 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1576 if ( const SMDS_MeshElement* face = SMDS_Mesh::FindFace( myCurFace.myNodes ))
1577 faces.push_back( face );
1580 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1581 const SMDS_MeshFace* face = 0;
1582 const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1583 switch ( NbFaceNodes( iF )) {
1585 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1587 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1589 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1590 nodes[3], nodes[4], nodes[5]); break;
1592 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1593 nodes[4], nodes[5], nodes[6], nodes[7]); break;
1596 faces.push_back( face );
1598 return faces.size();
1602 //================================================================================
1604 * \brief Fill vector with boundary edges existing in the mesh
1605 * \param edges - vector of found edges
1606 * \retval int - nb of found faces
1608 //================================================================================
1610 int SMDS_VolumeTool::GetAllExistingEdges(std::vector<const SMDS_MeshElement*> & edges) const
1613 edges.reserve( myVolumeNodes.size() * 2 );
1614 for ( size_t i = 0; i < myVolumeNodes.size()-1; ++i ) {
1615 for ( size_t j = i + 1; j < myVolumeNodes.size(); ++j ) {
1616 if ( IsLinked( i, j )) {
1617 const SMDS_MeshElement* edge =
1618 SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1620 edges.push_back( edge );
1624 return edges.size();
1627 //================================================================================
1629 * \brief Return minimal square distance between connected corner nodes
1631 //================================================================================
1633 double SMDS_VolumeTool::MinLinearSize2() const
1635 double minSize = 1e+100;
1636 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1638 SaveFacet savedFacet( myCurFace );
1640 // it seems that compute distance twice is faster than organization of a sole computing
1641 myCurFace.myIndex = -1;
1642 for ( int iF = 0; iF < myNbFaces; ++iF )
1645 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1647 XYZ n1( myCurFace.myNodes[ iN ]);
1648 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1649 minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
1656 //================================================================================
1658 * \brief Return maximal square distance between connected corner nodes
1660 //================================================================================
1662 double SMDS_VolumeTool::MaxLinearSize2() const
1664 double maxSize = -1e+100;
1665 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1667 SaveFacet savedFacet( myCurFace );
1669 // it seems that compute distance twice is faster than organization of a sole computing
1670 myCurFace.myIndex = -1;
1671 for ( int iF = 0; iF < myNbFaces; ++iF )
1674 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1676 XYZ n1( myCurFace.myNodes[ iN ]);
1677 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1678 maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude());
1685 //================================================================================
1687 * \brief Fast quickly check that only one volume is built on the face nodes
1688 * This check is valid for conformal meshes only
1690 //================================================================================
1692 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1694 const bool isFree = true;
1696 if ( !setFace( faceIndex ))
1699 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1701 const int di = myVolume->IsQuadratic() ? 2 : 1;
1702 const int nbN = ( myCurFace.myNbNodes/di <= 4 && !IsPoly()) ? 3 : myCurFace.myNbNodes/di; // nb nodes to check
1704 SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
1705 while ( eIt->more() )
1707 const SMDS_MeshElement* vol = eIt->next();
1708 if ( vol == myVolume )
1711 for ( iN = 1; iN < nbN; ++iN )
1712 if ( vol->GetNodeIndex( nodes[ iN*di ]) < 0 )
1714 if ( iN == nbN ) // nbN nodes are shared with vol
1716 // if ( vol->IsPoly() || vol->NbFaces() > 6 ) // vol is polyhed or hex prism
1718 // int nb = myCurFace.myNbNodes;
1719 // if ( myVolume->GetEntityType() != vol->GetEntityType() )
1720 // nb -= ( GetCenterNodeIndex(0) > 0 );
1721 // std::set<const SMDS_MeshNode*> faceNodes( nodes, nodes + nb );
1722 // if ( SMDS_VolumeTool( vol ).GetFaceIndex( faceNodes ) < 0 )
1725 if ( otherVol ) *otherVol = vol;
1729 if ( otherVol ) *otherVol = 0;
1733 //================================================================================
1735 * \brief Thorough check that only one volume is built on the face nodes
1737 //================================================================================
1739 bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1741 const bool isFree = true;
1743 if (!setFace( faceIndex ))
1746 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1747 const int nbFaceNodes = myCurFace.myNbNodes;
1749 // evaluate nb of face nodes shared by other volumes
1750 int maxNbShared = -1;
1751 typedef std::map< const SMDS_MeshElement*, int > TElemIntMap;
1752 TElemIntMap volNbShared;
1753 TElemIntMap::iterator vNbIt;
1754 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1755 const SMDS_MeshNode* n = nodes[ iNode ];
1756 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1757 while ( eIt->more() ) {
1758 const SMDS_MeshElement* elem = eIt->next();
1759 if ( elem != myVolume ) {
1760 vNbIt = volNbShared.insert( std::make_pair( elem, 0 )).first;
1762 if ( vNbIt->second > maxNbShared )
1763 maxNbShared = vNbIt->second;
1767 if ( maxNbShared < 3 )
1768 return isFree; // is free
1770 // find volumes laying on the opposite side of the face
1771 // and sharing all nodes
1772 XYZ intNormal; // internal normal
1773 GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
1774 if ( IsFaceExternal( faceIndex ))
1775 intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
1776 XYZ p0 ( nodes[0] ), baryCenter;
1777 for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); ) {
1778 const int& nbShared = (*vNbIt).second;
1779 if ( nbShared >= 3 ) {
1780 SMDS_VolumeTool volume( (*vNbIt).first );
1781 volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
1782 XYZ intNormal2( baryCenter - p0 );
1783 if ( intNormal.Dot( intNormal2 ) < 0 ) {
1785 if ( nbShared >= nbFaceNodes )
1787 // a volume shares the whole facet
1788 if ( otherVol ) *otherVol = vNbIt->first;
1795 // remove a volume from volNbShared map
1796 volNbShared.erase( vNbIt++ );
1799 // here volNbShared contains only volumes laying on the opposite side of
1800 // the face and sharing 3 or more but not all face nodes with myVolume
1801 if ( volNbShared.size() < 2 ) {
1802 return isFree; // is free
1805 // check if the whole area of a face is shared
1806 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1808 const SMDS_MeshNode* n = nodes[ iNode ];
1809 // check if n is shared by one of volumes of volNbShared
1810 bool isShared = false;
1811 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1812 while ( eIt->more() && !isShared )
1813 isShared = volNbShared.count( eIt->next() );
1817 if ( otherVol ) *otherVol = volNbShared.begin()->first;
1820 // if ( !myVolume->IsPoly() )
1822 // bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
1823 // for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1824 // SMDS_VolumeTool volume( (*vNbIt).first );
1825 // bool prevLinkShared = false;
1826 // int nbSharedLinks = 0;
1827 // for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1828 // bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
1829 // if ( linkShared )
1831 // if ( linkShared && prevLinkShared &&
1832 // volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
1833 // isShared[ iNode ] = true;
1834 // prevLinkShared = linkShared;
1836 // if ( nbSharedLinks == nbFaceNodes )
1837 // return !free; // is not free
1838 // if ( nbFaceNodes == 4 ) {
1839 // // check traingle parts 1 & 3
1840 // if ( isShared[1] && isShared[3] )
1841 // return !free; // is not free
1842 // // check triangle parts 0 & 2;
1843 // // 0 part could not be checked in the loop; check it here
1844 // if ( isShared[2] && prevLinkShared &&
1845 // volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
1846 // volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
1847 // return !free; // is not free
1854 //=======================================================================
1855 //function : GetFaceIndex
1856 //purpose : Return index of a face formed by theFaceNodes
1857 //=======================================================================
1859 int SMDS_VolumeTool::GetFaceIndex( const std::set<const SMDS_MeshNode*>& theFaceNodes,
1860 const int theFaceIndexHint ) const
1862 if ( theFaceIndexHint >= 0 )
1864 int nbNodes = NbFaceNodes( theFaceIndexHint );
1865 if ( nbNodes == (int) theFaceNodes.size() )
1867 const SMDS_MeshNode** nodes = GetFaceNodes( theFaceIndexHint );
1869 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1874 return theFaceIndexHint;
1877 for ( int iFace = 0; iFace < myNbFaces; iFace++ )
1879 if ( iFace == theFaceIndexHint )
1881 int nbNodes = NbFaceNodes( iFace );
1882 if ( nbNodes == (int) theFaceNodes.size() )
1884 const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
1886 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1897 //=======================================================================
1898 //function : GetFaceIndex
1899 //purpose : Return index of a face formed by theFaceNodes
1900 //=======================================================================
1902 /*int SMDS_VolumeTool::GetFaceIndex( const std::set<int>& theFaceNodesIndices )
1904 for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1905 const int* nodes = GetFaceNodesIndices( iFace );
1906 int nbFaceNodes = NbFaceNodes( iFace );
1907 std::set<int> nodeSet;
1908 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1909 nodeSet.insert( nodes[ iNode ] );
1910 if ( theFaceNodesIndices == nodeSet )
1916 //=======================================================================
1917 //function : setFace
1919 //=======================================================================
1921 bool SMDS_VolumeTool::setFace( int faceIndex ) const
1926 if ( myCurFace.myIndex == faceIndex )
1929 myCurFace.myIndex = -1;
1931 if ( faceIndex < 0 || faceIndex >= NbFaces() )
1934 if (myVolume->IsPoly())
1936 if ( !myPolyedre ) {
1937 MESSAGE("Warning: bad volumic element");
1942 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1943 if ( !myAllFacesNbNodes ) {
1944 me->myPolyQuantities = myPolyedre->GetQuantities();
1945 myAllFacesNbNodes = &myPolyQuantities[0];
1947 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
1948 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
1949 me->myPolyIndices.resize( myCurFace.myNbNodes + 1 );
1950 myCurFace.myNodeIndices = & me->myPolyIndices[0];
1951 int shift = std::accumulate( myAllFacesNbNodes, myAllFacesNbNodes+faceIndex, 0 );
1952 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
1954 myCurFace.myNodes [ iNode ] = myVolumeNodes[ shift + iNode ];
1955 myCurFace.myNodeIndices[ iNode ] = shift + iNode;
1957 myCurFace.myNodes [ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ]; // last = first
1958 myCurFace.myNodeIndices[ myCurFace.myNbNodes ] = myCurFace.myNodeIndices[ 0 ];
1960 // check orientation
1961 if (myExternalFaces)
1963 myCurFace.myIndex = faceIndex; // avoid infinite recursion in IsFaceExternal()
1964 myExternalFaces = false; // force normal computation by IsFaceExternal()
1965 if ( !IsFaceExternal( faceIndex ))
1966 std::reverse( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
1967 myExternalFaces = true;
1972 if ( !myAllFacesNodeIndices_F )
1974 // choose data for an element type
1975 switch ( myVolumeNodes.size() ) {
1977 myAllFacesNodeIndices_F = &Tetra_F [0][0];
1978 //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
1979 myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
1980 myAllFacesNbNodes = Tetra_nbN;
1981 myMaxFaceNbNodes = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
1984 myAllFacesNodeIndices_F = &Pyramid_F [0][0];
1985 //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
1986 myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
1987 myAllFacesNbNodes = Pyramid_nbN;
1988 myMaxFaceNbNodes = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
1991 myAllFacesNodeIndices_F = &Penta_F [0][0];
1992 //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
1993 myAllFacesNodeIndices_RE = &Penta_RE[0][0];
1994 myAllFacesNbNodes = Penta_nbN;
1995 myMaxFaceNbNodes = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
1998 myAllFacesNodeIndices_F = &Hexa_F [0][0];
1999 ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
2000 myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
2001 myAllFacesNbNodes = Hexa_nbN;
2002 myMaxFaceNbNodes = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
2005 myAllFacesNodeIndices_F = &QuadTetra_F [0][0];
2006 //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
2007 myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
2008 myAllFacesNbNodes = QuadTetra_nbN;
2009 myMaxFaceNbNodes = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
2012 myAllFacesNodeIndices_F = &QuadPyram_F [0][0];
2013 //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
2014 myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
2015 myAllFacesNbNodes = QuadPyram_nbN;
2016 myMaxFaceNbNodes = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
2019 myAllFacesNodeIndices_F = &QuadPenta_F [0][0];
2020 //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
2021 myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
2022 myAllFacesNbNodes = QuadPenta_nbN;
2023 myMaxFaceNbNodes = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
2027 myAllFacesNodeIndices_F = &QuadHexa_F [0][0];
2028 //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
2029 myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
2030 myAllFacesNbNodes = QuadHexa_nbN;
2031 myMaxFaceNbNodes = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
2032 if ( !myIgnoreCentralNodes && myVolumeNodes.size() == 27 )
2034 myAllFacesNodeIndices_F = &TriQuadHexa_F [0][0];
2035 //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
2036 myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
2037 myAllFacesNbNodes = TriQuadHexa_nbN;
2038 myMaxFaceNbNodes = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
2042 myAllFacesNodeIndices_F = &HexPrism_F [0][0];
2043 //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
2044 myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
2045 myAllFacesNbNodes = HexPrism_nbN;
2046 myMaxFaceNbNodes = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
2052 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2053 // if ( myExternalFaces )
2054 // myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2056 // myCurFace.myNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
2057 myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2060 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2061 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2062 myCurFace.myNodes[ iNode ] = myVolumeNodes[ myCurFace.myNodeIndices[ iNode ]];
2063 myCurFace.myNodes[ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ];
2066 myCurFace.myIndex = faceIndex;
2071 //=======================================================================
2072 //function : GetType
2073 //purpose : return VolumeType by nb of nodes in a volume
2074 //=======================================================================
2076 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
2078 switch ( nbNodes ) {
2079 case 4: return TETRA;
2080 case 5: return PYRAM;
2081 case 6: return PENTA;
2082 case 8: return HEXA;
2083 case 10: return QUAD_TETRA;
2084 case 13: return QUAD_PYRAM;
2085 case 15: return QUAD_PENTA;
2087 case 27: return QUAD_HEXA;
2088 case 12: return HEX_PRISM;
2089 default:return UNKNOWN;
2093 //=======================================================================
2094 //function : NbFaces
2095 //purpose : return nb of faces by volume type
2096 //=======================================================================
2098 int SMDS_VolumeTool::NbFaces( VolumeType type )
2102 case QUAD_TETRA: return 4;
2104 case QUAD_PYRAM: return 5;
2106 case QUAD_PENTA: return 5;
2108 case QUAD_HEXA : return 6;
2109 case HEX_PRISM : return 8;
2114 //================================================================================
2116 * \brief Useful to know nb of corner nodes of a quadratic volume
2117 * \param type - volume type
2118 * \retval int - nb of corner nodes
2120 //================================================================================
2122 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
2126 case QUAD_TETRA: return 4;
2128 case QUAD_PYRAM: return 5;
2130 case QUAD_PENTA: return 6;
2132 case QUAD_HEXA : return 8;
2133 case HEX_PRISM : return 12;
2140 //=======================================================================
2141 //function : GetFaceNodesIndices
2142 //purpose : Return the array of face nodes indices
2143 // To comfort link iteration, the array
2144 // length == NbFaceNodes( faceIndex ) + 1 and
2145 // the last node index == the first one.
2146 //=======================================================================
2148 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
2153 case TETRA: return Tetra_F[ faceIndex ];
2154 case PYRAM: return Pyramid_F[ faceIndex ];
2155 case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
2156 case HEXA: return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
2157 case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
2158 case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
2159 case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
2160 // what about SMDSEntity_TriQuad_Hexa?
2161 case QUAD_HEXA: return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
2162 case HEX_PRISM: return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
2168 //=======================================================================
2169 //function : NbFaceNodes
2170 //purpose : Return number of nodes in the array of face nodes
2171 //=======================================================================
2173 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
2177 case TETRA: return Tetra_nbN[ faceIndex ];
2178 case PYRAM: return Pyramid_nbN[ faceIndex ];
2179 case PENTA: return Penta_nbN[ faceIndex ];
2180 case HEXA: return Hexa_nbN[ faceIndex ];
2181 case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
2182 case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
2183 case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
2184 // what about SMDSEntity_TriQuad_Hexa?
2185 case QUAD_HEXA: return QuadHexa_nbN[ faceIndex ];
2186 case HEX_PRISM: return HexPrism_nbN[ faceIndex ];
2192 //=======================================================================
2193 //function : Element
2194 //purpose : return element
2195 //=======================================================================
2197 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
2199 return static_cast<const SMDS_MeshVolume*>( myVolume );
2202 //=======================================================================
2204 //purpose : return element ID
2205 //=======================================================================
2207 int SMDS_VolumeTool::ID() const
2209 return myVolume ? myVolume->GetID() : 0;