1 // Copyright (C) 2007-2023 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>
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 };
94 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
95 { 0, 1, 2, 3, 0 }, // All faces have external normals
101 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
102 { 0, 3, 2, 1, 0 }, // All faces but a bottom have external normals
107 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
118 // | / \ | PENTAHEDRON
124 static int Penta_F [5][5] = { // FORWARD
125 { 0, 1, 2, 0, 0 }, // All faces have external normals
126 { 3, 5, 4, 3, 3 }, // 0 is bottom, 1 is top face
130 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
136 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
143 // N4+----------+N7 |
144 // | | | | HEXAHEDRON
145 // | N1+------|---+N2
151 static int Hexa_F [6][5] = { // FORWARD
153 { 4, 7, 6, 5, 4 }, // all face normals are external
158 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
160 { 4, 5, 6, 7, 4 }, // all face normals are external
165 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
166 static int Hexa_oppF[] = { 1, 0, 4, 5, 2, 3 }; // oppopsite facet indices
185 static int HexPrism_F [8][7] = { // FORWARD
186 { 0, 1, 2, 3, 4, 5, 0 },
187 { 6,11,10, 9, 8, 7, 6 },
188 { 0, 6, 7, 1, 0, 0, 0 },
189 { 1, 7, 8, 2, 1, 1, 1 },
190 { 2, 8, 9, 3, 2, 2, 2 },
191 { 3, 9,10, 4, 3, 3, 3 },
192 { 4,10,11, 5, 4, 4, 4 },
193 { 5,11, 6, 0, 5, 5, 5 }};
194 static int HexPrism_RE [8][7] = { // REVERSED -> EXTERNAL
195 { 0, 5, 4, 3, 2, 1, 0 },
196 { 6,11,10, 9, 8, 7, 6 },
197 { 0, 6, 7, 1, 0, 0, 0 },
198 { 1, 7, 8, 2, 1, 1, 1 },
199 { 2, 8, 9, 3, 2, 2, 2 },
200 { 3, 9,10, 4, 3, 3, 3 },
201 { 4,10,11, 5, 4, 4, 4 },
202 { 5,11, 6, 0, 5, 5, 5 }};
203 static int HexPrism_nbN [] = { 6, 6, 4, 4, 4, 4, 4, 4 };
212 // N0 +---|---+ N1 TETRAHEDRON
220 static int QuadTetra_F [4][7] = { // FORWARD
221 { 0, 4, 1, 5, 2, 6, 0 }, // All faces have external normals
222 { 0, 7, 3, 8, 1, 4, 0 },
223 { 1, 8, 3, 9, 2, 5, 1 },
224 { 0, 6, 2, 9, 3, 7, 0 }};
225 static int QuadTetra_RE [4][7] = { // REVERSED -> FORWARD (EXTERNAL)
226 { 0, 6, 2, 5, 1, 4, 0 }, // All faces have external normals
227 { 0, 4, 1, 8, 3, 7, 0 },
228 { 1, 5, 2, 9, 3, 8, 1 },
229 { 0, 7, 3, 9, 2, 6, 0 }};
230 static int QuadTetra_nbN [] = { 6, 6, 6, 6 };
240 // | | 9 - middle point for (0,4) etc.
253 static int QuadPyram_F [5][9] = { // FORWARD
254 { 0, 5, 1, 6, 2, 7, 3, 8, 0 }, // All faces have external normals
255 { 0, 9, 4, 10,1, 5, 0, 4, 4 },
256 { 1, 10,4, 11,2, 6, 1, 4, 4 },
257 { 2, 11,4, 12,3, 7, 2, 4, 4 },
258 { 3, 12,4, 9, 0, 8, 3, 4, 4 }};
259 static int QuadPyram_RE [5][9] = { // REVERSED -> FORWARD (EXTERNAL)
260 { 0, 8, 3, 7, 2, 6, 1, 5, 0 }, // All faces but a bottom have external normals
261 { 0, 5, 1, 10,4, 9, 0, 4, 4 },
262 { 1, 6, 2, 11,4, 10,1, 4, 4 },
263 { 2, 7, 3, 12,4, 11,2, 4, 4 },
264 { 3, 8, 0, 9, 4, 12,3, 4, 4 }};
265 static int QuadPyram_nbN [] = { 8, 6, 6, 6, 6 };
289 static int QuadPenta_F [5][9] = { // FORWARD
290 { 0, 6, 1, 7, 2, 8, 0, 0, 0 },
291 { 3, 11,5, 10,4, 9, 3, 3, 3 },
292 { 0, 12,3, 9, 4, 13,1, 6, 0 },
293 { 1, 13,4, 10,5, 14,2, 7, 1 },
294 { 0, 8, 2, 14,5, 11,3, 12,0 }};
295 static int QuadPenta_RE [5][9] = { // REVERSED -> EXTERNAL
296 { 0, 8, 2, 7, 1, 6, 0, 0, 0 },
297 { 3, 9, 4, 10,5, 11,3, 3, 3 },
298 { 0, 6, 1, 13,4, 9, 3, 12,0 },
299 { 1, 7, 2, 14,5, 10,4, 13,1 },
300 { 0, 12,3, 11,5, 14,2, 8, 0 }};
301 static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 };
305 // N5+-----+-----+N6 +-----+-----+
307 // 12+ | 14+ | + | +25 + |
309 // N4+-----+-----+N7 | QUADRATIC +-----+-----+ | Central nodes
310 // | | 15 | | HEXAHEDRON | | | | of tri-quadratic
311 // | | | | | | | | HEXAHEDRON
312 // | 17+ | +18 | + 22+ | +
314 // | | | | | + | 26+ | + |
316 // 16+ | +19 | + | +24 + |
319 // | N1+-----+-|---+N2 | +-----+-|---+
321 // | +8 | +10 | + 20+ | +
323 // N0+-----+-----+N3 +-----+-----+
326 static int QuadHexa_F [6][9] = { // FORWARD
327 { 0, 8, 1, 9, 2, 10,3, 11,0 }, // all face normals are external,
328 { 4, 15,7, 14,6, 13,5, 12,4 },
329 { 0, 16,4, 12,5, 17,1, 8, 0 },
330 { 1, 17,5, 13,6, 18,2, 9, 1 },
331 { 3, 10,2, 18,6, 14,7, 19,3 },
332 { 0, 11,3, 19,7, 15,4, 16,0 }};
333 static int QuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL
334 { 0, 11,3, 10,2, 9, 1, 8, 0 }, // all face normals are external
335 { 4, 12,5, 13,6, 14,7, 15,4 },
336 { 0, 8, 1, 17,5, 12,4, 16,0 },
337 { 1, 9, 2, 18,6, 13,5, 17,1 },
338 { 3, 19,7, 14,6, 18,2, 10,3 },
339 { 0, 16,4, 15,7, 19,3, 11,0 }};
340 static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 };
342 static int TriQuadHexa_F [6][9] = { // FORWARD
343 { 0, 8, 1, 9, 2, 10,3, 11, 20 }, // all face normals are external
344 { 4, 15,7, 14,6, 13,5, 12, 25 },
345 { 0, 16,4, 12,5, 17,1, 8, 21 },
346 { 1, 17,5, 13,6, 18,2, 9, 22 },
347 { 3, 10,2, 18,6, 14,7, 19, 23 },
348 { 0, 11,3, 19,7, 15,4, 16, 24 }};
349 static int TriQuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL
350 { 0, 11,3, 10,2, 9, 1, 8, 20 }, // opposite faces are neighbouring,
351 { 4, 12,5, 13,6, 14,7, 15, 25 }, // all face normals are external
352 { 0, 8, 1, 17,5, 12,4, 16, 21 },
353 { 1, 9, 2, 18,6, 13,5, 17, 22 },
354 { 3, 19,7, 14,6, 18,2, 10, 23 },
355 { 0, 16,4, 15,7, 19,3, 11, 24 }};
356 static int TriQuadHexa_nbN [] = { 9, 9, 9, 9, 9, 9 };
359 // ========================================================
360 // to perform some calculations without linkage to CASCADE
361 // ========================================================
366 XYZ() { x = 0; y = 0; z = 0; }
367 XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
368 XYZ( const XYZ& other ) { x = other.x; y = other.y; z = other.z; }
369 XYZ( const SMDS_MeshNode* n ) { x = n->X(); y = n->Y(); z = n->Z(); }
370 double* data() { return &x; }
371 inline XYZ operator-( const XYZ& other );
372 inline XYZ operator+( const XYZ& other );
373 inline XYZ Crossed( const XYZ& other );
374 inline XYZ operator-();
375 inline double Dot( const XYZ& other );
376 inline double Magnitude();
377 inline double SquareMagnitude();
378 inline XYZ Normalize();
380 inline XYZ XYZ::operator-( const XYZ& Right ) {
381 return XYZ(x - Right.x, y - Right.y, z - Right.z);
383 inline XYZ XYZ::operator-() {
384 return XYZ(-x,-y,-z);
386 inline XYZ XYZ::operator+( const XYZ& Right ) {
387 return XYZ(x + Right.x, y + Right.y, z + Right.z);
389 inline XYZ XYZ::Crossed( const XYZ& Right ) {
390 return XYZ (y * Right.z - z * Right.y,
391 z * Right.x - x * Right.z,
392 x * Right.y - y * Right.x);
394 inline double XYZ::Dot( const XYZ& Other ) {
395 return(x * Other.x + y * Other.y + z * Other.z);
397 inline double XYZ::Magnitude() {
398 return sqrt (x * x + y * y + z * z);
400 inline double XYZ::SquareMagnitude() {
401 return (x * x + y * y + z * z);
403 inline XYZ XYZ::Normalize() {
404 double magnitude = Magnitude();
405 if ( magnitude != 0.0 )
406 return XYZ(x /= magnitude,y /= magnitude,z /= magnitude );
411 //================================================================================
413 * \brief Return linear type corresponding to a quadratic one
415 //================================================================================
417 SMDS_VolumeTool::VolumeType quadToLinear(SMDS_VolumeTool::VolumeType quadType)
419 SMDS_VolumeTool::VolumeType linType = SMDS_VolumeTool::VolumeType( int(quadType)-4 );
420 const int nbCornersByQuad = SMDS_VolumeTool::NbCornerNodes( quadType );
421 if ( SMDS_VolumeTool::NbCornerNodes( linType ) == nbCornersByQuad )
425 for ( ; iLin < SMDS_VolumeTool::NB_VOLUME_TYPES; ++iLin )
426 if ( SMDS_VolumeTool::NbCornerNodes( SMDS_VolumeTool::VolumeType( iLin )) == nbCornersByQuad)
427 return SMDS_VolumeTool::VolumeType( iLin );
429 return SMDS_VolumeTool::UNKNOWN;
434 //================================================================================
436 * \brief Saver/restorer of a SMDS_VolumeTool::myCurFace
438 //================================================================================
440 struct SMDS_VolumeTool::SaveFacet
442 SMDS_VolumeTool::Facet mySaved;
443 SMDS_VolumeTool::Facet& myToRestore;
444 SaveFacet( SMDS_VolumeTool::Facet& facet ): myToRestore( facet )
447 mySaved.myNodes.swap( facet.myNodes );
451 if ( myToRestore.myIndex != mySaved.myIndex )
452 myToRestore = mySaved;
453 myToRestore.myNodes.swap( mySaved.myNodes );
457 //=======================================================================
458 //function : SMDS_VolumeTool
460 //=======================================================================
462 SMDS_VolumeTool::SMDS_VolumeTool ()
467 //=======================================================================
468 //function : SMDS_VolumeTool
470 //=======================================================================
472 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume,
473 const bool ignoreCentralNodes)
475 Set( theVolume, ignoreCentralNodes );
478 //=======================================================================
479 //function : SMDS_VolumeTool
481 //=======================================================================
483 SMDS_VolumeTool::~SMDS_VolumeTool()
485 myCurFace.myNodeIndices = NULL;
488 //=======================================================================
489 //function : SetVolume
490 //purpose : Set volume to iterate on
491 //=======================================================================
493 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume,
494 const bool ignoreCentralNodes,
495 const std::vector<const SMDS_MeshNode*>* otherNodes)
500 myIgnoreCentralNodes = ignoreCentralNodes;
504 myVolumeNodes.clear();
505 myPolyIndices.clear();
506 myPolyQuantities.clear();
507 myPolyFacetOri.clear();
510 myExternalFaces = false;
512 myAllFacesNodeIndices_F = 0;
513 myAllFacesNodeIndices_RE = 0;
514 myAllFacesNbNodes = 0;
516 myCurFace.myIndex = -1;
517 myCurFace.myNodeIndices = NULL;
518 myCurFace.myNodes.clear();
521 if ( !theVolume || theVolume->GetType() != SMDSAbs_Volume )
524 myVolume = theVolume;
525 myNbFaces = theVolume->NbFaces();
526 if ( myVolume->IsPoly() )
528 myPolyedre = SMDS_Mesh::DownCast<SMDS_MeshVolume>( myVolume );
529 myPolyFacetOri.resize( myNbFaces, 0 );
533 myVolumeNodes.resize( myVolume->NbNodes() );
536 if ( otherNodes->size() != myVolumeNodes.size() )
537 return ( myVolume = 0 );
538 for ( size_t i = 0; i < otherNodes->size(); ++i )
539 if ( ! ( myVolumeNodes[i] = (*otherNodes)[0] ))
540 return ( myVolume = 0 );
544 myVolumeNodes.assign( myVolume->begin_nodes(), myVolume->end_nodes() );
549 return ( myVolume = 0 );
553 // define volume orientation
555 if ( GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z ))
557 const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
558 int topNodeIndex = myVolume->NbCornerNodes() - 1;
559 while ( !IsLinked( 0, topNodeIndex, /*ignoreMediumNodes=*/true )) --topNodeIndex;
560 const SMDS_MeshNode* topNode = myVolumeNodes[ topNodeIndex ];
561 XYZ upDir (topNode->X() - botNode->X(),
562 topNode->Y() - botNode->Y(),
563 topNode->Z() - botNode->Z() );
564 myVolForward = ( botNormal.Dot( upDir ) < 0 );
567 myCurFace.myIndex = -1; // previous setFace(0) didn't take myVolForward into account
572 //=======================================================================
574 //purpose : Inverse volume
575 //=======================================================================
577 #define SWAP_NODES(nodes,i1,i2) \
579 const SMDS_MeshNode* tmp = nodes[ i1 ]; \
580 nodes[ i1 ] = nodes[ i2 ]; \
583 void SMDS_VolumeTool::Inverse ()
585 if ( !myVolume ) return;
587 if (myVolume->IsPoly()) {
588 MESSAGE("Warning: attempt to inverse polyhedral volume");
592 myVolForward = !myVolForward;
593 myCurFace.myIndex = -1;
595 // inverse top and bottom faces
596 switch ( myVolumeNodes.size() ) {
598 SWAP_NODES( myVolumeNodes, 1, 2 );
601 SWAP_NODES( myVolumeNodes, 1, 3 );
604 SWAP_NODES( myVolumeNodes, 1, 2 );
605 SWAP_NODES( myVolumeNodes, 4, 5 );
608 SWAP_NODES( myVolumeNodes, 1, 3 );
609 SWAP_NODES( myVolumeNodes, 5, 7 );
612 SWAP_NODES( myVolumeNodes, 1, 5 );
613 SWAP_NODES( myVolumeNodes, 2, 4 );
614 SWAP_NODES( myVolumeNodes, 7, 11 );
615 SWAP_NODES( myVolumeNodes, 8, 10 );
619 SWAP_NODES( myVolumeNodes, 1, 2 );
620 SWAP_NODES( myVolumeNodes, 4, 6 );
621 SWAP_NODES( myVolumeNodes, 8, 9 );
624 SWAP_NODES( myVolumeNodes, 1, 3 );
625 SWAP_NODES( myVolumeNodes, 5, 8 );
626 SWAP_NODES( myVolumeNodes, 6, 7 );
627 SWAP_NODES( myVolumeNodes, 10, 12 );
630 SWAP_NODES( myVolumeNodes, 1, 2 );
631 SWAP_NODES( myVolumeNodes, 4, 5 );
632 SWAP_NODES( myVolumeNodes, 6, 8 );
633 SWAP_NODES( myVolumeNodes, 9, 11 );
634 SWAP_NODES( myVolumeNodes, 13, 14 );
637 SWAP_NODES( myVolumeNodes, 1, 3 );
638 SWAP_NODES( myVolumeNodes, 5, 7 );
639 SWAP_NODES( myVolumeNodes, 8, 11 );
640 SWAP_NODES( myVolumeNodes, 9, 10 );
641 SWAP_NODES( myVolumeNodes, 12, 15 );
642 SWAP_NODES( myVolumeNodes, 13, 14 );
643 SWAP_NODES( myVolumeNodes, 17, 19 );
646 SWAP_NODES( myVolumeNodes, 1, 3 );
647 SWAP_NODES( myVolumeNodes, 5, 7 );
648 SWAP_NODES( myVolumeNodes, 8, 11 );
649 SWAP_NODES( myVolumeNodes, 9, 10 );
650 SWAP_NODES( myVolumeNodes, 12, 15 );
651 SWAP_NODES( myVolumeNodes, 13, 14 );
652 SWAP_NODES( myVolumeNodes, 17, 19 );
653 SWAP_NODES( myVolumeNodes, 21, 24 );
654 SWAP_NODES( myVolumeNodes, 22, 23 );
660 //=======================================================================
661 //function : GetVolumeType
663 //=======================================================================
665 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
670 switch( myVolumeNodes.size() ) {
671 case 4: return TETRA;
672 case 5: return PYRAM;
673 case 6: return PENTA;
675 case 12: return HEX_PRISM;
676 case 10: return QUAD_TETRA;
677 case 13: return QUAD_PYRAM;
678 case 15: return QUAD_PENTA;
679 case 20: return QUAD_HEXA;
680 case 27: return QUAD_HEXA;
687 //=======================================================================
688 //function : getTetraVolume
690 //=======================================================================
692 static double getTetraVolume(const SMDS_MeshNode* n1,
693 const SMDS_MeshNode* n2,
694 const SMDS_MeshNode* n3,
695 const SMDS_MeshNode* n4)
697 double p1[3], p2[3], p3[3], p4[3];
703 double Q1 = -(p1[ 0 ]-p2[ 0 ])*(p3[ 1 ]*p4[ 2 ]-p4[ 1 ]*p3[ 2 ]);
704 double Q2 = (p1[ 0 ]-p3[ 0 ])*(p2[ 1 ]*p4[ 2 ]-p4[ 1 ]*p2[ 2 ]);
705 double R1 = -(p1[ 0 ]-p4[ 0 ])*(p2[ 1 ]*p3[ 2 ]-p3[ 1 ]*p2[ 2 ]);
706 double R2 = -(p2[ 0 ]-p3[ 0 ])*(p1[ 1 ]*p4[ 2 ]-p4[ 1 ]*p1[ 2 ]);
707 double S1 = (p2[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p3[ 2 ]-p3[ 1 ]*p1[ 2 ]);
708 double S2 = -(p3[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p2[ 2 ]-p2[ 1 ]*p1[ 2 ]);
710 return (Q1+Q2+R1+R2+S1+S2)/6.0;
713 //=======================================================================
715 //purpose : Return element volume
716 //=======================================================================
717 double SMDS_VolumeTool::GetSize() const
723 if ( myVolume->IsPoly() )
728 SaveFacet savedFacet( myCurFace );
730 // split a polyhedron into tetrahedrons
733 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
734 for ( int f = 0; f < NbFaces(); ++f )
737 XYZ area (0,0,0), p1( myCurFace.myNodes[0] );
738 for ( int n = 0; n < myCurFace.myNbNodes; ++n )
740 XYZ p2( myCurFace.myNodes[ n+1 ]);
741 area = area + p1.Crossed( p2 );
745 oriOk = oriOk && IsFaceExternal( f );
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 //=======================================================================
1599 //function : GetFaceArea
1600 //purpose : Return face area
1601 //=======================================================================
1603 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
1606 if ( !setFace( faceIndex ))
1609 XYZ p1 ( myCurFace.myNodes[0] );
1610 XYZ p2 ( myCurFace.myNodes[1] );
1611 XYZ p3 ( myCurFace.myNodes[2] );
1612 XYZ aVec12( p2 - p1 );
1613 XYZ aVec13( p3 - p1 );
1614 area += aVec12.Crossed( aVec13 ).Magnitude();
1616 if (myVolume->IsPoly())
1618 for ( int i = 3; i < myCurFace.myNbNodes; ++i )
1620 XYZ pI ( myCurFace.myNodes[i] );
1621 XYZ aVecI( pI - p1 );
1622 area += aVec13.Crossed( aVecI ).Magnitude();
1628 if ( myCurFace.myNbNodes == 4 ) {
1629 XYZ p4 ( myCurFace.myNodes[3] );
1630 XYZ aVec14( p4 - p1 );
1631 area += aVec14.Crossed( aVec13 ).Magnitude();
1637 //================================================================================
1639 * \brief Return index of the node located at face center of a quadratic element like HEX27
1641 //================================================================================
1643 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
1645 if ( myAllFacesNbNodes && myVolumeNodes.size() == 27 ) // classic element with 27 nodes
1647 switch ( faceIndex ) {
1651 return faceIndex + 19;
1657 //=======================================================================
1658 //function : GetOppFaceIndex
1659 //purpose : Return index of the opposite face if it exists, else -1.
1660 //=======================================================================
1662 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1666 MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1670 const int nbHoriFaces = 2;
1672 if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1673 switch ( myVolumeNodes.size() ) {
1676 if ( faceIndex == 0 || faceIndex == 1 )
1677 ind = 1 - faceIndex;
1681 if ( faceIndex <= 1 ) // top or bottom
1682 ind = 1 - faceIndex;
1684 const int nbSideFaces = myAllFacesNbNodes[0];
1685 ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1690 ind = GetOppFaceIndexOfHex( faceIndex );
1698 //=======================================================================
1699 //function : GetOppFaceIndexOfHex
1700 //purpose : Return index of the opposite face of the hexahedron
1701 //=======================================================================
1703 int SMDS_VolumeTool::GetOppFaceIndexOfHex( int faceIndex )
1705 return Hexa_oppF[ faceIndex ];
1708 //=======================================================================
1709 //function : IsLinked
1710 //purpose : return true if theNode1 is linked with theNode2
1711 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1712 //=======================================================================
1714 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1715 const SMDS_MeshNode* theNode2,
1716 const bool theIgnoreMediumNodes) const
1721 if (myVolume->IsPoly()) {
1723 MESSAGE("Warning: bad volumic element");
1726 if ( !myAllFacesNbNodes ) {
1727 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1728 me->myPolyQuantities = myPolyedre->GetQuantities();
1729 myAllFacesNbNodes = &myPolyQuantities[0];
1731 int from, to = 0, d1 = 1, d2 = 2;
1732 if ( myPolyedre->IsQuadratic() ) {
1733 if ( theIgnoreMediumNodes ) {
1739 std::vector<const SMDS_MeshNode*>::const_iterator i;
1740 for (int iface = 0; iface < myNbFaces; iface++)
1743 to += myPolyQuantities[iface];
1744 i = std::find( myVolumeNodes.begin() + from, myVolumeNodes.begin() + to, theNode1 );
1745 if ( i != myVolumeNodes.end() )
1747 if (( theNode2 == *( i-d1 ) ||
1748 theNode2 == *( i+d1 )))
1751 (( theNode2 == *( i-d2 ) ||
1752 theNode2 == *( i+d2 ))))
1759 // find nodes indices
1760 int i1 = -1, i2 = -1, nbFound = 0;
1761 for ( size_t i = 0; i < myVolumeNodes.size() && nbFound < 2; i++ )
1763 if ( myVolumeNodes[ i ] == theNode1 )
1765 else if ( myVolumeNodes[ i ] == theNode2 )
1768 return IsLinked( i1, i2 );
1771 //=======================================================================
1772 //function : IsLinked
1773 //purpose : return true if the node with theNode1Index is linked
1774 // with the node with theNode2Index
1775 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1776 //=======================================================================
1778 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1779 const int theNode2Index,
1780 bool theIgnoreMediumNodes) const
1782 if ( myVolume->IsPoly() ) {
1783 return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1786 int minInd = std::min( theNode1Index, theNode2Index );
1787 int maxInd = std::max( theNode1Index, theNode2Index );
1789 if ( minInd < 0 || maxInd > (int)myVolumeNodes.size() - 1 || maxInd == minInd )
1792 VolumeType type = GetVolumeType();
1793 if ( myVolume->IsQuadratic() )
1795 int firstMediumInd = myVolume->NbCornerNodes();
1796 if ( minInd >= firstMediumInd )
1797 return false; // both nodes are medium - not linked
1798 if ( maxInd < firstMediumInd ) // both nodes are corners
1800 if ( theIgnoreMediumNodes )
1801 type = quadToLinear(type); // to check linkage of corner nodes only
1803 return false; // corner nodes are not linked directly in a quadratic cell
1811 switch ( maxInd - minInd ) {
1812 case 1: return minInd != 3;
1813 case 3: return minInd == 0 || minInd == 4;
1814 case 4: return true;
1821 switch ( maxInd - minInd ) {
1823 case 3: return true;
1828 switch ( maxInd - minInd ) {
1829 case 1: return minInd != 2;
1830 case 2: return minInd == 0 || minInd == 3;
1831 case 3: return true;
1838 case 0: return ( maxInd==4 || maxInd==6 || maxInd==7 );
1839 case 1: return ( maxInd==4 || maxInd==5 || maxInd==8 );
1840 case 2: return ( maxInd==5 || maxInd==6 || maxInd==9 );
1841 case 3: return ( maxInd==7 || maxInd==8 || maxInd==9 );
1849 case 0: return ( maxInd==8 || maxInd==11 || maxInd==16 );
1850 case 1: return ( maxInd==8 || maxInd==9 || maxInd==17 );
1851 case 2: return ( maxInd==9 || maxInd==10 || maxInd==18 );
1852 case 3: return ( maxInd==10 || maxInd==11 || maxInd==19 );
1853 case 4: return ( maxInd==12 || maxInd==15 || maxInd==16 );
1854 case 5: return ( maxInd==12 || maxInd==13 || maxInd==17 );
1855 case 6: return ( maxInd==13 || maxInd==14 || maxInd==18 );
1856 case 7: return ( maxInd==14 || maxInd==15 || maxInd==19 );
1864 case 0: return ( maxInd==5 || maxInd==8 || maxInd==9 );
1865 case 1: return ( maxInd==5 || maxInd==6 || maxInd==10 );
1866 case 2: return ( maxInd==6 || maxInd==7 || maxInd==11 );
1867 case 3: return ( maxInd==7 || maxInd==8 || maxInd==12 );
1868 case 4: return ( maxInd==9 || maxInd==10 || maxInd==11 || maxInd==12 );
1876 case 0: return ( maxInd==6 || maxInd==8 || maxInd==12 );
1877 case 1: return ( maxInd==6 || maxInd==7 || maxInd==13 );
1878 case 2: return ( maxInd==7 || maxInd==8 || maxInd==14 );
1879 case 3: return ( maxInd==9 || maxInd==11 || maxInd==12 );
1880 case 4: return ( maxInd==9 || maxInd==10 || maxInd==13 );
1881 case 5: return ( maxInd==10 || maxInd==11 || maxInd==14 );
1888 const int diff = maxInd-minInd;
1889 if ( diff > 6 ) return false;// not linked top and bottom
1890 if ( diff == 6 ) return true; // linked top and bottom
1891 return diff == 1 || diff == 7;
1898 //=======================================================================
1899 //function : GetNodeIndex
1900 //purpose : Return an index of theNode
1901 //=======================================================================
1903 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1906 for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
1907 if ( myVolumeNodes[ i ] == theNode )
1914 //================================================================================
1916 * \brief Fill vector with boundary faces existing in the mesh
1917 * \param faces - vector of found nodes
1918 * \retval int - nb of found faces
1920 //================================================================================
1922 int SMDS_VolumeTool::GetAllExistingFaces(std::vector<const SMDS_MeshElement*> & faces) const
1925 SaveFacet savedFacet( myCurFace );
1927 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1929 if ( const SMDS_MeshElement* face = SMDS_Mesh::FindFace( myCurFace.myNodes ))
1930 faces.push_back( face );
1933 for ( int iF = 0; iF < NbFaces(); ++iF ) {
1934 const SMDS_MeshFace* face = 0;
1935 const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1936 switch ( NbFaceNodes( iF )) {
1938 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1940 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1942 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1943 nodes[3], nodes[4], nodes[5]); break;
1945 face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1946 nodes[4], nodes[5], nodes[6], nodes[7]); break;
1949 faces.push_back( face );
1951 return faces.size();
1955 //================================================================================
1957 * \brief Fill vector with boundary edges existing in the mesh
1958 * \param edges - vector of found edges
1959 * \retval int - nb of found faces
1961 //================================================================================
1963 int SMDS_VolumeTool::GetAllExistingEdges(std::vector<const SMDS_MeshElement*> & edges) const
1966 edges.reserve( myVolumeNodes.size() * 2 );
1967 for ( size_t i = 0; i < myVolumeNodes.size()-1; ++i ) {
1968 for ( size_t j = i + 1; j < myVolumeNodes.size(); ++j ) {
1969 if ( IsLinked( i, j )) {
1970 const SMDS_MeshElement* edge =
1971 SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1973 edges.push_back( edge );
1977 return edges.size();
1980 //================================================================================
1982 * \brief Return minimal square distance between connected corner nodes
1984 //================================================================================
1986 double SMDS_VolumeTool::MinLinearSize2() const
1988 double minSize = 1e+100;
1989 int iQ = myVolume->IsQuadratic() ? 2 : 1;
1991 SaveFacet savedFacet( myCurFace );
1993 // it seems that compute distance twice is faster than organization of a sole computing
1994 myCurFace.myIndex = -1;
1995 for ( int iF = 0; iF < myNbFaces; ++iF )
1998 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
2000 XYZ n1( myCurFace.myNodes[ iN ]);
2001 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
2002 minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
2009 //================================================================================
2011 * \brief Return maximal square distance between connected corner nodes
2013 //================================================================================
2015 double SMDS_VolumeTool::MaxLinearSize2() const
2017 double maxSize = -1e+100;
2018 int iQ = myVolume->IsQuadratic() ? 2 : 1;
2020 SaveFacet savedFacet( myCurFace );
2022 // it seems that compute distance twice is faster than organization of a sole computing
2023 myCurFace.myIndex = -1;
2024 for ( int iF = 0; iF < myNbFaces; ++iF )
2027 for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
2029 XYZ n1( myCurFace.myNodes[ iN ]);
2030 XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
2031 maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude());
2038 //================================================================================
2040 * \brief Fast quickly check that only one volume is built on the face nodes
2041 * This check is valid for conformal meshes only
2043 //================================================================================
2045 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
2047 const bool isFree = true;
2049 if ( !setFace( faceIndex ))
2052 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
2054 const int di = myVolume->IsQuadratic() ? 2 : 1;
2055 const int nbN = ( myCurFace.myNbNodes/di <= 4 && !IsPoly()) ? 3 : myCurFace.myNbNodes/di; // nb nodes to check
2057 SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
2058 while ( eIt->more() )
2060 const SMDS_MeshElement* vol = eIt->next();
2061 if ( vol == myVolume )
2064 for ( iN = 1; iN < nbN; ++iN )
2065 if ( vol->GetNodeIndex( nodes[ iN*di ]) < 0 )
2067 if ( iN == nbN ) // nbN nodes are shared with vol
2069 // if ( vol->IsPoly() || vol->NbFaces() > 6 ) // vol is polyhed or hex prism
2071 // int nb = myCurFace.myNbNodes;
2072 // if ( myVolume->GetEntityType() != vol->GetEntityType() )
2073 // nb -= ( GetCenterNodeIndex(0) > 0 );
2074 // std::set<const SMDS_MeshNode*> faceNodes( nodes, nodes + nb );
2075 // if ( SMDS_VolumeTool( vol ).GetFaceIndex( faceNodes ) < 0 )
2078 if ( otherVol ) *otherVol = vol;
2082 if ( otherVol ) *otherVol = 0;
2086 //================================================================================
2088 * \brief Thorough check that only one volume is built on the face nodes
2090 //================================================================================
2092 bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
2094 const bool isFree = true;
2096 if (!setFace( faceIndex ))
2099 const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
2100 const int nbFaceNodes = myCurFace.myNbNodes;
2102 // evaluate nb of face nodes shared by other volumes
2103 int maxNbShared = -1;
2104 typedef std::map< const SMDS_MeshElement*, int > TElemIntMap;
2105 TElemIntMap volNbShared;
2106 TElemIntMap::iterator vNbIt;
2107 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
2108 const SMDS_MeshNode* n = nodes[ iNode ];
2109 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
2110 while ( eIt->more() ) {
2111 const SMDS_MeshElement* elem = eIt->next();
2112 if ( elem != myVolume ) {
2113 vNbIt = volNbShared.insert( std::make_pair( elem, 0 )).first;
2115 if ( vNbIt->second > maxNbShared )
2116 maxNbShared = vNbIt->second;
2120 if ( maxNbShared < 3 )
2121 return isFree; // is free
2123 // find volumes laying on the opposite side of the face
2124 // and sharing all nodes
2125 XYZ intNormal; // internal normal
2126 GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
2127 if ( IsFaceExternal( faceIndex ))
2128 intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
2129 XYZ p0 ( nodes[0] ), baryCenter;
2130 for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); ) {
2131 const int& nbShared = (*vNbIt).second;
2132 if ( nbShared >= 3 ) {
2133 SMDS_VolumeTool volume( (*vNbIt).first );
2134 volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
2135 XYZ intNormal2( baryCenter - p0 );
2136 if ( intNormal.Dot( intNormal2 ) < 0 ) {
2138 if ( nbShared >= nbFaceNodes )
2140 // a volume shares the whole facet
2141 if ( otherVol ) *otherVol = vNbIt->first;
2148 // remove a volume from volNbShared map
2149 volNbShared.erase( vNbIt++ );
2152 // here volNbShared contains only volumes laying on the opposite side of
2153 // the face and sharing 3 or more but not all face nodes with myVolume
2154 if ( volNbShared.size() < 2 ) {
2155 return isFree; // is free
2158 // check if the whole area of a face is shared
2159 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
2161 const SMDS_MeshNode* n = nodes[ iNode ];
2162 // check if n is shared by one of volumes of volNbShared
2163 bool isShared = false;
2164 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
2165 while ( eIt->more() && !isShared )
2166 isShared = volNbShared.count( eIt->next() );
2170 if ( otherVol ) *otherVol = volNbShared.begin()->first;
2173 // if ( !myVolume->IsPoly() )
2175 // bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
2176 // for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
2177 // SMDS_VolumeTool volume( (*vNbIt).first );
2178 // bool prevLinkShared = false;
2179 // int nbSharedLinks = 0;
2180 // for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
2181 // bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
2182 // if ( linkShared )
2184 // if ( linkShared && prevLinkShared &&
2185 // volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
2186 // isShared[ iNode ] = true;
2187 // prevLinkShared = linkShared;
2189 // if ( nbSharedLinks == nbFaceNodes )
2190 // return !free; // is not free
2191 // if ( nbFaceNodes == 4 ) {
2192 // // check traingle parts 1 & 3
2193 // if ( isShared[1] && isShared[3] )
2194 // return !free; // is not free
2195 // // check triangle parts 0 & 2;
2196 // // 0 part could not be checked in the loop; check it here
2197 // if ( isShared[2] && prevLinkShared &&
2198 // volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
2199 // volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
2200 // return !free; // is not free
2207 //=======================================================================
2208 //function : GetFaceIndex
2209 //purpose : Return index of a face formed by theFaceNodes
2210 //=======================================================================
2212 int SMDS_VolumeTool::GetFaceIndex( const std::set<const SMDS_MeshNode*>& theFaceNodes,
2213 const int theFaceIndexHint ) const
2215 if ( theFaceIndexHint >= 0 )
2217 int nbNodes = NbFaceNodes( theFaceIndexHint );
2218 if ( nbNodes == (int) theFaceNodes.size() )
2220 const SMDS_MeshNode** nodes = GetFaceNodes( theFaceIndexHint );
2222 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
2227 return theFaceIndexHint;
2230 for ( int iFace = 0; iFace < myNbFaces; iFace++ )
2232 if ( iFace == theFaceIndexHint )
2234 int nbNodes = NbFaceNodes( iFace );
2235 if ( nbNodes == (int) theFaceNodes.size() )
2237 const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
2239 if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
2250 //=======================================================================
2251 //function : GetFaceIndex
2252 //purpose : Return index of a face formed by theFaceNodes
2253 //=======================================================================
2255 /*int SMDS_VolumeTool::GetFaceIndex( const std::set<int>& theFaceNodesIndices )
2257 for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
2258 const int* nodes = GetFaceNodesIndices( iFace );
2259 int nbFaceNodes = NbFaceNodes( iFace );
2260 std::set<int> nodeSet;
2261 for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
2262 nodeSet.insert( nodes[ iNode ] );
2263 if ( theFaceNodesIndices == nodeSet )
2269 //=======================================================================
2270 //function : setFace
2272 //=======================================================================
2274 bool SMDS_VolumeTool::setFace( int faceIndex ) const
2279 if ( myCurFace.myIndex == faceIndex )
2282 myCurFace.myIndex = -1;
2284 if ( faceIndex < 0 || faceIndex >= NbFaces() )
2287 if (myVolume->IsPoly())
2289 if ( !myPolyedre ) {
2290 MESSAGE("Warning: bad volumic element");
2295 SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
2296 if ( !myAllFacesNbNodes ) {
2297 me->myPolyQuantities = myPolyedre->GetQuantities();
2298 myAllFacesNbNodes = &myPolyQuantities[0];
2300 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2301 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2302 me->myPolyIndices.resize( myCurFace.myNbNodes + 1 );
2303 myCurFace.myNodeIndices = & me->myPolyIndices[0];
2304 int shift = std::accumulate( myAllFacesNbNodes, myAllFacesNbNodes+faceIndex, 0 );
2305 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2307 myCurFace.myNodes [ iNode ] = myVolumeNodes[ shift + iNode ];
2308 myCurFace.myNodeIndices[ iNode ] = shift + iNode;
2310 myCurFace.myNodes [ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ]; // last = first
2311 myCurFace.myNodeIndices[ myCurFace.myNbNodes ] = myCurFace.myNodeIndices[ 0 ];
2313 // check orientation
2314 if (myExternalFaces)
2316 myCurFace.myIndex = faceIndex; // avoid infinite recursion in IsFaceExternal()
2317 myExternalFaces = false; // force normal computation by IsFaceExternal()
2318 if ( !IsFaceExternal( faceIndex ))
2319 std::reverse( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
2320 myExternalFaces = true;
2325 if ( !myAllFacesNodeIndices_F )
2327 // choose data for an element type
2328 switch ( myVolumeNodes.size() ) {
2330 myAllFacesNodeIndices_F = &Tetra_F [0][0];
2331 //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
2332 myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
2333 myAllFacesNbNodes = Tetra_nbN;
2334 myMaxFaceNbNodes = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
2337 myAllFacesNodeIndices_F = &Pyramid_F [0][0];
2338 //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
2339 myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
2340 myAllFacesNbNodes = Pyramid_nbN;
2341 myMaxFaceNbNodes = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
2344 myAllFacesNodeIndices_F = &Penta_F [0][0];
2345 //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
2346 myAllFacesNodeIndices_RE = &Penta_RE[0][0];
2347 myAllFacesNbNodes = Penta_nbN;
2348 myMaxFaceNbNodes = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
2351 myAllFacesNodeIndices_F = &Hexa_F [0][0];
2352 ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
2353 myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
2354 myAllFacesNbNodes = Hexa_nbN;
2355 myMaxFaceNbNodes = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
2358 myAllFacesNodeIndices_F = &QuadTetra_F [0][0];
2359 //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
2360 myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
2361 myAllFacesNbNodes = QuadTetra_nbN;
2362 myMaxFaceNbNodes = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
2365 myAllFacesNodeIndices_F = &QuadPyram_F [0][0];
2366 //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
2367 myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
2368 myAllFacesNbNodes = QuadPyram_nbN;
2369 myMaxFaceNbNodes = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
2372 myAllFacesNodeIndices_F = &QuadPenta_F [0][0];
2373 //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
2374 myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
2375 myAllFacesNbNodes = QuadPenta_nbN;
2376 myMaxFaceNbNodes = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
2380 myAllFacesNodeIndices_F = &QuadHexa_F [0][0];
2381 //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
2382 myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
2383 myAllFacesNbNodes = QuadHexa_nbN;
2384 myMaxFaceNbNodes = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
2385 if ( !myIgnoreCentralNodes && myVolumeNodes.size() == 27 )
2387 myAllFacesNodeIndices_F = &TriQuadHexa_F [0][0];
2388 //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
2389 myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
2390 myAllFacesNbNodes = TriQuadHexa_nbN;
2391 myMaxFaceNbNodes = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
2395 myAllFacesNodeIndices_F = &HexPrism_F [0][0];
2396 //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
2397 myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
2398 myAllFacesNbNodes = HexPrism_nbN;
2399 myMaxFaceNbNodes = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
2405 myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2406 // if ( myExternalFaces )
2407 // myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2409 // myCurFace.myNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
2410 myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2413 myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2414 for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2415 myCurFace.myNodes[ iNode ] = myVolumeNodes[ myCurFace.myNodeIndices[ iNode ]];
2416 myCurFace.myNodes[ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ];
2419 myCurFace.myIndex = faceIndex;
2424 //=======================================================================
2425 //function : GetType
2426 //purpose : return VolumeType by nb of nodes in a volume
2427 //=======================================================================
2429 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
2431 switch ( nbNodes ) {
2432 case 4: return TETRA;
2433 case 5: return PYRAM;
2434 case 6: return PENTA;
2435 case 8: return HEXA;
2436 case 10: return QUAD_TETRA;
2437 case 13: return QUAD_PYRAM;
2438 case 15: return QUAD_PENTA;
2440 case 27: return QUAD_HEXA;
2441 case 12: return HEX_PRISM;
2442 default:return UNKNOWN;
2446 //=======================================================================
2447 //function : NbFaces
2448 //purpose : return nb of faces by volume type
2449 //=======================================================================
2451 int SMDS_VolumeTool::NbFaces( VolumeType type )
2455 case QUAD_TETRA: return 4;
2457 case QUAD_PYRAM: return 5;
2459 case QUAD_PENTA: return 5;
2461 case QUAD_HEXA : return 6;
2462 case HEX_PRISM : return 8;
2467 //================================================================================
2469 * \brief Useful to know nb of corner nodes of a quadratic volume
2470 * \param type - volume type
2471 * \retval int - nb of corner nodes
2473 //================================================================================
2475 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
2479 case QUAD_TETRA: return 4;
2481 case QUAD_PYRAM: return 5;
2483 case QUAD_PENTA: return 6;
2485 case QUAD_HEXA : return 8;
2486 case HEX_PRISM : return 12;
2493 //=======================================================================
2494 //function : GetFaceNodesIndices
2495 //purpose : Return the array of face nodes indices
2496 // To comfort link iteration, the array
2497 // length == NbFaceNodes( faceIndex ) + 1 and
2498 // the last node index == the first one.
2499 //=======================================================================
2501 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
2506 case TETRA: return Tetra_F[ faceIndex ];
2507 case PYRAM: return Pyramid_F[ faceIndex ];
2508 case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
2509 case HEXA: return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
2510 case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
2511 case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
2512 case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
2513 // what about SMDSEntity_TriQuad_Hexa?
2514 case QUAD_HEXA: return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
2515 case HEX_PRISM: return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
2521 //=======================================================================
2522 //function : NbFaceNodes
2523 //purpose : Return number of nodes in the array of face nodes
2524 //=======================================================================
2526 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
2530 case TETRA: return Tetra_nbN[ faceIndex ];
2531 case PYRAM: return Pyramid_nbN[ faceIndex ];
2532 case PENTA: return Penta_nbN[ faceIndex ];
2533 case HEXA: return Hexa_nbN[ faceIndex ];
2534 case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
2535 case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
2536 case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
2537 // what about SMDSEntity_TriQuad_Hexa?
2538 case QUAD_HEXA: return QuadHexa_nbN[ faceIndex ];
2539 case HEX_PRISM: return HexPrism_nbN[ faceIndex ];
2545 //=======================================================================
2546 //function : Element
2547 //purpose : return element
2548 //=======================================================================
2550 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
2552 return static_cast<const SMDS_MeshVolume*>( myVolume );
2555 //=======================================================================
2557 //purpose : return element ID
2558 //=======================================================================
2560 smIdType SMDS_VolumeTool::ID() const
2562 return myVolume ? myVolume->GetID() : 0;