1 // Copyright (C) 2007-2024 CEA, EDF, 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>
49 // ======================================================
50 // Node indices in faces depending on volume orientation
51 // making most faces normals external
52 // ======================================================
53 // For all elements, 0-th face is bottom based on the first nodes.
54 // For prismatic elements (tetra,hexa,prisms), 1-th face is a top one.
55 // For all elements, side faces follow order of bottom nodes
56 // ======================================================
64 // N0 +---|---+ N1 TETRAHEDRON
72 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
73 { 0, 1, 2, 0 }, // All faces have external normals
77 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
78 { 0, 2, 1, 0 }, // All faces have external normals
82 static int Tetra_nbN [] = { 3, 3, 3, 3 };
95 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
96 { 0, 1, 2, 3, 0 }, // All faces have external normals
102 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
103 { 0, 3, 2, 1, 0 }, // All faces but a bottom have external normals
108 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
119 // | / \ | PENTAHEDRON
125 static int Penta_F [5][5] = { // FORWARD
126 { 0, 1, 2, 0, 0 }, // All faces have external normals
127 { 3, 5, 4, 3, 3 }, // 0 is bottom, 1 is top face
131 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
137 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
144 // N4+----------+N7 |
145 // | | | | HEXAHEDRON
146 // | N1+------|---+N2
152 static int Hexa_F [6][5] = { // FORWARD
154 { 4, 7, 6, 5, 4 }, // all face normals are external
159 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
161 { 4, 5, 6, 7, 4 }, // all face normals are external
166 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
167 static int Hexa_oppF[] = { 1, 0, 4, 5, 2, 3 }; // oppopsite facet indices
186 static int HexPrism_F [8][7] = { // FORWARD
187 { 0, 1, 2, 3, 4, 5, 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_RE [8][7] = { // REVERSED -> EXTERNAL
196 { 0, 5, 4, 3, 2, 1, 0 },
197 { 6,11,10, 9, 8, 7, 6 },
198 { 0, 6, 7, 1, 0, 0, 0 },
199 { 1, 7, 8, 2, 1, 1, 1 },
200 { 2, 8, 9, 3, 2, 2, 2 },
201 { 3, 9,10, 4, 3, 3, 3 },
202 { 4,10,11, 5, 4, 4, 4 },
203 { 5,11, 6, 0, 5, 5, 5 }};
204 static int HexPrism_nbN [] = { 6, 6, 4, 4, 4, 4, 4, 4 };
213 // N0 +---|---+ N1 TETRAHEDRON
221 static int QuadTetra_F [4][7] = { // FORWARD
222 { 0, 4, 1, 5, 2, 6, 0 }, // All faces have external normals
223 { 0, 7, 3, 8, 1, 4, 0 },
224 { 1, 8, 3, 9, 2, 5, 1 },
225 { 0, 6, 2, 9, 3, 7, 0 }};
226 static int QuadTetra_RE [4][7] = { // REVERSED -> FORWARD (EXTERNAL)
227 { 0, 6, 2, 5, 1, 4, 0 }, // All faces have external normals
228 { 0, 4, 1, 8, 3, 7, 0 },
229 { 1, 5, 2, 9, 3, 8, 1 },
230 { 0, 7, 3, 9, 2, 6, 0 }};
231 static int QuadTetra_nbN [] = { 6, 6, 6, 6 };
241 // | | 9 - middle point for (0,4) etc.
254 static int QuadPyram_F [5][9] = { // FORWARD
255 { 0, 5, 1, 6, 2, 7, 3, 8, 0 }, // All faces have external normals
256 { 0, 9, 4, 10,1, 5, 0, 4, 4 },
257 { 1, 10,4, 11,2, 6, 1, 4, 4 },
258 { 2, 11,4, 12,3, 7, 2, 4, 4 },
259 { 3, 12,4, 9, 0, 8, 3, 4, 4 }};
260 static int QuadPyram_RE [5][9] = { // REVERSED -> FORWARD (EXTERNAL)
261 { 0, 8, 3, 7, 2, 6, 1, 5, 0 }, // All faces but a bottom have external normals
262 { 0, 5, 1, 10,4, 9, 0, 4, 4 },
263 { 1, 6, 2, 11,4, 10,1, 4, 4 },
264 { 2, 7, 3, 12,4, 11,2, 4, 4 },
265 { 3, 8, 0, 9, 4, 12,3, 4, 4 }};
266 static int QuadPyram_nbN [] = { 8, 6, 6, 6, 6 };
290 static int QuadPenta_F [5][9] = { // FORWARD
291 { 0, 6, 1, 7, 2, 8, 0, 0, 0 },
292 { 3, 11,5, 10,4, 9, 3, 3, 3 },
293 { 0, 12,3, 9, 4, 13,1, 6, 0 },
294 { 1, 13,4, 10,5, 14,2, 7, 1 },
295 { 0, 8, 2, 14,5, 11,3, 12,0 }};
296 static int QuadPenta_RE [5][9] = { // REVERSED -> EXTERNAL
297 { 0, 8, 2, 7, 1, 6, 0, 0, 0 },
298 { 3, 9, 4, 10,5, 11,3, 3, 3 },
299 { 0, 6, 1, 13,4, 9, 3, 12,0 },
300 { 1, 7, 2, 14,5, 10,4, 13,1 },
301 { 0, 12,3, 11,5, 14,2, 8, 0 }};
302 static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 };
306 // N5+-----+-----+N6 +-----+-----+
308 // 12+ | 14+ | + | +25 + |
310 // N4+-----+-----+N7 | QUADRATIC +-----+-----+ | Central nodes
311 // | | 15 | | HEXAHEDRON | | | | of tri-quadratic
312 // | | | | | | | | HEXAHEDRON
313 // | 17+ | +18 | + 22+ | +
315 // | | | | | + | 26+ | + |
317 // 16+ | +19 | + | +24 + |
320 // | N1+-----+-|---+N2 | +-----+-|---+
322 // | +8 | +10 | + 20+ | +
324 // N0+-----+-----+N3 +-----+-----+
327 static int QuadHexa_F [6][9] = { // FORWARD
328 { 0, 8, 1, 9, 2, 10,3, 11,0 }, // all face normals are external,
329 { 4, 15,7, 14,6, 13,5, 12,4 },
330 { 0, 16,4, 12,5, 17,1, 8, 0 },
331 { 1, 17,5, 13,6, 18,2, 9, 1 },
332 { 3, 10,2, 18,6, 14,7, 19,3 },
333 { 0, 11,3, 19,7, 15,4, 16,0 }};
334 static int QuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL
335 { 0, 11,3, 10,2, 9, 1, 8, 0 }, // all face normals are external
336 { 4, 12,5, 13,6, 14,7, 15,4 },
337 { 0, 8, 1, 17,5, 12,4, 16,0 },
338 { 1, 9, 2, 18,6, 13,5, 17,1 },
339 { 3, 19,7, 14,6, 18,2, 10,3 },
340 { 0, 16,4, 15,7, 19,3, 11,0 }};
341 static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 };
343 static int TriQuadHexa_F [6][9] = { // FORWARD
344 { 0, 8, 1, 9, 2, 10,3, 11, 20 }, // all face normals are external
345 { 4, 15,7, 14,6, 13,5, 12, 25 },
346 { 0, 16,4, 12,5, 17,1, 8, 21 },
347 { 1, 17,5, 13,6, 18,2, 9, 22 },
348 { 3, 10,2, 18,6, 14,7, 19, 23 },
349 { 0, 11,3, 19,7, 15,4, 16, 24 }};
350 static int TriQuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL
351 { 0, 11,3, 10,2, 9, 1, 8, 20 }, // opposite faces are neighbouring,
352 { 4, 12,5, 13,6, 14,7, 15, 25 }, // all face normals are external
353 { 0, 8, 1, 17,5, 12,4, 16, 21 },
354 { 1, 9, 2, 18,6, 13,5, 17, 22 },
355 { 3, 19,7, 14,6, 18,2, 10, 23 },
356 { 0, 16,4, 15,7, 19,3, 11, 24 }};
357 static int TriQuadHexa_nbN [] = { 9, 9, 9, 9, 9, 9 };
360 // ========================================================
361 // to perform some calculations without linkage to CASCADE
362 // ========================================================
367 XYZ() { x = 0; y = 0; z = 0; }
368 XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
369 XYZ( const XYZ& other ) { x = other.x; y = other.y; z = other.z; }
370 XYZ( const SMDS_MeshNode* n ) { x = n->X(); y = n->Y(); z = n->Z(); }
371 double* data() { return &x; }
372 inline XYZ operator-( const XYZ& other );
373 inline XYZ operator+( const XYZ& other );
374 inline XYZ Crossed( const XYZ& other );
375 inline XYZ operator-();
376 inline double Dot( const XYZ& other );
377 inline double Magnitude();
378 inline double SquareMagnitude();
379 inline XYZ Normalize();
381 inline XYZ XYZ::operator-( const XYZ& Right ) {
382 return XYZ(x - Right.x, y - Right.y, z - Right.z);
384 inline XYZ XYZ::operator-() {
385 return XYZ(-x,-y,-z);
387 inline XYZ XYZ::operator+( const XYZ& Right ) {
388 return XYZ(x + Right.x, y + Right.y, z + Right.z);
390 inline XYZ XYZ::Crossed( const XYZ& Right ) {
391 return XYZ (y * Right.z - z * Right.y,
392 z * Right.x - x * Right.z,
393 x * Right.y - y * Right.x);
395 inline double XYZ::Dot( const XYZ& Other ) {
396 return(x * Other.x + y * Other.y + z * Other.z);
398 inline double XYZ::Magnitude() {
399 return sqrt (x * x + y * y + z * z);
401 inline double XYZ::SquareMagnitude() {
402 return (x * x + y * y + z * z);
404 inline XYZ XYZ::Normalize() {
405 double magnitude = Magnitude();
406 if ( magnitude != 0.0 )
407 return XYZ(x /= magnitude,y /= magnitude,z /= magnitude );
412 //================================================================================
414 * \brief Return linear type corresponding to a quadratic one
416 //================================================================================
418 SMDS_VolumeTool::VolumeType quadToLinear(SMDS_VolumeTool::VolumeType quadType)
420 SMDS_VolumeTool::VolumeType linType = SMDS_VolumeTool::VolumeType( int(quadType)-4 );
421 const int nbCornersByQuad = SMDS_VolumeTool::NbCornerNodes( quadType );
422 if ( SMDS_VolumeTool::NbCornerNodes( linType ) == nbCornersByQuad )
426 for ( ; iLin < SMDS_VolumeTool::NB_VOLUME_TYPES; ++iLin )
427 if ( SMDS_VolumeTool::NbCornerNodes( SMDS_VolumeTool::VolumeType( iLin )) == nbCornersByQuad)
428 return SMDS_VolumeTool::VolumeType( iLin );
430 return SMDS_VolumeTool::UNKNOWN;
435 //================================================================================
437 * \brief Saver/restorer of a SMDS_VolumeTool::myCurFace
439 //================================================================================
441 struct SMDS_VolumeTool::SaveFacet
443 SMDS_VolumeTool::Facet mySaved;
444 SMDS_VolumeTool::Facet& myToRestore;
445 SaveFacet( SMDS_VolumeTool::Facet& facet ): myToRestore( facet )
448 mySaved.myNodes.swap( facet.myNodes );
452 if ( myToRestore.myIndex != mySaved.myIndex )
453 myToRestore = mySaved;
454 myToRestore.myNodes.swap( mySaved.myNodes );
458 //=======================================================================
459 //function : SMDS_VolumeTool
461 //=======================================================================
463 SMDS_VolumeTool::SMDS_VolumeTool ()
468 //=======================================================================
469 //function : SMDS_VolumeTool
471 //=======================================================================
473 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume,
474 const bool ignoreCentralNodes)
476 Set( theVolume, ignoreCentralNodes );
479 //=======================================================================
480 //function : SMDS_VolumeTool
482 //=======================================================================
484 SMDS_VolumeTool::~SMDS_VolumeTool()
486 myCurFace.myNodeIndices = NULL;
489 //=======================================================================
490 //function : SetVolume
491 //purpose : Set volume to iterate on
492 //=======================================================================
494 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume,
495 const bool ignoreCentralNodes,
496 const std::vector<const SMDS_MeshNode*>* otherNodes)
501 myIgnoreCentralNodes = ignoreCentralNodes;
505 myVolumeNodes.clear();
506 myPolyIndices.clear();
507 myPolyQuantities.clear();
508 myPolyFacetOri.clear();
511 myExternalFaces = false;
513 myAllFacesNodeIndices_F = 0;
514 myAllFacesNodeIndices_RE = 0;
515 myAllFacesNbNodes = 0;
517 myCurFace.myIndex = -1;
518 myCurFace.myNodeIndices = NULL;
519 myCurFace.myNodes.clear();
522 if ( !theVolume || theVolume->GetType() != SMDSAbs_Volume )
525 myVolume = theVolume;
526 myNbFaces = theVolume->NbFaces();
527 if ( myVolume->IsPoly() )
529 myPolyedre = SMDS_Mesh::DownCast<SMDS_MeshVolume>( myVolume );
530 myPolyFacetOri.resize( myNbFaces, 0 );
534 myVolumeNodes.resize( myVolume->NbNodes() );
537 if ( otherNodes->size() != myVolumeNodes.size() )
538 return ( myVolume = 0 );
539 for ( size_t i = 0; i < otherNodes->size(); ++i )
540 if ( ! ( myVolumeNodes[i] = (*otherNodes)[0] ))
541 return ( myVolume = 0 );
545 myVolumeNodes.assign( myVolume->begin_nodes(), myVolume->end_nodes() );
550 return ( myVolume = 0 );
554 // define volume orientation
556 if ( GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z ))
558 const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
559 int topNodeIndex = myVolume->NbCornerNodes() - 1;
560 while ( !IsLinked( 0, topNodeIndex, /*ignoreMediumNodes=*/true )) --topNodeIndex;
561 const SMDS_MeshNode* topNode = myVolumeNodes[ topNodeIndex ];
562 XYZ upDir (topNode->X() - botNode->X(),
563 topNode->Y() - botNode->Y(),
564 topNode->Z() - botNode->Z() );
565 myVolForward = ( botNormal.Dot( upDir ) < 0 );
568 myCurFace.myIndex = -1; // previous setFace(0) didn't take myVolForward into account
573 //=======================================================================
575 //purpose : Inverse volume
576 //=======================================================================
578 #define SWAP_NODES(nodes,i1,i2) \
580 const SMDS_MeshNode* tmp = nodes[ i1 ]; \
581 nodes[ i1 ] = nodes[ i2 ]; \
584 void SMDS_VolumeTool::Inverse ()
586 if ( !myVolume ) return;
588 if (myVolume->IsPoly()) {
589 MESSAGE("Warning: attempt to inverse polyhedral volume");
593 myVolForward = !myVolForward;
594 myCurFace.myIndex = -1;
596 // inverse top and bottom faces
597 switch ( myVolumeNodes.size() ) {
599 SWAP_NODES( myVolumeNodes, 1, 2 );
602 SWAP_NODES( myVolumeNodes, 1, 3 );
605 SWAP_NODES( myVolumeNodes, 1, 2 );
606 SWAP_NODES( myVolumeNodes, 4, 5 );
609 SWAP_NODES( myVolumeNodes, 1, 3 );
610 SWAP_NODES( myVolumeNodes, 5, 7 );
613 SWAP_NODES( myVolumeNodes, 1, 5 );
614 SWAP_NODES( myVolumeNodes, 2, 4 );
615 SWAP_NODES( myVolumeNodes, 7, 11 );
616 SWAP_NODES( myVolumeNodes, 8, 10 );
620 SWAP_NODES( myVolumeNodes, 1, 2 );
621 SWAP_NODES( myVolumeNodes, 4, 6 );
622 SWAP_NODES( myVolumeNodes, 8, 9 );
625 SWAP_NODES( myVolumeNodes, 1, 3 );
626 SWAP_NODES( myVolumeNodes, 5, 8 );
627 SWAP_NODES( myVolumeNodes, 6, 7 );
628 SWAP_NODES( myVolumeNodes, 10, 12 );
631 SWAP_NODES( myVolumeNodes, 1, 2 );
632 SWAP_NODES( myVolumeNodes, 4, 5 );
633 SWAP_NODES( myVolumeNodes, 6, 8 );
634 SWAP_NODES( myVolumeNodes, 9, 11 );
635 SWAP_NODES( myVolumeNodes, 13, 14 );
638 SWAP_NODES( myVolumeNodes, 1, 3 );
639 SWAP_NODES( myVolumeNodes, 5, 7 );
640 SWAP_NODES( myVolumeNodes, 8, 11 );
641 SWAP_NODES( myVolumeNodes, 9, 10 );
642 SWAP_NODES( myVolumeNodes, 12, 15 );
643 SWAP_NODES( myVolumeNodes, 13, 14 );
644 SWAP_NODES( myVolumeNodes, 17, 19 );
647 SWAP_NODES( myVolumeNodes, 1, 3 );
648 SWAP_NODES( myVolumeNodes, 5, 7 );
649 SWAP_NODES( myVolumeNodes, 8, 11 );
650 SWAP_NODES( myVolumeNodes, 9, 10 );
651 SWAP_NODES( myVolumeNodes, 12, 15 );
652 SWAP_NODES( myVolumeNodes, 13, 14 );
653 SWAP_NODES( myVolumeNodes, 17, 19 );
654 SWAP_NODES( myVolumeNodes, 21, 24 );
655 SWAP_NODES( myVolumeNodes, 22, 23 );
661 //=======================================================================
662 //function : GetVolumeType
664 //=======================================================================
666 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
671 switch( myVolumeNodes.size() ) {
672 case 4: return TETRA;
673 case 5: return PYRAM;
674 case 6: return PENTA;
676 case 12: return HEX_PRISM;
677 case 10: return QUAD_TETRA;
678 case 13: return QUAD_PYRAM;
679 case 15: return QUAD_PENTA;
680 case 20: return QUAD_HEXA;
681 case 27: return QUAD_HEXA;
688 //=======================================================================
689 //function : getTetraVolume
691 //=======================================================================
693 static double getTetraVolume(const SMDS_MeshNode* n1,
694 const SMDS_MeshNode* n2,
695 const SMDS_MeshNode* n3,
696 const SMDS_MeshNode* n4)
698 double p1[3], p2[3], p3[3], p4[3];
704 double Q1 = -(p1[ 0 ]-p2[ 0 ])*(p3[ 1 ]*p4[ 2 ]-p4[ 1 ]*p3[ 2 ]);
705 double Q2 = (p1[ 0 ]-p3[ 0 ])*(p2[ 1 ]*p4[ 2 ]-p4[ 1 ]*p2[ 2 ]);
706 double R1 = -(p1[ 0 ]-p4[ 0 ])*(p2[ 1 ]*p3[ 2 ]-p3[ 1 ]*p2[ 2 ]);
707 double R2 = -(p2[ 0 ]-p3[ 0 ])*(p1[ 1 ]*p4[ 2 ]-p4[ 1 ]*p1[ 2 ]);
708 double S1 = (p2[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p3[ 2 ]-p3[ 1 ]*p1[ 2 ]);
709 double S2 = -(p3[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p2[ 2 ]-p2[ 1 ]*p1[ 2 ]);
711 return (Q1+Q2+R1+R2+S1+S2)/6.0;
714 //=======================================================================
716 //purpose : Return element volume
717 //=======================================================================
718 double SMDS_VolumeTool::GetSize() const
724 if ( myVolume->IsPoly() )
729 SaveFacet savedFacet( myCurFace );
731 // split a polyhedron into tetrahedrons
733 bool oriOk = AllFacesSameOriented();
734 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
735 for ( int f = 0; f < NbFaces(); ++f )
738 XYZ area (0,0,0), p1( myCurFace.myNodes[0] );
739 for ( int n = 0; n < myCurFace.myNbNodes; ++n )
741 XYZ p2( myCurFace.myNodes[ n+1 ]);
742 area = area + p1.Crossed( p2 );
748 if ( !oriOk && V > 0 )
753 const static int ind[] = {
754 0, 1, 3, 6, 11, 23, 31, 44, 58, 78 };
755 const static int vtab[][4] = { // decomposition into tetra in the order of enum VolumeType
788 // quadratic tetrahedron
813 // quadratic pentahedron
830 // quadratic hexahedron
855 int type = GetVolumeType();
857 int n2 = ind[type+1];
859 for (int i = n1; i < n2; i++) {
860 V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
861 myVolumeNodes[ vtab[i][1] ],
862 myVolumeNodes[ vtab[i][2] ],
863 myVolumeNodes[ vtab[i][3] ]);
870 //=======================================================================
871 //function : getTetraScaledJacobian
872 //purpose : Given the smesh nodes in the canonical order of the tetrahedron, return the scaled jacobian
873 //=======================================================================
874 static double getTetraScaledJacobian(const SMDS_MeshNode* n0,
875 const SMDS_MeshNode* n1,
876 const SMDS_MeshNode* n2,
877 const SMDS_MeshNode* n3)
879 const double sqrt = std::sqrt(2.0);
880 // Get the coordinates
885 // Define the edges connecting the nodes
888 XYZ L2 = p2-p0; // invert the definition of doc to get the proper orientation of the crossed product
892 double Jacobian = L2.Crossed( L0 ).Dot( L3 );
893 double norm0 = L0.Magnitude();
894 double norm1 = L1.Magnitude();
895 double norm2 = L2.Magnitude();
896 double norm3 = L3.Magnitude();
897 double norm4 = L4.Magnitude();
898 double norm5 = L5.Magnitude();
900 std::array<double, 5> norms{};
902 norms[1] = norm3*norm4*norm5;
903 norms[2] = norm1*norm2*norm5;
904 norms[3] = norm0*norm1*norm4;
905 norms[4] = norm0*norm2*norm3;
907 auto findMaxNorm = std::max_element(norms.begin(), norms.end());
908 double maxNorm = *findMaxNorm;
910 if ( std::fabs( maxNorm ) < std::numeric_limits<double>::min() )
911 maxNorm = std::numeric_limits<double>::max();
913 return Jacobian * sqrt / maxNorm;
916 //=======================================================================
917 //function : getPyramidScaledJacobian
918 //purpose : Given the pyramid, compute the scaled jacobian of the four tetrahedrons and return the minimun value.
919 //=======================================================================
920 static double getPyramidScaledJacobian(const SMDS_MeshNode* n0,
921 const SMDS_MeshNode* n1,
922 const SMDS_MeshNode* n2,
923 const SMDS_MeshNode* n3,
924 const SMDS_MeshNode* n4)
926 const double sqrt = std::sqrt(2.0);
927 std::array<double, 4> tetScaledJacobian{};
928 tetScaledJacobian[0] = getTetraScaledJacobian(n0, n1, n3, n4);
929 tetScaledJacobian[1] = getTetraScaledJacobian(n1, n2, n0, n4);
930 tetScaledJacobian[2] = getTetraScaledJacobian(n2, n3, n1, n4);
931 tetScaledJacobian[3] = getTetraScaledJacobian(n3, n0, n2, n4);
933 auto minEntry = std::min_element(tetScaledJacobian.begin(), tetScaledJacobian.end());
935 double scaledJacobian = (*minEntry) * 2.0/sqrt;
936 return scaledJacobian < 1.0 ? scaledJacobian : 1.0 - (scaledJacobian - 1.0);
941 //=======================================================================
942 //function : getHexaScaledJacobian
943 //purpose : Evaluate the scaled jacobian on the eight vertices of the hexahedron and return the minimal registered value
944 //remark : Follow the reference numeration described at the top of the class.
945 //=======================================================================
946 static double getHexaScaledJacobian(const SMDS_MeshNode* n0,
947 const SMDS_MeshNode* n1,
948 const SMDS_MeshNode* n2,
949 const SMDS_MeshNode* n3,
950 const SMDS_MeshNode* n4,
951 const SMDS_MeshNode* n5,
952 const SMDS_MeshNode* n6,
953 const SMDS_MeshNode* n7)
955 // Scaled jacobian is an scalar quantity measuring the deviation of the geometry from the perfect geometry
956 // Get the coordinates
966 // Define the edges connecting the nodes
967 XYZ L0 = (p1-p0).Normalize();
968 XYZ L1 = (p2-p1).Normalize();
969 XYZ L2 = (p3-p2).Normalize();
970 XYZ L3 = (p3-p0).Normalize();
971 XYZ L4 = (p4-p0).Normalize();
972 XYZ L5 = (p5-p1).Normalize();
973 XYZ L6 = (p6-p2).Normalize();
974 XYZ L7 = (p7-p3).Normalize();
975 XYZ L8 = (p5-p4).Normalize();
976 XYZ L9 = (p6-p5).Normalize();
977 XYZ L10 = (p7-p6).Normalize();
978 XYZ L11 = (p7-p4).Normalize();
979 XYZ X0 = (p1-p0+p2-p3+p6-p7+p5-p4).Normalize();
980 XYZ X1 = (p3-p0+p2-p1+p7-p4+p6-p5).Normalize();
981 XYZ X2 = (p4-p0+p7-p3+p5-p1+p6-p2).Normalize();
983 std::array<double, 9> scaledJacobian{};
984 //Scaled jacobian of nodes following their numeration
985 scaledJacobian[0] = L4.Crossed( L3).Dot( L0 ); // For L0
986 scaledJacobian[1] = L5.Crossed(-L0).Dot( L1 ); // For L1
987 scaledJacobian[2] = L6.Crossed(-L1).Dot( L2 ); // For L2
988 scaledJacobian[3] = L7.Crossed(-L2).Dot(-L3 ); // For L3
989 scaledJacobian[4] = -L4.Crossed( L8).Dot( L11 ); // For L11
990 scaledJacobian[5] = -L5.Crossed( L9).Dot(-L8 ); // For L8
991 scaledJacobian[6] = -L6.Crossed(L10).Dot(-L9 ); // For L9
992 scaledJacobian[7] = -L7.Crossed(-L11).Dot(-L10 ); // For L10
993 scaledJacobian[8] = X2.Crossed( X1).Dot( X0 ); // For principal axes
995 auto minScaledJacobian = std::min_element(scaledJacobian.begin(), scaledJacobian.end());
996 return *minScaledJacobian;
1000 //=======================================================================
1001 //function : getTetraNormalizedJacobian
1002 //purpose : Return the jacobian of the tetrahedron based on normalized vectors
1003 //=======================================================================
1004 static double getTetraNormalizedJacobian(const SMDS_MeshNode* n0,
1005 const SMDS_MeshNode* n1,
1006 const SMDS_MeshNode* n2,
1007 const SMDS_MeshNode* n3)
1009 const double sqrt = std::sqrt(2.0);
1010 // Get the coordinates
1015 // Define the normalized edges connecting the nodes
1016 XYZ L0 = (p1-p0).Normalize();
1017 XYZ L2 = (p2-p0).Normalize(); // invert the definition of doc to get the proper orientation of the crossed product
1018 XYZ L3 = (p3-p0).Normalize();
1019 return L2.Crossed( L0 ).Dot( L3 );
1022 //=======================================================================
1023 //function : getPentaScaledJacobian
1024 //purpose : Evaluate the scaled jacobian on the pentahedron based on decomposed tetrahedrons
1025 //=======================================================================
1032 // N0 +---------+ N2
1033 // | | | NUMERATION RERENCE FOLLOWING POSSITIVE RIGHT HAND RULE
1035 // | / \ | PENTAHEDRON
1039 // N3 +---------+ N5
1046 // N0 +--|---+ N2 TETRAHEDRON ASSOCIATED TO N0
1047 // \ | / Numeration passed to getTetraScaledJacobian
1048 // \ | / N0=N0; N1=N2; N2=N3; N3=N1
1059 // N2 +---|---+ N5 TETRAHEDRON ASSOCIATED TO N2
1060 // \ | / Numeration passed to getTetraScaledJacobian
1061 // \ | / N0=N2; N1=N5; N2=N0; N3=N1
1072 // N3 +---|---+ N0 TETRAHEDRON ASSOCIATED TO N3
1073 // \ | / Numeration passed to getTetraScaledJacobian
1074 // \ | / N0=N3; N1=N0; N2=N5; N3=N4
1085 // N1 +---|---+ N2 TETRAHEDRON ASSOCIATED TO N1
1086 // \ | / Numeration passed to getTetraScaledJacobian
1087 // \ | / N0=N1; N1=N2; N2=N0; N3=N3
1095 static double getPentaScaledJacobian(const SMDS_MeshNode* n0,
1096 const SMDS_MeshNode* n1,
1097 const SMDS_MeshNode* n2,
1098 const SMDS_MeshNode* n3,
1099 const SMDS_MeshNode* n4,
1100 const SMDS_MeshNode* n5)
1102 std::array<double, 6> scaledJacobianOfReferenceTetra{};
1103 scaledJacobianOfReferenceTetra[0] = getTetraNormalizedJacobian(n0, n2, n3, n1); // For n0
1104 scaledJacobianOfReferenceTetra[1] = getTetraNormalizedJacobian(n2, n5, n0, n1); // For n2
1105 scaledJacobianOfReferenceTetra[2] = getTetraNormalizedJacobian(n3, n0, n5, n4); // For n3
1106 scaledJacobianOfReferenceTetra[3] = getTetraNormalizedJacobian(n5, n3, n2, n4); // For n5
1107 scaledJacobianOfReferenceTetra[4] = getTetraNormalizedJacobian(n1, n2, n0, n3); // For n1
1108 scaledJacobianOfReferenceTetra[5] = getTetraNormalizedJacobian(n4, n3, n5, n2); // For n4
1110 auto minScaledJacobian = std::min_element(scaledJacobianOfReferenceTetra.begin(), scaledJacobianOfReferenceTetra.end());
1111 double minScalJac = (*minScaledJacobian)* 2.0 / std::sqrt(3.0);
1114 return std::min(minScalJac, std::numeric_limits<double>::max());
1116 return std::max(minScalJac, -std::numeric_limits<double>::max());
1119 //=======================================================================
1120 //function : getHexaPrismScaledJacobian
1121 //purpose : Evaluate the scaled jacobian on the hexaprism by decomposing the goemetry into three 1hexa + 2 pentahedrons
1122 //=======================================================================
1123 static double getHexaPrismScaledJacobian(const SMDS_MeshNode* n0,
1124 const SMDS_MeshNode* n1,
1125 const SMDS_MeshNode* n2,
1126 const SMDS_MeshNode* n3,
1127 const SMDS_MeshNode* n4,
1128 const SMDS_MeshNode* n5,
1129 const SMDS_MeshNode* n6,
1130 const SMDS_MeshNode* n7,
1131 const SMDS_MeshNode* n8,
1132 const SMDS_MeshNode* n9,
1133 const SMDS_MeshNode* n10,
1134 const SMDS_MeshNode* n11)
1136 // The Pentahedron from the left
1137 // n0=n0; n1=n1; n2=n2; n3=n6; n4=n7, n5=n8;
1138 double scaledJacobianPentleft = getPentaScaledJacobian( n0, n1, n2, n6, n7, n8 );
1139 // The core Hexahedron
1140 // n0=n0; n1=n2, n2=n3; n3=n5; n4=n6; n5=n8; n6=n9; n7=n11
1141 double scaledJacobianHexa = getHexaScaledJacobian( n0, n2, n3, n5, n6, n8, n9, n11 );
1142 // The Pentahedron from the right
1143 // n0=n5; n1=n4; n2=n3; n3=n11; n4=n10; n5=n9
1144 double scaledJacobianPentright = getPentaScaledJacobian( n5, n4, n3, n11, n10, n9 );
1146 return std::min( scaledJacobianHexa, std::min( scaledJacobianPentleft, scaledJacobianPentright ) );
1150 //=======================================================================
1151 //function : GetScaledJacobian
1152 //purpose : Return element Scaled Jacobian using the generic definition given
1153 // in https://gitlab.kitware.com/third-party/verdict/-/blob/master/SAND2007-2853p.pdf
1154 //=======================================================================
1156 double SMDS_VolumeTool::GetScaledJacobian() const
1159 // For Tetra, call directly the getTetraScaledJacobian
1160 double scaledJacobian = 0.;
1162 VolumeType type = GetVolumeType();
1167 scaledJacobian = getTetraScaledJacobian( myVolumeNodes[0], myVolumeNodes[1], myVolumeNodes[2], myVolumeNodes[3] );
1171 scaledJacobian = getHexaScaledJacobian( myVolumeNodes[0], myVolumeNodes[1], myVolumeNodes[2], myVolumeNodes[3],
1172 myVolumeNodes[4], myVolumeNodes[5], myVolumeNodes[6], myVolumeNodes[7] );
1176 scaledJacobian = getPyramidScaledJacobian( myVolumeNodes[0], myVolumeNodes[1], myVolumeNodes[2], myVolumeNodes[3], myVolumeNodes[4] );
1180 scaledJacobian = getPentaScaledJacobian( myVolumeNodes[0], myVolumeNodes[1],
1181 myVolumeNodes[2], myVolumeNodes[3],
1182 myVolumeNodes[4], myVolumeNodes[5] );
1185 scaledJacobian = getHexaPrismScaledJacobian( myVolumeNodes[0], myVolumeNodes[1], myVolumeNodes[2], myVolumeNodes[3],
1186 myVolumeNodes[4], myVolumeNodes[5], myVolumeNodes[6], myVolumeNodes[7],
1187 myVolumeNodes[8], myVolumeNodes[9], myVolumeNodes[10], myVolumeNodes[11]);
1193 return scaledJacobian;
1197 //=======================================================================
1198 //function : GetBaryCenter
1200 //=======================================================================
1202 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
1208 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
1209 X += myVolumeNodes[ i ]->X();
1210 Y += myVolumeNodes[ i ]->Y();
1211 Z += myVolumeNodes[ i ]->Z();
1213 X /= myVolumeNodes.size();
1214 Y /= myVolumeNodes.size();
1215 Z /= myVolumeNodes.size();
1220 //================================================================================
1222 * \brief Classify a point
1223 * \param tol - thickness of faces
1225 //================================================================================
1227 bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
1229 // LIMITATION: for convex volumes only
1231 for ( int iF = 0; iF < myNbFaces; ++iF )
1234 if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z ))
1236 if ( !IsFaceExternal( iF ))
1237 faceNormal = XYZ() - faceNormal; // reverse
1239 XYZ face2p( p - XYZ( myCurFace.myNodes[0] ));
1240 if ( face2p.Dot( faceNormal ) > tol )
1246 //=======================================================================
1247 //function : SetExternalNormal
1248 //purpose : Node order will be so that faces normals are external
1249 //=======================================================================
1251 void SMDS_VolumeTool::SetExternalNormal ()
1253 myExternalFaces = true;
1254 myCurFace.myIndex = -1;
1257 //=======================================================================
1258 //function : NbFaceNodes
1259 //purpose : Return number of nodes in the array of face nodes
1260 //=======================================================================
1262 int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const
1264 if ( !setFace( faceIndex ))
1266 return myCurFace.myNbNodes;
1269 //=======================================================================
1270 //function : GetFaceNodes
1271 //purpose : Return pointer to the array of face nodes.
1272 // To comfort link iteration, the array
1273 // length == NbFaceNodes( faceIndex ) + 1 and
1274 // the last node == the first one.
1275 //=======================================================================
1277 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const
1279 if ( !setFace( faceIndex ))
1281 return &myCurFace.myNodes[0];
1284 //=======================================================================
1285 //function : GetFaceNodesIndices
1286 //purpose : Return pointer to the array of face nodes indices
1287 // To comfort link iteration, the array
1288 // length == NbFaceNodes( faceIndex ) + 1 and
1289 // the last node index == the first one.
1290 //=======================================================================
1292 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const
1294 if ( !setFace( faceIndex ))
1297 return myCurFace.myNodeIndices;
1300 //=======================================================================
1301 //function : GetFaceNodes
1302 //purpose : Return a set of face nodes.
1303 //=======================================================================
1305 bool SMDS_VolumeTool::GetFaceNodes (int faceIndex,
1306 std::set<const SMDS_MeshNode*>& theFaceNodes ) const
1308 if ( !setFace( faceIndex ))
1311 theFaceNodes.clear();
1312 theFaceNodes.insert( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
1319 struct NLink : public std::pair<smIdType,smIdType>
1322 NLink(const SMDS_MeshNode* n1=0, const SMDS_MeshNode* n2=0, int ori=1 )
1326 if (( myOri = ( n1->GetID() < n2->GetID() )))
1328 first = n1->GetID();
1329 second = n2->GetID();
1334 first = n2->GetID();
1335 second = n1->GetID();
1341 myOri = first = second = 0;
1344 //int Node1() const { return myOri == -1 ? second : first; }
1346 //bool IsSameOri( const std::pair<int,int>& link ) const { return link.first == Node1(); }
1350 //=======================================================================
1351 //function : IsFaceExternal
1352 //purpose : Check normal orientation of a given face
1353 //=======================================================================
1355 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const
1357 if ( myExternalFaces || !myVolume )
1360 if ( !myPolyedre ) // all classical volumes have external facet normals
1363 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1365 if ( myPolyFacetOri[ faceIndex ])
1366 return myPolyFacetOri[ faceIndex ] > 0;
1368 int ori = 0; // -1-in, +1-out, 0-undef
1369 double minProj, maxProj;
1370 if ( projectNodesToNormal( faceIndex, minProj, maxProj ))
1372 // all nodes are on the same side of the facet
1373 ori = ( minProj < 0 ? +1 : -1 );
1374 me->myPolyFacetOri[ faceIndex ] = ori;
1376 if ( !myFwdLinks.empty() ) // concave polyhedron; collect oriented links
1377 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1379 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1380 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1385 SaveFacet savedFacet( myCurFace );
1387 // concave polyhedron
1389 if ( myFwdLinks.empty() ) // get links of the least ambiguously oriented facet
1391 for ( size_t i = 0; i < myPolyFacetOri.size() && !ori; ++i )
1392 ori = myPolyFacetOri[ i ];
1394 if ( !ori ) // none facet is oriented yet
1396 // find the least ambiguously oriented facet
1397 int faceMostConvex = -1;
1398 std::map< double, int > convexity2face;
1399 for ( size_t iF = 0; iF < myPolyFacetOri.size() && faceMostConvex < 0; ++iF )
1401 if ( projectNodesToNormal( iF, minProj, maxProj ))
1403 // all nodes are on the same side of the facet
1404 me->myPolyFacetOri[ iF ] = ( minProj < 0 ? +1 : -1 );
1405 faceMostConvex = iF;
1409 ori = ( -minProj < maxProj ? -1 : +1 );
1410 double convexity = std::min( -minProj, maxProj ) / std::max( -minProj, maxProj );
1411 convexity2face.insert( std::make_pair( convexity, iF * ori ));
1414 if ( faceMostConvex < 0 ) // none facet has nodes on the same side
1416 // use the least ambiguous facet
1417 faceMostConvex = convexity2face.begin()->second;
1418 ori = ( faceMostConvex < 0 ? -1 : +1 );
1419 faceMostConvex = std::abs( faceMostConvex );
1420 me->myPolyFacetOri[ faceMostConvex ] = ori;
1423 // collect links of the oriented facets in myFwdLinks
1424 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1426 ori = myPolyFacetOri[ iF ];
1427 if ( !ori ) continue;
1429 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1431 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1432 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1437 // compare orientation of links of the facet with myFwdLinks
1439 setFace( faceIndex );
1440 std::vector< NLink > links( myCurFace.myNbNodes ), links2;
1441 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1443 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1444 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1445 if ( l2o != myFwdLinks.end() )
1446 ori = link.myOri * l2o->second * -1;
1449 while ( !ori ) // the facet has no common links with already oriented facets
1451 // orient and collect links of other non-oriented facets
1452 for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1454 if ( myPolyFacetOri[ iF ] ) continue; // already oriented
1458 for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1460 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1461 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1462 if ( l2o != myFwdLinks.end() )
1463 ori = link.myOri * l2o->second * -1;
1464 links2.push_back( link );
1466 if ( ori ) // one more facet oriented
1468 me->myPolyFacetOri[ iF ] = ori;
1469 for ( size_t i = 0; i < links2.size(); ++i )
1470 me->myFwdLinks.insert( make_pair( links2[i], links2[i].myOri * ori ));
1475 return false; // error in algorithm: infinite loop
1477 // try to orient the facet again
1479 for ( size_t i = 0; i < links.size() && !ori; ++i )
1481 std::map<Link, int>::const_iterator l2o = myFwdLinks.find( links[i] );
1482 if ( l2o != myFwdLinks.end() )
1483 ori = links[i].myOri * l2o->second * -1;
1485 me->myPolyFacetOri[ faceIndex ] = ori;
1491 //=======================================================================
1492 //function : projectNodesToNormal
1493 //purpose : compute min and max projections of all nodes to normal of a facet.
1494 //=======================================================================
1496 bool SMDS_VolumeTool::projectNodesToNormal( int faceIndex,
1499 double* normalXYZ ) const
1501 minProj = std::numeric_limits<double>::max();
1502 maxProj = std::numeric_limits<double>::min();
1505 if ( !GetFaceNormal( faceIndex, normal.x, normal.y, normal.z ))
1508 memcpy( normalXYZ, normal.data(), 3*sizeof(double));
1510 XYZ p0 ( myCurFace.myNodes[0] );
1511 for ( size_t i = 0; i < myVolumeNodes.size(); ++i )
1513 if ( std::find( myCurFace.myNodes.begin() + 1,
1514 myCurFace.myNodes.end(),
1515 myVolumeNodes[ i ] ) != myCurFace.myNodes.end() )
1516 continue; // node of the faceIndex-th facet
1518 double proj = normal.Dot( XYZ( myVolumeNodes[ i ]) - p0 );
1519 if ( proj < minProj ) minProj = proj;
1520 if ( proj > maxProj ) maxProj = proj;
1522 const double tol = 1e-7;
1525 bool diffSize = ( minProj * maxProj < 0 );
1528 // minProj = -minProj;
1530 // else if ( minProj < 0 )
1532 // minProj = -minProj;
1533 // maxProj = -maxProj;
1536 return !diffSize; // ? 0 : (minProj >= 0);
1539 //=======================================================================
1540 //function : GetFaceNormal
1541 //purpose : Return a normal to a face
1542 //=======================================================================
1544 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) const
1546 if ( !setFace( faceIndex ))
1549 const int iQuad = ( !myPolyedre && myCurFace.myNbNodes > 6 ) ? 2 : 1;
1550 XYZ p1 ( myCurFace.myNodes[0*iQuad] );
1551 XYZ p2 ( myCurFace.myNodes[1*iQuad] );
1552 XYZ p3 ( myCurFace.myNodes[2*iQuad] );
1553 XYZ aVec12( p2 - p1 );
1554 XYZ aVec13( p3 - p1 );
1555 XYZ cross = aVec12.Crossed( aVec13 );
1557 for ( int i = 3*iQuad; i < myCurFace.myNbNodes; i += iQuad )
1559 XYZ p4 ( myCurFace.myNodes[i] );
1560 XYZ aVec14( p4 - p1 );
1561 XYZ cross2 = aVec13.Crossed( aVec14 );
1562 cross = cross + cross2;
1566 double size = cross.Magnitude();
1567 if ( size <= std::numeric_limits<double>::min() )
1577 //================================================================================
1579 * \brief Return barycenter of a face
1581 //================================================================================
1583 bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const
1585 if ( !setFace( faceIndex ))
1589 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1591 X += myCurFace.myNodes[i]->X() / myCurFace.myNbNodes;
1592 Y += myCurFace.myNodes[i]->Y() / myCurFace.myNbNodes;
1593 Z += myCurFace.myNodes[i]->Z() / myCurFace.myNbNodes;
1598 //================================================================================
1600 * \brief Check that all the faces in a polyhedron follow the same orientation
1601 * \remark there is no differentiation for outward and inward face orientation.
1603 //================================================================================
1604 bool SMDS_VolumeTool::AllFacesSameOriented() const
1606 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
1607 bool validOrientation = true;
1608 std::map<Link, std::vector<int>> collectLinksOrientations;
1609 me->myFwdLinks.clear();
1610 for ( int faceId = 0; faceId < NbFaces(); ++faceId )
1612 me->setFace( faceId );
1613 myExternalFaces = false;
1616 for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1618 NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1619 std::map<Link, int>::const_iterator foundLink = myFwdLinks.find( link );
1621 if ( foundLink == myFwdLinks.end() )
1622 me->myFwdLinks.insert( make_pair( link, link.myOri ));
1624 collectLinksOrientations[ link ].push_back( link.myOri );
1628 // Check duality of the orientations
1629 std::map<Link, std::vector<int>>::const_iterator theLinks;
1630 for ( theLinks = collectLinksOrientations.begin(); theLinks != collectLinksOrientations.end(); theLinks++ )
1632 if ( theLinks->second.size() == 2 ) // 99% of the cases there are 2 faces per link
1634 if ( 1 != -1*theLinks->second[0]*theLinks->second[1] )
1639 if ( theLinks->second.size() % 2 != 0 )// Dont treat uneven number of links
1642 // In the other 1% of the cases we count the number occurrence and check that they match
1643 int minusOne = std::count( theLinks->second.begin(), theLinks->second.end(), -1 );
1644 int plusOne = std::count( theLinks->second.begin(), theLinks->second.end(), 1 );
1645 if ( minusOne != plusOne )
1649 return validOrientation;
1652 //=======================================================================
1653 //function : GetFaceArea
1654 //purpose : Return face area
1655 //=======================================================================
1657 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
1660 if ( !setFace( faceIndex ))
1663 XYZ p1 ( myCurFace.myNodes[0] );
1664 XYZ p2 ( myCurFace.myNodes[1] );
1665 XYZ p3 ( myCurFace.myNodes[2] );
1666 XYZ aVec12( p2 - p1 );
1667 XYZ aVec13( p3 - p1 );
1668 area += aVec12.Crossed( aVec13 ).Magnitude();
1670 if (myVolume->IsPoly())
1672 for ( int i = 3; i < myCurFace.myNbNodes; ++i )
1674 XYZ pI ( myCurFace.myNodes[i] );
1675 XYZ aVecI( pI - p1 );
1676 area += aVec13.Crossed( aVecI ).Magnitude();
1682 if ( myCurFace.myNbNodes == 4 ) {
1683 XYZ p4 ( myCurFace.myNodes[3] );
1684 XYZ aVec14( p4 - p1 );
1685 area += aVec14.Crossed( aVec13 ).Magnitude();
1691 //================================================================================
1693 * \brief Return index of the node located at face center of a quadratic element like HEX27
1695 //================================================================================
1697 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
1699 if ( myAllFacesNbNodes && myVolumeNodes.size() == 27 ) // classic element with 27 nodes
1701 switch ( faceIndex ) {
1705 return faceIndex + 19;
1711 //=======================================================================
1712 //function : GetOppFaceIndex
1713 //purpose : Return index of the opposite face if it exists, else -1.
1714 //=======================================================================
1716 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1720 MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1724 const int nbHoriFaces = 2;
1726 if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1727 switch ( myVolumeNodes.size() ) {
1730 if ( faceIndex == 0 || faceIndex == 1 )
1731 ind = 1 - faceIndex;
1735 if ( faceIndex <= 1 ) // top or bottom
1736 ind = 1 - faceIndex;
1738 const int nbSideFaces = myAllFacesNbNodes[0];
1739 ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1744 ind = GetOppFaceIndexOfHex( faceIndex );
1752 //=======================================================================
1753 //function : GetOppFaceIndexOfHex
1754 //purpose : Return index of the opposite face of the hexahedron
1755 //=======================================================================
1757 int SMDS_VolumeTool::GetOppFaceIndexOfHex( int faceIndex )
1759 return Hexa_oppF[ faceIndex ];
1762 //=======================================================================
1763 //function : IsLinked
1764 //purpose : return true if theNode1 is linked with theNode2
1765 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1766 //=======================================================================
1768 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1769 const SMDS_MeshNode* theNode2,
1770 const bool theIgnoreMediumNodes) const
1775 if (myVolume->IsPoly()) {
1777 MESSAGE("Warning: bad volumic element");
1780 if ( !myAllFacesNbNodes ) {
1781 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1782 me->myPolyQuantities = myPolyedre->GetQuantities();
1783 myAllFacesNbNodes = &myPolyQuantities[0];
1785 int from, to = 0, d1 = 1, d2 = 2;
1786 if ( myPolyedre->IsQuadratic() ) {
1787 if ( theIgnoreMediumNodes ) {
1793 std::vector<const SMDS_MeshNode*>::const_iterator i;
1794 for (int iface = 0; iface < myNbFaces; iface++)
1797 to += myPolyQuantities[iface];
1798 i = std::find( myVolumeNodes.begin() + from, myVolumeNodes.begin() + to, theNode1 );
1799 if ( i != myVolumeNodes.end() )
1801 if (( theNode2 == *( i-d1 ) ||
1802 theNode2 == *( i+d1 )))
1805 (( theNode2 == *( i-d2 ) ||
1806 theNode2 == *( i+d2 ))))
1813 // find nodes indices
1814 int i1 = -1, i2 = -1, nbFound = 0;
1815 for ( size_t i = 0; i < myVolumeNodes.size() && nbFound < 2; i++ )
1817 if ( myVolumeNodes[ i ] == theNode1 )
1819 else if ( myVolumeNodes[ i ] == theNode2 )
1822 return IsLinked( i1, i2 );
1825 //=======================================================================
1826 //function : IsLinked
1827 //purpose : return true if the node with theNode1Index is linked
1828 // with the node with theNode2Index
1829 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1830 //=======================================================================
1832 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1833 const int theNode2Index,
1834 bool theIgnoreMediumNodes) const
1836 if ( myVolume->IsPoly() ) {
1837 return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1840 int minInd = std::min( theNode1Index, theNode2Index );
1841 int maxInd = std::max( theNode1Index, theNode2Index );
1843 if ( minInd < 0 || maxInd > (int)myVolumeNodes.size() - 1 || maxInd == minInd )
1846 VolumeType type = GetVolumeType();
1847 if ( myVolume->IsQuadratic() )
1849 int firstMediumInd = myVolume->NbCornerNodes();
1850 if ( minInd >= firstMediumInd )
1851 return false; // both nodes are medium - not linked
1852 if ( maxInd < firstMediumInd ) // both nodes are corners
1854 if ( theIgnoreMediumNodes )
1855 type = quadToLinear(type); // to check linkage of corner nodes only
1857 return false; // corner nodes are not linked directly in a quadratic cell
1865 switch ( maxInd - minInd ) {
1866 case 1: return minInd != 3;
1867 case 3: return minInd == 0 || minInd == 4;
1868 case 4: return true;
1875 switch ( maxInd - minInd ) {
1877 case 3: return true;
1882 switch ( maxInd - minInd ) {
1883 case 1: return minInd != 2;
1884 case 2: return minInd == 0 || minInd == 3;
1885 case 3: return true;
1892 case 0: return ( maxInd==4 || maxInd==6 || maxInd==7 );
1893 case 1: return ( maxInd==4 || maxInd==5 || maxInd==8 );
1894 case 2: return ( maxInd==5 || maxInd==6 || maxInd==9 );
1895 case 3: return ( maxInd==7 || maxInd==8 || maxInd==9 );
1903 case 0: return ( maxInd==8 || maxInd==11 || maxInd==16 );
1904 case 1: return ( maxInd==8 || maxInd==9 || maxInd==17 );
1905 case 2: return ( maxInd==9 || maxInd==10 || maxInd==18 );
1906 case 3: return ( maxInd==10 || maxInd==11 || maxInd==19 );
1907 case 4: return ( maxInd==12 || maxInd==15 || maxInd==16 );
1908 case 5: return ( maxInd==12 || maxInd==13 || maxInd==17 );
1909 case 6: return ( maxInd==13 || maxInd==14 || maxInd==18 );
1910 case 7: return ( maxInd==14 || maxInd==15 || maxInd==19 );
1918 case 0: return ( maxInd==5 || maxInd==8 || maxInd==9 );
1919 case 1: return ( maxInd==5 || maxInd==6 || maxInd==10 );
1920 case 2: return ( maxInd==6 || maxInd==7 || maxInd==11 );
1921 case 3: return ( maxInd==7 || maxInd==8 || maxInd==12 );
1922 case 4: return ( maxInd==9 || maxInd==10 || maxInd==11 || maxInd==12 );
1930 case 0: return ( maxInd==6 || maxInd==8 || maxInd==12 );
1931 case 1: return ( maxInd==6 || maxInd==7 || maxInd==13 );
1932 case 2: return ( maxInd==7 || maxInd==8 || maxInd==14 );
1933 case 3: return ( maxInd==9 || maxInd==11 || maxInd==12 );
1934 case 4: return ( maxInd==9 || maxInd==10 || maxInd==13 );
1935 case 5: return ( maxInd==10 || maxInd==11 || maxInd==14 );
1942 const int diff = maxInd-minInd;
1943 if ( diff > 6 ) return false;// not linked top and bottom
1944 if ( diff == 6 ) return true; // linked top and bottom
1945 return diff == 1 || diff == 7;
1952 //=======================================================================
1953 //function : GetNodeIndex
1954 //purpose : Return an index of theNode
1955 //=======================================================================
1957 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1960 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
1961 if ( myVolumeNodes[ i ] == theNode )
1968 //================================================================================
1970 * \brief Fill vector with boundary faces existing in the mesh
1971 * \param faces - vector of found nodes
1972 * \retval int - nb of found faces
1974 //================================================================================
1976 int SMDS_VolumeTool::GetAllExistingFaces(std::vector<const SMDS_MeshElement*> & faces) const
1979 SaveFacet savedFacet( myCurFace );
1981 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1983 if ( const SMDS_MeshElement* face = SMDS_Mesh::FindFace( myCurFace.myNodes ))
1984 faces.push_back( face );
1987 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1988 const SMDS_MeshFace* face = 0;
1989 const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1990 switch ( NbFaceNodes( iF )) {
1992 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1994 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1996 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1997 nodes[3], nodes[4], nodes[5]); break;
1999 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
2000 nodes[4], nodes[5], nodes[6], nodes[7]); break;
2003 faces.push_back( face );
2005 return faces.size();
2009 //================================================================================
2011 * \brief Fill vector with boundary edges existing in the mesh
2012 * \param edges - vector of found edges
2013 * \retval int - nb of found faces
2015 //================================================================================
2017 int SMDS_VolumeTool::GetAllExistingEdges(std::vector<const SMDS_MeshElement*> & edges) const
2020 edges.reserve( myVolumeNodes.size() * 2 );
2021 for ( size_t i = 0; i < myVolumeNodes.size()-1; ++i ) {
2022 for ( size_t j = i + 1; j < myVolumeNodes.size(); ++j ) {
2023 if ( IsLinked( i, j )) {
2024 const SMDS_MeshElement* edge =
2025 SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
2027 edges.push_back( edge );
2031 return edges.size();
2034 //================================================================================
2036 * \brief Return minimal square distance between connected corner nodes
2038 //================================================================================
2040 double SMDS_VolumeTool::MinLinearSize2() const
2042 double minSize = 1e+100;
2043 int iQ = myVolume->IsQuadratic() ? 2 : 1;
2045 SaveFacet savedFacet( myCurFace );
2047 // it seems that compute distance twice is faster than organization of a sole computing
2048 myCurFace.myIndex = -1;
2049 for ( int iF = 0; iF < myNbFaces; ++iF )
2052 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
2054 XYZ n1( myCurFace.myNodes[ iN ]);
2055 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
2056 minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
2063 //================================================================================
2065 * \brief Return maximal square distance between connected corner nodes
2067 //================================================================================
2069 double SMDS_VolumeTool::MaxLinearSize2() const
2071 double maxSize = -1e+100;
2072 int iQ = myVolume->IsQuadratic() ? 2 : 1;
2074 SaveFacet savedFacet( myCurFace );
2076 // it seems that compute distance twice is faster than organization of a sole computing
2077 myCurFace.myIndex = -1;
2078 for ( int iF = 0; iF < myNbFaces; ++iF )
2081 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
2083 XYZ n1( myCurFace.myNodes[ iN ]);
2084 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
2085 maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude());
2092 //================================================================================
2094 * \brief Fast quickly check that only one volume is built on the face nodes
2095 * This check is valid for conformal meshes only
2097 //================================================================================
2099 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
2101 const bool isFree = true;
2103 if ( !setFace( faceIndex ))
2106 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
2108 const int di = myVolume->IsQuadratic() ? 2 : 1;
2109 const int nbN = ( myCurFace.myNbNodes/di <= 4 && !IsPoly()) ? 3 : myCurFace.myNbNodes/di; // nb nodes to check
2111 SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
2112 while ( eIt->more() )
2114 const SMDS_MeshElement* vol = eIt->next();
2115 if ( vol == myVolume )
2118 for ( iN = 1; iN < nbN; ++iN )
2119 if ( vol->GetNodeIndex( nodes[ iN*di ]) < 0 )
2121 if ( iN == nbN ) // nbN nodes are shared with vol
2123 // if ( vol->IsPoly() || vol->NbFaces() > 6 ) // vol is polyhed or hex prism
2125 // int nb = myCurFace.myNbNodes;
2126 // if ( myVolume->GetEntityType() != vol->GetEntityType() )
2127 // nb -= ( GetCenterNodeIndex(0) > 0 );
2128 // std::set<const SMDS_MeshNode*> faceNodes( nodes, nodes + nb );
2129 // if ( SMDS_VolumeTool( vol ).GetFaceIndex( faceNodes ) < 0 )
2132 if ( otherVol ) *otherVol = vol;
2136 if ( otherVol ) *otherVol = 0;
2140 //================================================================================
2142 * \brief Check that only one volume is built on the face nodes
2143 * Different to IsFreeFace function, all nodes of the face are checked.
2144 * For non conforming meshes, the face that is not conform with the neighbor
2145 * will be identify as free.
2147 //================================================================================
2149 bool SMDS_VolumeTool::IsFreeFaceCheckAllNodes( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
2151 const bool isFree = true;
2153 if ( !setFace( faceIndex ))
2156 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
2158 const int di = myVolume->IsQuadratic() ? 2 : 1;
2159 const int nbN = myCurFace.myNbNodes/di;
2160 std::vector<bool> allNodesCoincideWithNeighbor(nbN,false);
2162 for (int nodeId = 0; nodeId < nbN; nodeId++)
2164 SMDS_ElemIteratorPtr eIt = nodes[nodeId]->GetInverseElementIterator( SMDSAbs_Volume );
2166 while ( eIt->more() )
2168 const SMDS_MeshElement* vol = eIt->next();
2169 if ( vol == myVolume )
2177 if ( count==0 /*free corner in the face means free face*/)
2179 if ( otherVol ) *otherVol = 0;
2183 return IsFreeFace( faceIndex, otherVol );
2186 //================================================================================
2188 * \brief Thorough check that only one volume is built on the face nodes
2190 //================================================================================
2192 bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
2194 const bool isFree = true;
2196 if (!setFace( faceIndex ))
2199 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
2200 const int nbFaceNodes = myCurFace.myNbNodes;
2202 // evaluate nb of face nodes shared by other volumes
2203 int maxNbShared = -1;
2204 typedef std::map< const SMDS_MeshElement*, int > TElemIntMap;
2205 TElemIntMap volNbShared;
2206 TElemIntMap::iterator vNbIt;
2207 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
2208 const SMDS_MeshNode* n = nodes[ iNode ];
2209 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
2210 while ( eIt->more() ) {
2211 const SMDS_MeshElement* elem = eIt->next();
2212 if ( elem != myVolume ) {
2213 vNbIt = volNbShared.insert( std::make_pair( elem, 0 )).first;
2215 if ( vNbIt->second > maxNbShared )
2216 maxNbShared = vNbIt->second;
2220 if ( maxNbShared < 3 )
2221 return isFree; // is free
2223 // find volumes laying on the opposite side of the face
2224 // and sharing all nodes
2225 XYZ intNormal; // internal normal
2226 GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
2227 if ( IsFaceExternal( faceIndex ))
2228 intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
2229 XYZ p0 ( nodes[0] ), baryCenter;
2230 for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); ) {
2231 const int& nbShared = (*vNbIt).second;
2232 if ( nbShared >= 3 ) {
2233 SMDS_VolumeTool volume( (*vNbIt).first );
2234 volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
2235 XYZ intNormal2( baryCenter - p0 );
2236 if ( intNormal.Dot( intNormal2 ) < 0 ) {
2238 if ( nbShared >= nbFaceNodes )
2240 // a volume shares the whole facet
2241 if ( otherVol ) *otherVol = vNbIt->first;
2248 // remove a volume from volNbShared map
2249 volNbShared.erase( vNbIt++ );
2252 // here volNbShared contains only volumes laying on the opposite side of
2253 // the face and sharing 3 or more but not all face nodes with myVolume
2254 if ( volNbShared.size() < 2 ) {
2255 return isFree; // is free
2258 // check if the whole area of a face is shared
2259 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
2261 const SMDS_MeshNode* n = nodes[ iNode ];
2262 // check if n is shared by one of volumes of volNbShared
2263 bool isShared = false;
2264 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
2265 while ( eIt->more() && !isShared )
2266 isShared = volNbShared.count( eIt->next() );
2270 if ( otherVol ) *otherVol = volNbShared.begin()->first;
2273 // if ( !myVolume->IsPoly() )
2275 // bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
2276 // for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
2277 // SMDS_VolumeTool volume( (*vNbIt).first );
2278 // bool prevLinkShared = false;
2279 // int nbSharedLinks = 0;
2280 // for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
2281 // bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
2282 // if ( linkShared )
2284 // if ( linkShared && prevLinkShared &&
2285 // volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
2286 // isShared[ iNode ] = true;
2287 // prevLinkShared = linkShared;
2289 // if ( nbSharedLinks == nbFaceNodes )
2290 // return !free; // is not free
2291 // if ( nbFaceNodes == 4 ) {
2292 // // check traingle parts 1 & 3
2293 // if ( isShared[1] && isShared[3] )
2294 // return !free; // is not free
2295 // // check triangle parts 0 & 2;
2296 // // 0 part could not be checked in the loop; check it here
2297 // if ( isShared[2] && prevLinkShared &&
2298 // volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
2299 // volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
2300 // return !free; // is not free
2307 //=======================================================================
2308 //function : GetFaceIndex
2309 //purpose : Return index of a face formed by theFaceNodes
2310 //=======================================================================
2312 int SMDS_VolumeTool::GetFaceIndex( const std::set<const SMDS_MeshNode*>& theFaceNodes,
2313 const int theFaceIndexHint ) const
2315 if ( theFaceIndexHint >= 0 )
2317 int nbNodes = NbFaceNodes( theFaceIndexHint );
2318 if ( nbNodes == (int) theFaceNodes.size() )
2320 const SMDS_MeshNode** nodes = GetFaceNodes( theFaceIndexHint );
2322 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
2327 return theFaceIndexHint;
2330 for ( int iFace = 0; iFace < myNbFaces; iFace++ )
2332 if ( iFace == theFaceIndexHint )
2334 int nbNodes = NbFaceNodes( iFace );
2335 if ( nbNodes == (int) theFaceNodes.size() )
2337 const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
2339 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
2350 //=======================================================================
2351 //function : GetFaceIndex
2352 //purpose : Return index of a face formed by theFaceNodes
2353 //=======================================================================
2355 /*int SMDS_VolumeTool::GetFaceIndex( const std::set<int>& theFaceNodesIndices )
2357 for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
2358 const int* nodes = GetFaceNodesIndices( iFace );
2359 int nbFaceNodes = NbFaceNodes( iFace );
2360 std::set<int> nodeSet;
2361 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
2362 nodeSet.insert( nodes[ iNode ] );
2363 if ( theFaceNodesIndices == nodeSet )
2369 //=======================================================================
2370 //function : setFace
2372 //=======================================================================
2374 bool SMDS_VolumeTool::setFace( int faceIndex ) const
2379 if ( myCurFace.myIndex == faceIndex )
2382 myCurFace.myIndex = -1;
2384 if ( faceIndex < 0 || faceIndex >= NbFaces() )
2387 if (myVolume->IsPoly())
2389 if ( !myPolyedre ) {
2390 MESSAGE("Warning: bad volumic element");
2395 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
2396 if ( !myAllFacesNbNodes ) {
2397 me->myPolyQuantities = myPolyedre->GetQuantities();
2398 myAllFacesNbNodes = &myPolyQuantities[0];
2400 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2401 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2402 me->myPolyIndices.resize( myCurFace.myNbNodes + 1 );
2403 myCurFace.myNodeIndices = & me->myPolyIndices[0];
2404 int shift = std::accumulate( myAllFacesNbNodes, myAllFacesNbNodes+faceIndex, 0 );
2405 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2407 myCurFace.myNodes [ iNode ] = myVolumeNodes[ shift + iNode ];
2408 myCurFace.myNodeIndices[ iNode ] = shift + iNode;
2410 myCurFace.myNodes [ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ]; // last = first
2411 myCurFace.myNodeIndices[ myCurFace.myNbNodes ] = myCurFace.myNodeIndices[ 0 ];
2413 // check orientation
2414 if (myExternalFaces)
2416 myCurFace.myIndex = faceIndex; // avoid infinite recursion in IsFaceExternal()
2417 myExternalFaces = false; // force normal computation by IsFaceExternal()
2418 if ( !IsFaceExternal( faceIndex ))
2419 std::reverse( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
2420 myExternalFaces = true;
2425 if ( !myAllFacesNodeIndices_F )
2427 // choose data for an element type
2428 switch ( myVolumeNodes.size() ) {
2430 myAllFacesNodeIndices_F = &Tetra_F [0][0];
2431 //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
2432 myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
2433 myAllFacesNbNodes = Tetra_nbN;
2434 myMaxFaceNbNodes = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
2437 myAllFacesNodeIndices_F = &Pyramid_F [0][0];
2438 //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
2439 myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
2440 myAllFacesNbNodes = Pyramid_nbN;
2441 myMaxFaceNbNodes = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
2444 myAllFacesNodeIndices_F = &Penta_F [0][0];
2445 //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
2446 myAllFacesNodeIndices_RE = &Penta_RE[0][0];
2447 myAllFacesNbNodes = Penta_nbN;
2448 myMaxFaceNbNodes = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
2451 myAllFacesNodeIndices_F = &Hexa_F [0][0];
2452 ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
2453 myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
2454 myAllFacesNbNodes = Hexa_nbN;
2455 myMaxFaceNbNodes = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
2458 myAllFacesNodeIndices_F = &QuadTetra_F [0][0];
2459 //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
2460 myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
2461 myAllFacesNbNodes = QuadTetra_nbN;
2462 myMaxFaceNbNodes = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
2465 myAllFacesNodeIndices_F = &QuadPyram_F [0][0];
2466 //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
2467 myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
2468 myAllFacesNbNodes = QuadPyram_nbN;
2469 myMaxFaceNbNodes = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
2472 myAllFacesNodeIndices_F = &QuadPenta_F [0][0];
2473 //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
2474 myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
2475 myAllFacesNbNodes = QuadPenta_nbN;
2476 myMaxFaceNbNodes = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
2480 myAllFacesNodeIndices_F = &QuadHexa_F [0][0];
2481 //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
2482 myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
2483 myAllFacesNbNodes = QuadHexa_nbN;
2484 myMaxFaceNbNodes = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
2485 if ( !myIgnoreCentralNodes && myVolumeNodes.size() == 27 )
2487 myAllFacesNodeIndices_F = &TriQuadHexa_F [0][0];
2488 //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
2489 myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
2490 myAllFacesNbNodes = TriQuadHexa_nbN;
2491 myMaxFaceNbNodes = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
2495 myAllFacesNodeIndices_F = &HexPrism_F [0][0];
2496 //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
2497 myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
2498 myAllFacesNbNodes = HexPrism_nbN;
2499 myMaxFaceNbNodes = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
2505 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2506 // if ( myExternalFaces )
2507 // myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2509 // myCurFace.myNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
2510 myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2513 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2514 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2515 myCurFace.myNodes[ iNode ] = myVolumeNodes[ myCurFace.myNodeIndices[ iNode ]];
2516 myCurFace.myNodes[ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ];
2519 myCurFace.myIndex = faceIndex;
2524 //=======================================================================
2525 //function : GetType
2526 //purpose : return VolumeType by nb of nodes in a volume
2527 //=======================================================================
2529 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
2531 switch ( nbNodes ) {
2532 case 4: return TETRA;
2533 case 5: return PYRAM;
2534 case 6: return PENTA;
2535 case 8: return HEXA;
2536 case 10: return QUAD_TETRA;
2537 case 13: return QUAD_PYRAM;
2538 case 15: return QUAD_PENTA;
2540 case 27: return QUAD_HEXA;
2541 case 12: return HEX_PRISM;
2542 default:return UNKNOWN;
2546 //=======================================================================
2547 //function : NbFaces
2548 //purpose : return nb of faces by volume type
2549 //=======================================================================
2551 int SMDS_VolumeTool::NbFaces( VolumeType type )
2555 case QUAD_TETRA: return 4;
2557 case QUAD_PYRAM: return 5;
2559 case QUAD_PENTA: return 5;
2561 case QUAD_HEXA : return 6;
2562 case HEX_PRISM : return 8;
2567 //================================================================================
2569 * \brief Useful to know nb of corner nodes of a quadratic volume
2570 * \param type - volume type
2571 * \retval int - nb of corner nodes
2573 //================================================================================
2575 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
2579 case QUAD_TETRA: return 4;
2581 case QUAD_PYRAM: return 5;
2583 case QUAD_PENTA: return 6;
2585 case QUAD_HEXA : return 8;
2586 case HEX_PRISM : return 12;
2593 //=======================================================================
2594 //function : GetFaceNodesIndices
2595 //purpose : Return the array of face nodes indices
2596 // To comfort link iteration, the array
2597 // length == NbFaceNodes( faceIndex ) + 1 and
2598 // the last node index == the first one.
2599 //=======================================================================
2601 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
2606 case TETRA: return Tetra_F[ faceIndex ];
2607 case PYRAM: return Pyramid_F[ faceIndex ];
2608 case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
2609 case HEXA: return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
2610 case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
2611 case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
2612 case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
2613 // what about SMDSEntity_TriQuad_Hexa?
2614 case QUAD_HEXA: return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
2615 case HEX_PRISM: return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
2621 //=======================================================================
2622 //function : NbFaceNodes
2623 //purpose : Return number of nodes in the array of face nodes
2624 //=======================================================================
2626 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
2630 case TETRA: return Tetra_nbN[ faceIndex ];
2631 case PYRAM: return Pyramid_nbN[ faceIndex ];
2632 case PENTA: return Penta_nbN[ faceIndex ];
2633 case HEXA: return Hexa_nbN[ faceIndex ];
2634 case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
2635 case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
2636 case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
2637 // what about SMDSEntity_TriQuad_Hexa?
2638 case QUAD_HEXA: return QuadHexa_nbN[ faceIndex ];
2639 case HEX_PRISM: return HexPrism_nbN[ faceIndex ];
2645 //=======================================================================
2646 //function : Element
2647 //purpose : return element
2648 //=======================================================================
2650 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
2652 return static_cast<const SMDS_MeshVolume*>( myVolume );
2655 //=======================================================================
2657 //purpose : return element ID
2658 //=======================================================================
2660 smIdType SMDS_VolumeTool::ID() const
2662 return myVolume ? myVolume->GetID() : 0;