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 for ( int i = 3*iQuad; i < myCurFace.myNbNodes; i += iQuad )
1208 XYZ p4 ( myCurFace.myNodes[i] );
1209 XYZ aVec14( p4 - p1 );
1210 XYZ cross2 = aVec13.Crossed( aVec14 );
1211 cross = cross + cross2;
1215 double size = cross.Magnitude();
1216 if ( size <= std::numeric_limits<double>::min() )
1226 //================================================================================
1228 * \brief Return barycenter of a face
1230 //================================================================================
1232 bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const
1234 if ( !setFace( faceIndex ))
1238 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1240 X += myCurFace.myNodes[i]->X() / myCurFace.myNbNodes;
1241 Y += myCurFace.myNodes[i]->Y() / myCurFace.myNbNodes;
1242 Z += myCurFace.myNodes[i]->Z() / myCurFace.myNbNodes;
1247 //=======================================================================
1248 //function : GetFaceArea
1249 //purpose : Return face area
1250 //=======================================================================
1252 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
1255 if ( !setFace( faceIndex ))
1258 XYZ p1 ( myCurFace.myNodes[0] );
1259 XYZ p2 ( myCurFace.myNodes[1] );
1260 XYZ p3 ( myCurFace.myNodes[2] );
1261 XYZ aVec12( p2 - p1 );
1262 XYZ aVec13( p3 - p1 );
1263 area += aVec12.Crossed( aVec13 ).Magnitude();
1265 if (myVolume->IsPoly())
1267 for ( int i = 3; i < myCurFace.myNbNodes; ++i )
1269 XYZ pI ( myCurFace.myNodes[i] );
1270 XYZ aVecI( pI - p1 );
1271 area += aVec13.Crossed( aVecI ).Magnitude();
1277 if ( myCurFace.myNbNodes == 4 ) {
1278 XYZ p4 ( myCurFace.myNodes[3] );
1279 XYZ aVec14( p4 - p1 );
1280 area += aVec14.Crossed( aVec13 ).Magnitude();
1286 //================================================================================
1288 * \brief Return index of the node located at face center of a quadratic element like HEX27
1290 //================================================================================
1292 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
1294 if ( myAllFacesNbNodes && myVolumeNodes.size() == 27 ) // classic element with 27 nodes
1296 switch ( faceIndex ) {
1300 return faceIndex + 19;
1306 //=======================================================================
1307 //function : GetOppFaceIndex
1308 //purpose : Return index of the opposite face if it exists, else -1.
1309 //=======================================================================
1311 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1315 MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1319 const int nbHoriFaces = 2;
1321 if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1322 switch ( myVolumeNodes.size() ) {
1325 if ( faceIndex == 0 || faceIndex == 1 )
1326 ind = 1 - faceIndex;
1330 if ( faceIndex <= 1 ) // top or bottom
1331 ind = 1 - faceIndex;
1333 const int nbSideFaces = myAllFacesNbNodes[0];
1334 ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1339 ind = GetOppFaceIndexOfHex( faceIndex );
1347 //=======================================================================
1348 //function : GetOppFaceIndexOfHex
1349 //purpose : Return index of the opposite face of the hexahedron
1350 //=======================================================================
1352 int SMDS_VolumeTool::GetOppFaceIndexOfHex( int faceIndex )
1354 return Hexa_oppF[ faceIndex ];
1357 //=======================================================================
1358 //function : IsLinked
1359 //purpose : return true if theNode1 is linked with theNode2
1360 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1361 //=======================================================================
1363 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1364 const SMDS_MeshNode* theNode2,
1365 const bool theIgnoreMediumNodes) const
1370 if (myVolume->IsPoly()) {
1372 MESSAGE("Warning: bad volumic element");
1375 if ( !myAllFacesNbNodes ) {
1376 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1377 me->myPolyQuantities = myPolyedre->GetQuantities();
1378 myAllFacesNbNodes = &myPolyQuantities[0];
1380 int from, to = 0, d1 = 1, d2 = 2;
1381 if ( myPolyedre->IsQuadratic() ) {
1382 if ( theIgnoreMediumNodes ) {
1388 std::vector<const SMDS_MeshNode*>::const_iterator i;
1389 for (int iface = 0; iface < myNbFaces; iface++)
1392 to += myPolyQuantities[iface];
1393 i = std::find( myVolumeNodes.begin() + from, myVolumeNodes.begin() + to, theNode1 );
1394 if ( i != myVolumeNodes.end() )
1396 if (( theNode2 == *( i-d1 ) ||
1397 theNode2 == *( i+d1 )))
1400 (( theNode2 == *( i-d2 ) ||
1401 theNode2 == *( i+d2 ))))
1408 // find nodes indices
1409 int i1 = -1, i2 = -1, nbFound = 0;
1410 for ( size_t i = 0; i < myVolumeNodes.size() && nbFound < 2; i++ )
1412 if ( myVolumeNodes[ i ] == theNode1 )
1414 else if ( myVolumeNodes[ i ] == theNode2 )
1417 return IsLinked( i1, i2 );
1420 //=======================================================================
1421 //function : IsLinked
1422 //purpose : return true if the node with theNode1Index is linked
1423 // with the node with theNode2Index
1424 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1425 //=======================================================================
1427 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1428 const int theNode2Index,
1429 bool theIgnoreMediumNodes) const
1431 if ( myVolume->IsPoly() ) {
1432 return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1435 int minInd = std::min( theNode1Index, theNode2Index );
1436 int maxInd = std::max( theNode1Index, theNode2Index );
1438 if ( minInd < 0 || maxInd > (int)myVolumeNodes.size() - 1 || maxInd == minInd )
1441 VolumeType type = GetVolumeType();
1442 if ( myVolume->IsQuadratic() )
1444 int firstMediumInd = myVolume->NbCornerNodes();
1445 if ( minInd >= firstMediumInd )
1446 return false; // both nodes are medium - not linked
1447 if ( maxInd < firstMediumInd ) // both nodes are corners
1449 if ( theIgnoreMediumNodes )
1450 type = quadToLinear(type); // to check linkage of corner nodes only
1452 return false; // corner nodes are not linked directly in a quadratic cell
1460 switch ( maxInd - minInd ) {
1461 case 1: return minInd != 3;
1462 case 3: return minInd == 0 || minInd == 4;
1463 case 4: return true;
1470 switch ( maxInd - minInd ) {
1472 case 3: return true;
1477 switch ( maxInd - minInd ) {
1478 case 1: return minInd != 2;
1479 case 2: return minInd == 0 || minInd == 3;
1480 case 3: return true;
1487 case 0: if( maxInd==4 || maxInd==6 || maxInd==7 ) return true;
1488 case 1: if( maxInd==4 || maxInd==5 || maxInd==8 ) return true;
1489 case 2: if( maxInd==5 || maxInd==6 || maxInd==9 ) return true;
1490 case 3: if( maxInd==7 || maxInd==8 || maxInd==9 ) return true;
1498 case 0: if( maxInd==8 || maxInd==11 || maxInd==16 ) return true;
1499 case 1: if( maxInd==8 || maxInd==9 || maxInd==17 ) return true;
1500 case 2: if( maxInd==9 || maxInd==10 || maxInd==18 ) return true;
1501 case 3: if( maxInd==10 || maxInd==11 || maxInd==19 ) return true;
1502 case 4: if( maxInd==12 || maxInd==15 || maxInd==16 ) return true;
1503 case 5: if( maxInd==12 || maxInd==13 || maxInd==17 ) return true;
1504 case 6: if( maxInd==13 || maxInd==14 || maxInd==18 ) return true;
1505 case 7: if( maxInd==14 || maxInd==15 || maxInd==19 ) return true;
1513 case 0: if( maxInd==5 || maxInd==8 || maxInd==9 ) return true;
1514 case 1: if( maxInd==5 || maxInd==6 || maxInd==10 ) return true;
1515 case 2: if( maxInd==6 || maxInd==7 || maxInd==11 ) return true;
1516 case 3: if( maxInd==7 || maxInd==8 || maxInd==12 ) return true;
1517 case 4: if( maxInd==9 || maxInd==10 || maxInd==11 || maxInd==12 ) return true;
1525 case 0: if( maxInd==6 || maxInd==8 || maxInd==12 ) return true;
1526 case 1: if( maxInd==6 || maxInd==7 || maxInd==13 ) return true;
1527 case 2: if( maxInd==7 || maxInd==8 || maxInd==14 ) return true;
1528 case 3: if( maxInd==9 || maxInd==11 || maxInd==12 ) return true;
1529 case 4: if( maxInd==9 || maxInd==10 || maxInd==13 ) return true;
1530 case 5: if( maxInd==10 || maxInd==11 || maxInd==14 ) return true;
1537 const int diff = maxInd-minInd;
1538 if ( diff > 6 ) return false;// not linked top and bottom
1539 if ( diff == 6 ) return true; // linked top and bottom
1540 return diff == 1 || diff == 7;
1547 //=======================================================================
1548 //function : GetNodeIndex
1549 //purpose : Return an index of theNode
1550 //=======================================================================
1552 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1555 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
1556 if ( myVolumeNodes[ i ] == theNode )
1563 //================================================================================
1565 * \brief Fill vector with boundary faces existing in the mesh
1566 * \param faces - vector of found nodes
1567 * \retval int - nb of found faces
1569 //================================================================================
1571 int SMDS_VolumeTool::GetAllExistingFaces(std::vector<const SMDS_MeshElement*> & faces) const
1574 SaveFacet savedFacet( myCurFace );
1576 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1578 if ( const SMDS_MeshElement* face = SMDS_Mesh::FindFace( myCurFace.myNodes ))
1579 faces.push_back( face );
1582 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1583 const SMDS_MeshFace* face = 0;
1584 const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1585 switch ( NbFaceNodes( iF )) {
1587 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1589 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1591 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1592 nodes[3], nodes[4], nodes[5]); break;
1594 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1595 nodes[4], nodes[5], nodes[6], nodes[7]); break;
1598 faces.push_back( face );
1600 return faces.size();
1604 //================================================================================
1606 * \brief Fill vector with boundary edges existing in the mesh
1607 * \param edges - vector of found edges
1608 * \retval int - nb of found faces
1610 //================================================================================
1612 int SMDS_VolumeTool::GetAllExistingEdges(std::vector<const SMDS_MeshElement*> & edges) const
1615 edges.reserve( myVolumeNodes.size() * 2 );
1616 for ( size_t i = 0; i < myVolumeNodes.size()-1; ++i ) {
1617 for ( size_t j = i + 1; j < myVolumeNodes.size(); ++j ) {
1618 if ( IsLinked( i, j )) {
1619 const SMDS_MeshElement* edge =
1620 SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1622 edges.push_back( edge );
1626 return edges.size();
1629 //================================================================================
1631 * \brief Return minimal square distance between connected corner nodes
1633 //================================================================================
1635 double SMDS_VolumeTool::MinLinearSize2() const
1637 double minSize = 1e+100;
1638 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1640 SaveFacet savedFacet( myCurFace );
1642 // it seems that compute distance twice is faster than organization of a sole computing
1643 myCurFace.myIndex = -1;
1644 for ( int iF = 0; iF < myNbFaces; ++iF )
1647 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1649 XYZ n1( myCurFace.myNodes[ iN ]);
1650 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1651 minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
1658 //================================================================================
1660 * \brief Return maximal square distance between connected corner nodes
1662 //================================================================================
1664 double SMDS_VolumeTool::MaxLinearSize2() const
1666 double maxSize = -1e+100;
1667 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1669 SaveFacet savedFacet( myCurFace );
1671 // it seems that compute distance twice is faster than organization of a sole computing
1672 myCurFace.myIndex = -1;
1673 for ( int iF = 0; iF < myNbFaces; ++iF )
1676 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
1678 XYZ n1( myCurFace.myNodes[ iN ]);
1679 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
1680 maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude());
1687 //================================================================================
1689 * \brief Fast quickly check that only one volume is built on the face nodes
1690 * This check is valid for conformal meshes only
1692 //================================================================================
1694 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1696 const bool isFree = true;
1698 if ( !setFace( faceIndex ))
1701 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1703 const int di = myVolume->IsQuadratic() ? 2 : 1;
1704 const int nbN = ( myCurFace.myNbNodes/di <= 4 && !IsPoly()) ? 3 : myCurFace.myNbNodes/di; // nb nodes to check
1706 SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
1707 while ( eIt->more() )
1709 const SMDS_MeshElement* vol = eIt->next();
1710 if ( vol == myVolume )
1713 for ( iN = 1; iN < nbN; ++iN )
1714 if ( vol->GetNodeIndex( nodes[ iN*di ]) < 0 )
1716 if ( iN == nbN ) // nbN nodes are shared with vol
1718 // if ( vol->IsPoly() || vol->NbFaces() > 6 ) // vol is polyhed or hex prism
1720 // int nb = myCurFace.myNbNodes;
1721 // if ( myVolume->GetEntityType() != vol->GetEntityType() )
1722 // nb -= ( GetCenterNodeIndex(0) > 0 );
1723 // std::set<const SMDS_MeshNode*> faceNodes( nodes, nodes + nb );
1724 // if ( SMDS_VolumeTool( vol ).GetFaceIndex( faceNodes ) < 0 )
1727 if ( otherVol ) *otherVol = vol;
1731 if ( otherVol ) *otherVol = 0;
1735 //================================================================================
1737 * \brief Thorough check that only one volume is built on the face nodes
1739 //================================================================================
1741 bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1743 const bool isFree = true;
1745 if (!setFace( faceIndex ))
1748 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1749 const int nbFaceNodes = myCurFace.myNbNodes;
1751 // evaluate nb of face nodes shared by other volumes
1752 int maxNbShared = -1;
1753 typedef std::map< const SMDS_MeshElement*, int > TElemIntMap;
1754 TElemIntMap volNbShared;
1755 TElemIntMap::iterator vNbIt;
1756 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1757 const SMDS_MeshNode* n = nodes[ iNode ];
1758 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1759 while ( eIt->more() ) {
1760 const SMDS_MeshElement* elem = eIt->next();
1761 if ( elem != myVolume ) {
1762 vNbIt = volNbShared.insert( std::make_pair( elem, 0 )).first;
1764 if ( vNbIt->second > maxNbShared )
1765 maxNbShared = vNbIt->second;
1769 if ( maxNbShared < 3 )
1770 return isFree; // is free
1772 // find volumes laying on the opposite side of the face
1773 // and sharing all nodes
1774 XYZ intNormal; // internal normal
1775 GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
1776 if ( IsFaceExternal( faceIndex ))
1777 intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
1778 XYZ p0 ( nodes[0] ), baryCenter;
1779 for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); ) {
1780 const int& nbShared = (*vNbIt).second;
1781 if ( nbShared >= 3 ) {
1782 SMDS_VolumeTool volume( (*vNbIt).first );
1783 volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
1784 XYZ intNormal2( baryCenter - p0 );
1785 if ( intNormal.Dot( intNormal2 ) < 0 ) {
1787 if ( nbShared >= nbFaceNodes )
1789 // a volume shares the whole facet
1790 if ( otherVol ) *otherVol = vNbIt->first;
1797 // remove a volume from volNbShared map
1798 volNbShared.erase( vNbIt++ );
1801 // here volNbShared contains only volumes laying on the opposite side of
1802 // the face and sharing 3 or more but not all face nodes with myVolume
1803 if ( volNbShared.size() < 2 ) {
1804 return isFree; // is free
1807 // check if the whole area of a face is shared
1808 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1810 const SMDS_MeshNode* n = nodes[ iNode ];
1811 // check if n is shared by one of volumes of volNbShared
1812 bool isShared = false;
1813 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1814 while ( eIt->more() && !isShared )
1815 isShared = volNbShared.count( eIt->next() );
1819 if ( otherVol ) *otherVol = volNbShared.begin()->first;
1822 // if ( !myVolume->IsPoly() )
1824 // bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
1825 // for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1826 // SMDS_VolumeTool volume( (*vNbIt).first );
1827 // bool prevLinkShared = false;
1828 // int nbSharedLinks = 0;
1829 // for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1830 // bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
1831 // if ( linkShared )
1833 // if ( linkShared && prevLinkShared &&
1834 // volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
1835 // isShared[ iNode ] = true;
1836 // prevLinkShared = linkShared;
1838 // if ( nbSharedLinks == nbFaceNodes )
1839 // return !free; // is not free
1840 // if ( nbFaceNodes == 4 ) {
1841 // // check traingle parts 1 & 3
1842 // if ( isShared[1] && isShared[3] )
1843 // return !free; // is not free
1844 // // check triangle parts 0 & 2;
1845 // // 0 part could not be checked in the loop; check it here
1846 // if ( isShared[2] && prevLinkShared &&
1847 // volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
1848 // volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
1849 // return !free; // is not free
1856 //=======================================================================
1857 //function : GetFaceIndex
1858 //purpose : Return index of a face formed by theFaceNodes
1859 //=======================================================================
1861 int SMDS_VolumeTool::GetFaceIndex( const std::set<const SMDS_MeshNode*>& theFaceNodes,
1862 const int theFaceIndexHint ) const
1864 if ( theFaceIndexHint >= 0 )
1866 int nbNodes = NbFaceNodes( theFaceIndexHint );
1867 if ( nbNodes == (int) theFaceNodes.size() )
1869 const SMDS_MeshNode** nodes = GetFaceNodes( theFaceIndexHint );
1871 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1876 return theFaceIndexHint;
1879 for ( int iFace = 0; iFace < myNbFaces; iFace++ )
1881 if ( iFace == theFaceIndexHint )
1883 int nbNodes = NbFaceNodes( iFace );
1884 if ( nbNodes == (int) theFaceNodes.size() )
1886 const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
1888 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1899 //=======================================================================
1900 //function : GetFaceIndex
1901 //purpose : Return index of a face formed by theFaceNodes
1902 //=======================================================================
1904 /*int SMDS_VolumeTool::GetFaceIndex( const std::set<int>& theFaceNodesIndices )
1906 for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1907 const int* nodes = GetFaceNodesIndices( iFace );
1908 int nbFaceNodes = NbFaceNodes( iFace );
1909 std::set<int> nodeSet;
1910 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1911 nodeSet.insert( nodes[ iNode ] );
1912 if ( theFaceNodesIndices == nodeSet )
1918 //=======================================================================
1919 //function : setFace
1921 //=======================================================================
1923 bool SMDS_VolumeTool::setFace( int faceIndex ) const
1928 if ( myCurFace.myIndex == faceIndex )
1931 myCurFace.myIndex = -1;
1933 if ( faceIndex < 0 || faceIndex >= NbFaces() )
1936 if (myVolume->IsPoly())
1938 if ( !myPolyedre ) {
1939 MESSAGE("Warning: bad volumic element");
1944 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1945 if ( !myAllFacesNbNodes ) {
1946 me->myPolyQuantities = myPolyedre->GetQuantities();
1947 myAllFacesNbNodes = &myPolyQuantities[0];
1949 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
1950 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
1951 me->myPolyIndices.resize( myCurFace.myNbNodes + 1 );
1952 myCurFace.myNodeIndices = & me->myPolyIndices[0];
1953 int shift = std::accumulate( myAllFacesNbNodes, myAllFacesNbNodes+faceIndex, 0 );
1954 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
1956 myCurFace.myNodes [ iNode ] = myVolumeNodes[ shift + iNode ];
1957 myCurFace.myNodeIndices[ iNode ] = shift + iNode;
1959 myCurFace.myNodes [ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ]; // last = first
1960 myCurFace.myNodeIndices[ myCurFace.myNbNodes ] = myCurFace.myNodeIndices[ 0 ];
1962 // check orientation
1963 if (myExternalFaces)
1965 myCurFace.myIndex = faceIndex; // avoid infinite recursion in IsFaceExternal()
1966 myExternalFaces = false; // force normal computation by IsFaceExternal()
1967 if ( !IsFaceExternal( faceIndex ))
1968 std::reverse( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
1969 myExternalFaces = true;
1974 if ( !myAllFacesNodeIndices_F )
1976 // choose data for an element type
1977 switch ( myVolumeNodes.size() ) {
1979 myAllFacesNodeIndices_F = &Tetra_F [0][0];
1980 //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
1981 myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
1982 myAllFacesNbNodes = Tetra_nbN;
1983 myMaxFaceNbNodes = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
1986 myAllFacesNodeIndices_F = &Pyramid_F [0][0];
1987 //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
1988 myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
1989 myAllFacesNbNodes = Pyramid_nbN;
1990 myMaxFaceNbNodes = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
1993 myAllFacesNodeIndices_F = &Penta_F [0][0];
1994 //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
1995 myAllFacesNodeIndices_RE = &Penta_RE[0][0];
1996 myAllFacesNbNodes = Penta_nbN;
1997 myMaxFaceNbNodes = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
2000 myAllFacesNodeIndices_F = &Hexa_F [0][0];
2001 ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
2002 myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
2003 myAllFacesNbNodes = Hexa_nbN;
2004 myMaxFaceNbNodes = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
2007 myAllFacesNodeIndices_F = &QuadTetra_F [0][0];
2008 //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
2009 myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
2010 myAllFacesNbNodes = QuadTetra_nbN;
2011 myMaxFaceNbNodes = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
2014 myAllFacesNodeIndices_F = &QuadPyram_F [0][0];
2015 //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
2016 myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
2017 myAllFacesNbNodes = QuadPyram_nbN;
2018 myMaxFaceNbNodes = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
2021 myAllFacesNodeIndices_F = &QuadPenta_F [0][0];
2022 //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
2023 myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
2024 myAllFacesNbNodes = QuadPenta_nbN;
2025 myMaxFaceNbNodes = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
2029 myAllFacesNodeIndices_F = &QuadHexa_F [0][0];
2030 //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
2031 myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
2032 myAllFacesNbNodes = QuadHexa_nbN;
2033 myMaxFaceNbNodes = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
2034 if ( !myIgnoreCentralNodes && myVolumeNodes.size() == 27 )
2036 myAllFacesNodeIndices_F = &TriQuadHexa_F [0][0];
2037 //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
2038 myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
2039 myAllFacesNbNodes = TriQuadHexa_nbN;
2040 myMaxFaceNbNodes = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
2044 myAllFacesNodeIndices_F = &HexPrism_F [0][0];
2045 //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
2046 myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
2047 myAllFacesNbNodes = HexPrism_nbN;
2048 myMaxFaceNbNodes = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
2054 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2055 // if ( myExternalFaces )
2056 // myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2058 // myCurFace.myNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
2059 myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2062 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2063 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2064 myCurFace.myNodes[ iNode ] = myVolumeNodes[ myCurFace.myNodeIndices[ iNode ]];
2065 myCurFace.myNodes[ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ];
2068 myCurFace.myIndex = faceIndex;
2073 //=======================================================================
2074 //function : GetType
2075 //purpose : return VolumeType by nb of nodes in a volume
2076 //=======================================================================
2078 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
2080 switch ( nbNodes ) {
2081 case 4: return TETRA;
2082 case 5: return PYRAM;
2083 case 6: return PENTA;
2084 case 8: return HEXA;
2085 case 10: return QUAD_TETRA;
2086 case 13: return QUAD_PYRAM;
2087 case 15: return QUAD_PENTA;
2089 case 27: return QUAD_HEXA;
2090 case 12: return HEX_PRISM;
2091 default:return UNKNOWN;
2095 //=======================================================================
2096 //function : NbFaces
2097 //purpose : return nb of faces by volume type
2098 //=======================================================================
2100 int SMDS_VolumeTool::NbFaces( VolumeType type )
2104 case QUAD_TETRA: return 4;
2106 case QUAD_PYRAM: return 5;
2108 case QUAD_PENTA: return 5;
2110 case QUAD_HEXA : return 6;
2111 case HEX_PRISM : return 8;
2116 //================================================================================
2118 * \brief Useful to know nb of corner nodes of a quadratic volume
2119 * \param type - volume type
2120 * \retval int - nb of corner nodes
2122 //================================================================================
2124 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
2128 case QUAD_TETRA: return 4;
2130 case QUAD_PYRAM: return 5;
2132 case QUAD_PENTA: return 6;
2134 case QUAD_HEXA : return 8;
2135 case HEX_PRISM : return 12;
2142 //=======================================================================
2143 //function : GetFaceNodesIndices
2144 //purpose : Return the array of face nodes indices
2145 // To comfort link iteration, the array
2146 // length == NbFaceNodes( faceIndex ) + 1 and
2147 // the last node index == the first one.
2148 //=======================================================================
2150 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
2155 case TETRA: return Tetra_F[ faceIndex ];
2156 case PYRAM: return Pyramid_F[ faceIndex ];
2157 case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
2158 case HEXA: return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
2159 case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
2160 case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
2161 case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
2162 // what about SMDSEntity_TriQuad_Hexa?
2163 case QUAD_HEXA: return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
2164 case HEX_PRISM: return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
2170 //=======================================================================
2171 //function : NbFaceNodes
2172 //purpose : Return number of nodes in the array of face nodes
2173 //=======================================================================
2175 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
2179 case TETRA: return Tetra_nbN[ faceIndex ];
2180 case PYRAM: return Pyramid_nbN[ faceIndex ];
2181 case PENTA: return Penta_nbN[ faceIndex ];
2182 case HEXA: return Hexa_nbN[ faceIndex ];
2183 case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
2184 case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
2185 case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
2186 // what about SMDSEntity_TriQuad_Hexa?
2187 case QUAD_HEXA: return QuadHexa_nbN[ faceIndex ];
2188 case HEX_PRISM: return HexPrism_nbN[ faceIndex ];
2194 //=======================================================================
2195 //function : Element
2196 //purpose : return element
2197 //=======================================================================
2199 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
2201 return static_cast<const SMDS_MeshVolume*>( myVolume );
2204 //=======================================================================
2206 //purpose : return element ID
2207 //=======================================================================
2209 int SMDS_VolumeTool::ID() const
2211 return myVolume ? myVolume->GetID() : 0;