X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMDS%2FSMDS_VolumeTool.cxx;h=d773413e34c28da55946b08c8d9f61b067462e06;hp=22fbb84202d6f739d30da59cb9eabad8588813c8;hb=b9b3af112da57064b63a2b2c9644edcc808c1116;hpb=52f87e4dea822cded355f9eaa3b96df054904d23 diff --git a/src/SMDS/SMDS_VolumeTool.cxx b/src/SMDS/SMDS_VolumeTool.cxx index 22fbb8420..d773413e3 100644 --- a/src/SMDS/SMDS_VolumeTool.cxx +++ b/src/SMDS/SMDS_VolumeTool.cxx @@ -1,8 +1,29 @@ +// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + // File : SMDS_VolumeTool.cxx // Created : Tue Jul 13 12:22:13 2004 // Author : Edward AGAPOV (eap) -// Copyright : Open CASCADE - +// #ifdef _MSC_VER #pragma warning(disable:4786) #endif @@ -11,13 +32,14 @@ #include "SMDS_MeshElement.hxx" #include "SMDS_MeshNode.hxx" -#include "SMDS_PolyhedralVolumeOfNodes.hxx" +#include "SMDS_VtkVolume.hxx" +#include "SMDS_Mesh.hxx" #include "utilities.h" #include -#include -#include +#include +#include using namespace std; @@ -25,6 +47,10 @@ using namespace std; // Node indices in faces depending on volume orientation // making most faces normals external // ====================================================== +// For all elements, 0-th face is bottom based on the first nodes. +// For prismatic elements (tetra,hexa,prisms), 1-th face is a top one. +// For all elements, side faces follow order of bottom nodes +// ====================================================== /* // N3 @@ -45,11 +71,6 @@ static int Tetra_F [4][4] = { // FORWARD == EXTERNAL { 0, 3, 1, 0 }, { 1, 3, 2, 1 }, { 0, 2, 3, 0 }}; -static int Tetra_R [4][4] = { // REVERSED - { 0, 1, 2, 0 }, // All faces but a bottom have external normals - { 0, 1, 3, 0 }, - { 1, 2, 3, 1 }, - { 0, 3, 2, 0 }}; static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL) { 0, 2, 1, 0 }, // All faces have external normals { 0, 1, 3, 0 }, @@ -65,13 +86,8 @@ static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL { 0, 4, 1, 0, 4 }, { 1, 4, 2, 1, 4 }, { 2, 4, 3, 2, 4 }, - { 3, 4, 0, 3, 4 }}; -static int Pyramid_R [5][5] = { // REVERSED - { 0, 1, 2, 3, 0 }, // All faces but a bottom have external normals - { 0, 1, 4, 0, 4 }, - { 1, 2, 4, 1, 4 }, - { 2, 3, 4, 2, 4 }, - { 3, 0, 4, 3, 4 }}; + { 3, 4, 0, 3, 4 } +}; static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL) { 0, 3, 2, 1, 0 }, // All faces but a bottom have external normals { 0, 1, 4, 0, 4 }, @@ -96,29 +112,17 @@ static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 }; // N0 +---------+ N2 */ static int Penta_F [5][5] = { // FORWARD - { 0, 1, 2, 0, 0 }, // Top face has an internal normal, other - external - { 3, 4, 5, 3, 3 }, // 0 is bottom, 1 is top face - { 0, 2, 5, 3, 0 }, - { 1, 4, 5, 2, 1 }, - { 0, 3, 4, 1, 0 }}; -static int Penta_R [5][5] = { // REVERSED - { 0, 1, 2, 0, 0 }, // Bottom face has an internal normal, other - external - { 3, 4, 5, 3, 3 }, // 0 is bottom, 1 is top face - { 0, 3, 5, 2, 0 }, - { 1, 2, 5, 4, 1 }, - { 0, 1, 4, 3, 0 }}; -static int Penta_FE [5][5] = { // FORWARD -> EXTERNAL - { 0, 1, 2, 0, 0 }, - { 3, 5, 4, 3, 3 }, - { 0, 2, 5, 3, 0 }, + { 0, 1, 2, 0, 0 }, // All faces have external normals + { 3, 5, 4, 3, 3 }, // 0 is bottom, 1 is top face + { 0, 3, 4, 1, 0 }, { 1, 4, 5, 2, 1 }, - { 0, 3, 4, 1, 0 }}; + { 0, 2, 5, 3, 0 }}; static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL { 0, 2, 1, 0, 0 }, { 3, 4, 5, 3, 3 }, - { 0, 3, 5, 2, 0 }, + { 0, 1, 4, 3, 0 }, { 1, 2, 5, 4, 1 }, - { 0, 1, 4, 3, 0 }}; + { 0, 3, 5, 2, 0 }}; static int Penta_nbN [] = { 3, 3, 4, 4, 4 }; /* @@ -128,8 +132,6 @@ static int Penta_nbN [] = { 3, 3, 4, 4, 4 }; // / | / | // N4+----------+N7 | // | | | | HEXAHEDRON -// | | | | -// | | | | // | N1+------|---+N2 // | / | / // | / | / @@ -137,38 +139,218 @@ static int Penta_nbN [] = { 3, 3, 4, 4, 4 }; // N0+----------+N3 */ static int Hexa_F [6][5] = { // FORWARD - { 0, 1, 2, 3, 0 }, // opposite faces are neighbouring, - { 4, 5, 6, 7, 4 }, // odd face(1,3,5) normal is internal, even(0,2,4) - external - { 1, 0, 4, 5, 1 }, // same index nodes of opposite faces are linked - { 2, 3, 7, 6, 2 }, - { 0, 3, 7, 4, 0 }, - { 1, 2, 6, 5, 1 }}; -// static int Hexa_R [6][5] = { // REVERSED -// { 0, 3, 2, 1, 0 }, // opposite faces are neighbouring, -// { 4, 7, 6, 5, 4 }, // odd face(1,3,5) normal is external, even(0,2,4) - internal -// { 1, 5, 4, 0, 1 }, // same index nodes of opposite faces are linked -// { 2, 6, 7, 3, 2 }, -// { 0, 4, 7, 3, 0 }, -// { 1, 5, 6, 2, 1 }}; -static int Hexa_FE [6][5] = { // FORWARD -> EXTERNAL - { 0, 1, 2, 3, 0 } , // opposite faces are neighbouring, - { 4, 7, 6, 5, 4 }, // all face normals are external, - { 0, 4, 5, 1, 0 }, // links in opposite faces: 0-0, 1-3, 2-2, 3-1 + { 0, 1, 2, 3, 0 }, + { 4, 7, 6, 5, 4 }, // all face normals are external + { 0, 4, 5, 1, 0 }, + { 1, 5, 6, 2, 1 }, { 3, 2, 6, 7, 3 }, - { 0, 3, 7, 4, 0 }, - { 1, 5, 6, 2, 1 }}; + { 0, 3, 7, 4, 0 }}; static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL - { 0, 3, 2, 1, 0 }, // opposite faces are neighbouring, - { 4, 5, 6, 7, 4 }, // all face normals are external, - { 0, 1, 5, 4, 0 }, // links in opposite faces: 0-0, 1-3, 2-2, 3-1 + { 0, 3, 2, 1, 0 }, + { 4, 5, 6, 7, 4 }, // all face normals are external + { 0, 1, 5, 4, 0 }, + { 1, 2, 6, 5, 1 }, { 3, 7, 6, 2, 3 }, - { 0, 4, 7, 3, 0 }, - { 1, 2, 6, 5, 1 }}; + { 0, 4, 7, 3, 0 }}; static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 }; +static int Hexa_oppF[] = { 1, 0, 4, 5, 2, 3 }; // oppopsite facet indices + +/* +// N8 +------+ N9 +// / \ +// / \ +// N7 + + N10 +// \ / +// \ / +// N6 +------+ N11 +// HEXAGONAL PRISM +// N2 +------+ N3 +// / \ +// / \ +// N1 + + N4 +// \ / +// \ / +// N0 +------+ N5 +*/ +static int HexPrism_F [8][7] = { // FORWARD + { 0, 1, 2, 3, 4, 5, 0 }, + { 6,11,10, 9, 8, 7, 6 }, + { 0, 6, 7, 1, 0, 0, 0 }, + { 1, 7, 8, 2, 1, 1, 1 }, + { 2, 8, 9, 3, 2, 2, 2 }, + { 3, 9,10, 4, 3, 3, 3 }, + { 4,10,11, 5, 4, 4, 4 }, + { 5,11, 6, 0, 5, 5, 5 }}; +static int HexPrism_RE [8][7] = { // REVERSED -> EXTERNAL + { 0, 5, 4, 3, 2, 1, 0 }, + { 6,11,10, 9, 8, 7, 6 }, + { 0, 6, 7, 1, 0, 0, 0 }, + { 1, 7, 8, 2, 1, 1, 1 }, + { 2, 8, 9, 3, 2, 2, 2 }, + { 3, 9,10, 4, 3, 3, 3 }, + { 4,10,11, 5, 4, 4, 4 }, + { 5,11, 6, 0, 5, 5, 5 }}; +static int HexPrism_nbN [] = { 6, 6, 4, 4, 4, 4, 4, 4 }; + + +/* +// N3 +// + +// /|\ +// 7/ | \8 +// / |4 \ QUADRATIC +// N0 +---|---+ N1 TETRAHEDRON +// \ +9 / +// \ | / +// 6\ | /5 +// \|/ +// + +// N2 +*/ +static int QuadTetra_F [4][7] = { // FORWARD + { 0, 4, 1, 5, 2, 6, 0 }, // All faces have external normals + { 0, 7, 3, 8, 1, 4, 0 }, + { 1, 8, 3, 9, 2, 5, 1 }, + { 0, 6, 2, 9, 3, 7, 0 }}; +static int QuadTetra_RE [4][7] = { // REVERSED -> FORWARD (EXTERNAL) + { 0, 6, 2, 5, 1, 4, 0 }, // All faces have external normals + { 0, 4, 1, 8, 3, 7, 0 }, + { 1, 5, 2, 9, 3, 8, 1 }, + { 0, 7, 3, 9, 2, 6, 0 }}; +static int QuadTetra_nbN [] = { 6, 6, 6, 6 }; + +// +// QUADRATIC +// PYRAMID +// +// +4 +// +// +// 10+-----+11 +// | | 9 - middle point for (0,4) etc. +// | | +// 9+-----+12 +// +// 6 +// 1+----+----+2 +// | | +// | | +// 5+ +7 +// | | +// | | +// 0+----+----+3 +// 8 +static int QuadPyram_F [5][9] = { // FORWARD + { 0, 5, 1, 6, 2, 7, 3, 8, 0 }, // All faces have external normals + { 0, 9, 4, 10,1, 5, 0, 4, 4 }, + { 1, 10,4, 11,2, 6, 1, 4, 4 }, + { 2, 11,4, 12,3, 7, 2, 4, 4 }, + { 3, 12,4, 9, 0, 8, 3, 4, 4 }}; +static int QuadPyram_RE [5][9] = { // REVERSED -> FORWARD (EXTERNAL) + { 0, 8, 3, 7, 2, 6, 1, 5, 0 }, // All faces but a bottom have external normals + { 0, 5, 1, 10,4, 9, 0, 4, 4 }, + { 1, 6, 2, 11,4, 10,1, 4, 4 }, + { 2, 7, 3, 12,4, 11,2, 4, 4 }, + { 3, 8, 0, 9, 4, 12,3, 4, 4 }}; +static int QuadPyram_nbN [] = { 8, 6, 6, 6, 6 }; + +/* +// + N4 +// /|\ +// 9/ | \10 +// / | \ +// / | \ +// N3 +----+----+ N5 +// | |11 | +// | | | +// | +13 | QUADRATIC +// | | | PENTAHEDRON +// 12+ | +14 +// | | | +// | | | +// | + N1 | +// | / \ | +// | 6/ \7 | +// | / \ | +// |/ \| +// N0 +---------+ N2 +// 8 +*/ +static int QuadPenta_F [5][9] = { // FORWARD + { 0, 6, 1, 7, 2, 8, 0, 0, 0 }, + { 3, 11,5, 10,4, 9, 3, 3, 3 }, + { 0, 12,3, 9, 4, 13,1, 6, 0 }, + { 1, 13,4, 10,5, 14,2, 7, 1 }, + { 0, 8, 2, 14,5, 11,3, 12,0 }}; +static int QuadPenta_RE [5][9] = { // REVERSED -> EXTERNAL + { 0, 8, 2, 7, 1, 6, 0, 0, 0 }, + { 3, 9, 4, 10,5, 11,3, 3, 3 }, + { 0, 6, 1, 13,4, 9, 3, 12,0 }, + { 1, 7, 2, 14,5, 10,4, 13,1 }, + { 0, 12,3, 11,5, 14,2, 8, 0 }}; +static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 }; + +/* +// 13 +// N5+-----+-----+N6 +-----+-----+ +// /| /| /| /| +// 12+ | 14+ | + | +25 + | +// / | / | / | / | +// N4+-----+-----+N7 | QUADRATIC +-----+-----+ | Central nodes +// | | 15 | | HEXAHEDRON | | | | of tri-quadratic +// | | | | | | | | HEXAHEDRON +// | 17+ | +18 | + 22+ | + +// | | | | |21 | | | +// | | | | | + | 26+ | + | +// | | | | | | |23 | +// 16+ | +19 | + | +24 + | +// | | | | | | | | +// | | 9 | | | | | | +// | N1+-----+-|---+N2 | +-----+-|---+ +// | / | / | / | / +// | +8 | +10 | + 20+ | + +// |/ |/ |/ |/ +// N0+-----+-----+N3 +-----+-----+ +// 11 +*/ +static int QuadHexa_F [6][9] = { // FORWARD + { 0, 8, 1, 9, 2, 10,3, 11,0 }, // all face normals are external, + { 4, 15,7, 14,6, 13,5, 12,4 }, + { 0, 16,4, 12,5, 17,1, 8, 0 }, + { 1, 17,5, 13,6, 18,2, 9, 1 }, + { 3, 10,2, 18,6, 14,7, 19,3 }, + { 0, 11,3, 19,7, 15,4, 16,0 }}; +static int QuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL + { 0, 11,3, 10,2, 9, 1, 8, 0 }, // all face normals are external + { 4, 12,5, 13,6, 14,7, 15,4 }, + { 0, 8, 1, 17,5, 12,4, 16,0 }, + { 1, 9, 2, 18,6, 13,5, 17,1 }, + { 3, 19,7, 14,6, 18,2, 10,3 }, + { 0, 16,4, 15,7, 19,3, 11,0 }}; +static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 }; + +static int TriQuadHexa_F [6][9] = { // FORWARD + { 0, 8, 1, 9, 2, 10,3, 11, 20 }, // all face normals are external + { 4, 15,7, 14,6, 13,5, 12, 25 }, + { 0, 16,4, 12,5, 17,1, 8, 21 }, + { 1, 17,5, 13,6, 18,2, 9, 22 }, + { 3, 10,2, 18,6, 14,7, 19, 23 }, + { 0, 11,3, 19,7, 15,4, 16, 24 }}; +static int TriQuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL + { 0, 11,3, 10,2, 9, 1, 8, 20 }, // opposite faces are neighbouring, + { 4, 12,5, 13,6, 14,7, 15, 25 }, // all face normals are external + { 0, 8, 1, 17,5, 12,4, 16, 21 }, + { 1, 9, 2, 18,6, 13,5, 17, 22 }, + { 3, 19,7, 14,6, 18,2, 10, 23 }, + { 0, 16,4, 15,7, 19,3, 11, 24 }}; +static int TriQuadHexa_nbN [] = { 9, 9, 9, 9, 9, 9 }; + // ======================================================== // to perform some calculations without linkage to CASCADE // ======================================================== +namespace +{ struct XYZ { double x; double y; @@ -177,25 +359,56 @@ struct XYZ { XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; } XYZ( const XYZ& other ) { x = other.x; y = other.y; z = other.z; } XYZ( const SMDS_MeshNode* n ) { x = n->X(); y = n->Y(); z = n->Z(); } - XYZ operator-( const XYZ& other ); - XYZ Crossed( const XYZ& other ); - double Dot( const XYZ& other ); - double Magnitude(); + inline XYZ operator-( const XYZ& other ); + inline XYZ operator+( const XYZ& other ); + inline XYZ Crossed( const XYZ& other ); + inline double Dot( const XYZ& other ); + inline double Magnitude(); + inline double SquareMagnitude(); }; -XYZ XYZ::operator-( const XYZ& Right ) { +inline XYZ XYZ::operator-( const XYZ& Right ) { return XYZ(x - Right.x, y - Right.y, z - Right.z); } -XYZ XYZ::Crossed( const XYZ& Right ) { +inline XYZ XYZ::operator+( const XYZ& Right ) { + return XYZ(x + Right.x, y + Right.y, z + Right.z); +} +inline XYZ XYZ::Crossed( const XYZ& Right ) { return XYZ (y * Right.z - z * Right.y, z * Right.x - x * Right.z, x * Right.y - y * Right.x); } -double XYZ::Dot( const XYZ& Other ) { +inline double XYZ::Dot( const XYZ& Other ) { return(x * Other.x + y * Other.y + z * Other.z); } -double XYZ::Magnitude() { +inline double XYZ::Magnitude() { return sqrt (x * x + y * y + z * z); } +inline double XYZ::SquareMagnitude() { + return (x * x + y * y + z * z); +} + + //================================================================================ + /*! + * \brief Return linear type corresponding to a quadratic one + */ + //================================================================================ + + SMDS_VolumeTool::VolumeType quadToLinear(SMDS_VolumeTool::VolumeType quadType) + { + SMDS_VolumeTool::VolumeType linType = SMDS_VolumeTool::VolumeType( int(quadType)-4 ); + const int nbCornersByQuad = SMDS_VolumeTool::NbCornerNodes( quadType ); + if ( SMDS_VolumeTool::NbCornerNodes( linType ) == nbCornersByQuad ) + return linType; + + int iLin = 0; + for ( ; iLin < SMDS_VolumeTool::NB_VOLUME_TYPES; ++iLin ) + if ( SMDS_VolumeTool::NbCornerNodes( SMDS_VolumeTool::VolumeType( iLin )) == nbCornersByQuad) + return SMDS_VolumeTool::VolumeType( iLin ); + + return SMDS_VolumeTool::UNKNOWN; + } + +} // namespace //======================================================================= //function : SMDS_VolumeTool @@ -203,18 +416,10 @@ double XYZ::Magnitude() { //======================================================================= SMDS_VolumeTool::SMDS_VolumeTool () - : myVolume( 0 ), - myPolyedre( 0 ), - myVolForward( true ), - myNbFaces( 0 ), - myVolumeNbNodes( 0 ), - myVolumeNodes( NULL ), - myExternalFaces( false ), - myFaceNbNodes( 0 ), - myCurFace( -1 ), - myFaceNodeIndices( NULL ), - myFaceNodes( NULL ) + : myVolumeNodes( NULL ), + myFaceNodes( NULL ) { + Set( 0 ); } //======================================================================= @@ -222,20 +427,12 @@ SMDS_VolumeTool::SMDS_VolumeTool () //purpose : //======================================================================= -SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume) - : myVolume( 0 ), - myPolyedre( 0 ), - myVolForward( true ), - myNbFaces( 0 ), - myVolumeNbNodes( 0 ), - myVolumeNodes( NULL ), - myExternalFaces( false ), - myFaceNbNodes( 0 ), - myCurFace( -1 ), - myFaceNodeIndices( NULL ), - myFaceNodes( NULL ) +SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume, + const bool ignoreCentralNodes) + : myVolumeNodes( NULL ), + myFaceNodes( NULL ) { - Set( theVolume ); + Set( theVolume, ignoreCentralNodes ); } //======================================================================= @@ -245,14 +442,11 @@ SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume) SMDS_VolumeTool::~SMDS_VolumeTool() { - if (myVolumeNodes != NULL) { - delete [] myVolumeNodes; - myVolumeNodes = NULL; - } - if (myFaceNodes != NULL) { - delete [] myFaceNodes; - myFaceNodes = NULL; - } + if ( myVolumeNodes != NULL ) delete [] myVolumeNodes; + if ( myFaceNodes != NULL ) delete [] myFaceNodes; + + myFaceNodeIndices = NULL; + myVolumeNodes = myFaceNodes = NULL; } //======================================================================= @@ -260,10 +454,13 @@ SMDS_VolumeTool::~SMDS_VolumeTool() //purpose : Set volume to iterate on //======================================================================= -bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume) +bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume, + const bool ignoreCentralNodes) { + // reset fields myVolume = 0; myPolyedre = 0; + myIgnoreCentralNodes = ignoreCentralNodes; myVolForward = true; myNbFaces = 0; @@ -272,61 +469,63 @@ bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume) delete [] myVolumeNodes; myVolumeNodes = NULL; } + myPolyIndices.clear(); myExternalFaces = false; - myFaceNbNodes = 0; + + myAllFacesNodeIndices_F = 0; + //myAllFacesNodeIndices_FE = 0; + myAllFacesNodeIndices_RE = 0; + myAllFacesNbNodes = 0; myCurFace = -1; + myFaceNbNodes = 0; myFaceNodeIndices = NULL; if (myFaceNodes != NULL) { delete [] myFaceNodes; myFaceNodes = NULL; } - if ( theVolume && theVolume->GetType() == SMDSAbs_Volume ) - { - myVolume = theVolume; + // set volume data + if ( !theVolume || theVolume->GetType() != SMDSAbs_Volume ) + return false; - myNbFaces = theVolume->NbFaces(); - myVolumeNbNodes = theVolume->NbNodes(); + myVolume = theVolume; + if (myVolume->IsPoly()) + myPolyedre = dynamic_cast( myVolume ); - // set volume nodes - int iNode = 0; - myVolumeNodes = new const SMDS_MeshNode* [myVolumeNbNodes]; - SMDS_ElemIteratorPtr nodeIt = myVolume->nodesIterator(); - while ( nodeIt->more() ) { - myVolumeNodes[ iNode++ ] = static_cast( nodeIt->next() ); - } + myNbFaces = theVolume->NbFaces(); + myVolumeNbNodes = theVolume->NbNodes(); - if (myVolume->IsPoly()) { - myPolyedre = static_cast( myVolume ); - if (!myPolyedre) { - MESSAGE("Warning: bad volumic element"); - return false; - } - } else { - switch ( myVolumeNbNodes ) { - case 4: - case 5: - case 6: - case 8: { - // define volume orientation - XYZ botNormal; - GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z ); - const SMDS_MeshNode* topNode = myVolumeNodes[ myVolumeNbNodes - 1 ]; - const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ]; - XYZ upDir (topNode->X() - botNode->X(), - topNode->Y() - botNode->Y(), - topNode->Z() - botNode->Z() ); - myVolForward = ( botNormal.Dot( upDir ) < 0 ); - break; - } - default: - break; - } - } + // set nodes + int iNode = 0; + myVolumeNodes = new const SMDS_MeshNode* [myVolumeNbNodes]; + SMDS_ElemIteratorPtr nodeIt = myVolume->nodesIterator(); + while ( nodeIt->more() ) + myVolumeNodes[ iNode++ ] = static_cast( nodeIt->next() ); + + // check validity + if ( !setFace(0) ) + return ( myVolume = 0 ); + + if ( !myPolyedre ) + { + // define volume orientation + XYZ botNormal; + GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z ); + const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ]; + int topNodeIndex = myVolume->NbCornerNodes() - 1; + while ( !IsLinked( 0, topNodeIndex, /*ignoreMediumNodes=*/true )) --topNodeIndex; + const SMDS_MeshNode* topNode = myVolumeNodes[ topNodeIndex ]; + XYZ upDir (topNode->X() - botNode->X(), + topNode->Y() - botNode->Y(), + topNode->Z() - botNode->Z() ); + myVolForward = ( botNormal.Dot( upDir ) < 0 ); + + if ( !myVolForward ) + myCurFace = -1; // previous setFace(0) didn't take myVolForward into account } - return ( myVolume != 0 ); + return true; } //======================================================================= @@ -368,6 +567,51 @@ void SMDS_VolumeTool::Inverse () SWAP_NODES( myVolumeNodes, 1, 3 ); SWAP_NODES( myVolumeNodes, 5, 7 ); break; + case 12: + SWAP_NODES( myVolumeNodes, 1, 5 ); + SWAP_NODES( myVolumeNodes, 2, 4 ); + SWAP_NODES( myVolumeNodes, 7, 11 ); + SWAP_NODES( myVolumeNodes, 8, 10 ); + break; + + case 10: + SWAP_NODES( myVolumeNodes, 1, 2 ); + SWAP_NODES( myVolumeNodes, 4, 6 ); + SWAP_NODES( myVolumeNodes, 8, 9 ); + break; + case 13: + SWAP_NODES( myVolumeNodes, 1, 3 ); + SWAP_NODES( myVolumeNodes, 5, 8 ); + SWAP_NODES( myVolumeNodes, 6, 7 ); + SWAP_NODES( myVolumeNodes, 10, 12 ); + break; + case 15: + SWAP_NODES( myVolumeNodes, 1, 2 ); + SWAP_NODES( myVolumeNodes, 4, 5 ); + SWAP_NODES( myVolumeNodes, 6, 8 ); + SWAP_NODES( myVolumeNodes, 9, 11 ); + SWAP_NODES( myVolumeNodes, 13, 14 ); + break; + case 20: + SWAP_NODES( myVolumeNodes, 1, 3 ); + SWAP_NODES( myVolumeNodes, 5, 7 ); + SWAP_NODES( myVolumeNodes, 8, 11 ); + SWAP_NODES( myVolumeNodes, 9, 10 ); + SWAP_NODES( myVolumeNodes, 12, 15 ); + SWAP_NODES( myVolumeNodes, 13, 14 ); + SWAP_NODES( myVolumeNodes, 17, 19 ); + break; + case 27: + SWAP_NODES( myVolumeNodes, 1, 3 ); + SWAP_NODES( myVolumeNodes, 5, 7 ); + SWAP_NODES( myVolumeNodes, 8, 11 ); + SWAP_NODES( myVolumeNodes, 9, 10 ); + SWAP_NODES( myVolumeNodes, 12, 15 ); + SWAP_NODES( myVolumeNodes, 13, 14 ); + SWAP_NODES( myVolumeNodes, 17, 19 ); + SWAP_NODES( myVolumeNodes, 21, 24 ); + SWAP_NODES( myVolumeNodes, 22, 23 ); + break; default:; } } @@ -382,15 +626,18 @@ SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const if ( myPolyedre ) return POLYHEDA; - if ( myVolume ) { - static const VolumeType types[] = { - TETRA, // myVolumeNbNodes = 4 - PYRAM, // myVolumeNbNodes = 5 - PENTA, // myVolumeNbNodes = 6 - UNKNOWN, // myVolumeNbNodes = 7 - HEXA // myVolumeNbNodes = 8 - }; - return types[ myVolumeNbNodes - 4 ]; + switch( myVolumeNbNodes ) { + case 4: return TETRA; + case 5: return PYRAM; + case 6: return PENTA; + case 8: return HEXA; + case 12: return HEX_PRISM; + case 10: return QUAD_TETRA; + case 13: return QUAD_PYRAM; + case 15: return QUAD_PENTA; + case 20: return QUAD_HEXA; + case 27: return QUAD_HEXA; + default: break; } return UNKNOWN; @@ -450,30 +697,29 @@ double SMDS_VolumeTool::GetSize() const // split a polyhedron into tetrahedrons + int saveCurFace = myCurFace; SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this ); - XYZ baryCenter; - me->GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z); - SMDS_MeshNode bcNode ( baryCenter.x, baryCenter.y, baryCenter.z ); - for ( int f = 0; f < NbFaces(); ++f ) { - bool externalFace = me->IsFaceExternal( f ); // it calls setFace() - for ( int n = 2; n < myFaceNbNodes; ++n ) + me->setFace( f ); + XYZ area (0,0,0), p1( myFaceNodes[0] ); + for ( int n = 0; n < myFaceNbNodes; ++n ) { - double Vn = getTetraVolume( myFaceNodes[ 0 ], - myFaceNodes[ n-1 ], - myFaceNodes[ n ], - & bcNode ); -/// cout <<"++++ " << Vn << " nodes " <GetID() << " " <GetID() << " " <GetID() << " < " << V << endl; - V += externalFace ? -Vn : Vn; + XYZ p2( myFaceNodes[ n+1 ]); + area = area + p1.Crossed( p2 ); + p1 = p2; } + V += p1.Dot( area ); } + V /= 6; + if ( saveCurFace > -1 && saveCurFace != myCurFace ) + me->setFace( myCurFace ); } else { const static int ind[] = { - 0, 1, 3, 6, 11 }; - const static int vtab[][4] = { + 0, 1, 3, 6, 11, 23, 31, 44, 58, 78 }; + const static int vtab[][4] = { // decomposition into tetra in the order of enum VolumeType // tetrahedron { 0, 1, 2, 3 }, // pyramid @@ -488,15 +734,96 @@ double SMDS_VolumeTool::GetSize() const { 4, 1, 6, 5 }, { 1, 3, 6, 2 }, { 4, 6, 3, 7 }, - { 1, 4, 6, 3 } + { 1, 4, 6, 3 }, + // hexagonal prism + { 0, 1, 2, 7 }, + { 0, 7, 8, 6 }, + { 2, 7, 8, 0 }, + + { 0, 3, 4, 9 }, + { 0, 9, 10, 6 }, + { 4, 9, 10, 0 }, + + { 0, 3, 4, 9 }, + { 0, 9, 10, 6 }, + { 4, 9, 10, 0 }, + + { 0, 4, 5, 10 }, + { 0, 10, 11, 6 }, + { 5, 10, 11, 0 }, + + // quadratic tetrahedron + { 0, 4, 6, 7 }, + { 1, 5, 4, 8 }, + { 2, 6, 5, 9 }, + { 7, 8, 9, 3 }, + { 4, 6, 7, 9 }, + { 4, 5, 6, 9 }, + { 4, 7, 8, 9 }, + { 4, 5, 9, 8 }, + + // quadratic pyramid + { 0, 5, 8, 9 }, + { 1, 5,10, 6 }, + { 2, 6,11, 7 }, + { 3, 7,12, 8 }, + { 4, 9,11,10 }, + { 4, 9,12,11 }, + { 10, 5, 9, 8 }, + { 10, 8, 9,12 }, + { 10, 8,12, 7 }, + { 10, 7,12,11 }, + { 10, 7,11, 6 }, + { 10, 5, 8, 6 }, + { 10, 6, 8, 7 }, + + // quadratic pentahedron + { 12, 0, 8, 6 }, + { 12, 8, 7, 6 }, + { 12, 8, 2, 7 }, + { 12, 6, 7, 1 }, + { 12, 1, 7,13 }, + { 12, 7, 2,13 }, + { 12, 2,14,13 }, + + { 12, 3, 9,11 }, + { 12,11, 9,10 }, + { 12,11,10, 5 }, + { 12, 9, 4,10 }, + { 12,14, 5,10 }, + { 12,14,10, 4 }, + { 12,14, 4,13 }, + + // quadratic hexahedron + { 16, 0,11, 8 }, + { 16,11, 9, 8 }, + { 16, 8, 9, 1 }, + { 16,11, 3,10 }, + { 16,11,10, 9 }, + { 16,10, 2, 9 }, + { 16, 3,19, 2 }, + { 16, 2,19,18 }, + { 16, 2,18,17 }, + { 16, 2,17, 1 }, + + { 16, 4,12,15 }, + { 16,12, 5,13 }, + { 16,12,13,15 }, + { 16,13, 6,14 }, + { 16,13,14,15 }, + { 16,14, 7,15 }, + { 16, 6, 5,17 }, + { 16,18, 6,17 }, + { 16,18, 7, 6 }, + { 16,18,19, 7 }, + }; int type = GetVolumeType(); int n1 = ind[type]; int n2 = ind[type+1]; - for (int i = n1; i < n2; i++) - { + for (int i = n1; i < n2; i++) { V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ], myVolumeNodes[ vtab[i][1] ], myVolumeNodes[ vtab[i][2] ], @@ -529,6 +856,32 @@ bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const return true; } +//================================================================================ +/*! + * \brief Classify a point + * \param tol - thickness of faces + */ +//================================================================================ + +bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const +{ + // LIMITATION: for convex volumes only + XYZ p( X,Y,Z ); + for ( int iF = 0; iF < myNbFaces; ++iF ) + { + XYZ faceNormal; + if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z )) + continue; + if ( !IsFaceExternal( iF )) + faceNormal = XYZ() - faceNormal; // reverse + + XYZ face2p( p - XYZ( myFaceNodes[0] )); + if ( face2p.Dot( faceNormal ) > tol ) + return true; + } + return false; +} + //======================================================================= //function : SetExternalNormal //purpose : Node order will be so that faces normals are external @@ -545,7 +898,7 @@ void SMDS_VolumeTool::SetExternalNormal () //purpose : Return number of nodes in the array of face nodes //======================================================================= -int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) +int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const { if ( !setFace( faceIndex )) return 0; @@ -560,7 +913,7 @@ int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) // the last node == the first one. //======================================================================= -const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) +const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const { if ( !setFace( faceIndex )) return 0; @@ -575,14 +928,19 @@ const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) // the last node index == the first one. //======================================================================= -const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) +const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const { - if (myVolume->IsPoly()) { - MESSAGE("Warning: attempt to obtain FaceNodesIndices of polyhedral volume"); - return NULL; - } if ( !setFace( faceIndex )) return 0; + + if (myPolyedre) + { + SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this ); + me->myPolyIndices.resize( myFaceNbNodes + 1 ); + me->myFaceNodeIndices = & me->myPolyIndices[0]; + for ( int i = 0; i <= myFaceNbNodes; ++i ) + me->myFaceNodeIndices[i] = myVolume->GetNodeIndex( myFaceNodes[i] ); + } return myFaceNodeIndices; } @@ -592,25 +950,23 @@ const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) //======================================================================= bool SMDS_VolumeTool::GetFaceNodes (int faceIndex, - set& theFaceNodes ) + set& theFaceNodes ) const { if ( !setFace( faceIndex )) return false; theFaceNodes.clear(); - int iNode, nbNode = myFaceNbNodes; - for ( iNode = 0; iNode < nbNode; iNode++ ) - theFaceNodes.insert( myFaceNodes[ iNode ]); + theFaceNodes.insert( myFaceNodes, myFaceNodes + myFaceNbNodes ); return true; } //======================================================================= //function : IsFaceExternal -//purpose : Check normal orientation of a returned face +//purpose : Check normal orientation of a given face //======================================================================= -bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) +bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const { if ( myExternalFaces || !myVolume ) return true; @@ -625,22 +981,29 @@ bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) return true; } - switch ( myVolumeNbNodes ) { - case 4: - case 5: - // only the bottom of a reversed tetrahedron can be internal - return ( myVolForward || faceIndex != 0 ); - case 6: - // in a forward pentahedron, the top is internal, in a reversed one - bottom - return ( myVolForward ? faceIndex != 1 : faceIndex != 0 ); - case 8: { - // in a forward hexahedron, even face normal is external, odd - internal - bool odd = faceIndex % 2; - return ( myVolForward ? !odd : odd ); - } - default:; - } - return false; + // switch ( myVolumeNbNodes ) { + // case 4: + // case 5: + // case 10: + // case 13: + // // only the bottom of a reversed tetrahedron can be internal + // return ( myVolForward || faceIndex != 0 ); + // case 6: + // case 15: + // case 12: + // // in a forward prism, the top is internal, in a reversed one - bottom + // return ( myVolForward ? faceIndex != 1 : faceIndex != 0 ); + // case 8: + // case 20: + // case 27: { + // // in a forward hexahedron, even face normal is external, odd - internal + // bool odd = faceIndex % 2; + // return ( myVolForward ? !odd : odd ); + // } + // default:; + // } + // return false; + return true; } //======================================================================= @@ -648,29 +1011,28 @@ bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) //purpose : Return a normal to a face //======================================================================= -bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) +bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) const { if ( !setFace( faceIndex )) return false; - XYZ p1 ( myFaceNodes[0] ); - XYZ p2 ( myFaceNodes[1] ); - XYZ p3 ( myFaceNodes[2] ); + const int iQuad = ( myFaceNbNodes > 6 && !myPolyedre ) ? 2 : 1; + XYZ p1 ( myFaceNodes[0*iQuad] ); + XYZ p2 ( myFaceNodes[1*iQuad] ); + XYZ p3 ( myFaceNodes[2*iQuad] ); XYZ aVec12( p2 - p1 ); XYZ aVec13( p3 - p1 ); XYZ cross = aVec12.Crossed( aVec13 ); - if ( myFaceNbNodes == 4 ) { - XYZ p4 ( myFaceNodes[3] ); + if ( myFaceNbNodes >3*iQuad ) { + XYZ p4 ( myFaceNodes[3*iQuad] ); XYZ aVec14( p4 - p1 ); XYZ cross2 = aVec13.Crossed( aVec14 ); - cross.x += cross2.x; - cross.y += cross2.y; - cross.z += cross2.z; + cross = cross + cross2; } double size = cross.Magnitude(); - if ( size <= DBL_MIN ) + if ( size <= numeric_limits::min() ) return false; X = cross.x / size; @@ -680,12 +1042,33 @@ bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, doub return true; } +//================================================================================ +/*! + * \brief Return barycenter of a face + */ +//================================================================================ + +bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const +{ + if ( !setFace( faceIndex )) + return false; + + X = Y = Z = 0.0; + for ( int i = 0; i < myFaceNbNodes; ++i ) + { + X += myFaceNodes[i]->X() / myFaceNbNodes; + Y += myFaceNodes[i]->Y() / myFaceNbNodes; + Z += myFaceNodes[i]->Z() / myFaceNbNodes; + } + return true; +} + //======================================================================= //function : GetFaceArea //purpose : Return face area //======================================================================= -double SMDS_VolumeTool::GetFaceArea( int faceIndex ) +double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const { if (myVolume->IsPoly()) { MESSAGE("Warning: attempt to obtain area of a face of polyhedral volume"); @@ -710,6 +1093,26 @@ double SMDS_VolumeTool::GetFaceArea( int faceIndex ) return area; } +//================================================================================ +/*! + * \brief Return index of the node located at face center of a quadratic element like HEX27 + */ +//================================================================================ + +int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const +{ + if ( myAllFacesNbNodes && myVolumeNbNodes == 27 ) // classic element with 27 nodes + { + switch ( faceIndex ) { + case 0: return 20; + case 1: return 25; + default: + return faceIndex + 19; + } + } + return -1; +} + //======================================================================= //function : GetOppFaceIndex //purpose : Return index of the opposite face if it exists, else -1. @@ -718,19 +1121,32 @@ double SMDS_VolumeTool::GetFaceArea( int faceIndex ) int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const { int ind = -1; - if (myVolume->IsPoly()) { + if (myPolyedre) { MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume"); return ind; } + const int nbHoriFaces = 2; + if ( faceIndex >= 0 && faceIndex < NbFaces() ) { switch ( myVolumeNbNodes ) { case 6: + case 15: if ( faceIndex == 0 || faceIndex == 1 ) ind = 1 - faceIndex; break; case 8: - ind = faceIndex + ( faceIndex % 2 ? -1 : 1 ); + case 12: + if ( faceIndex <= 1 ) // top or bottom + ind = 1 - faceIndex; + else { + const int nbSideFaces = myAllFacesNbNodes[0]; + ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces; + } + break; + case 20: + case 27: + ind = GetOppFaceIndexOfHex( faceIndex ); break; default:; } @@ -738,13 +1154,25 @@ int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const return ind; } +//======================================================================= +//function : GetOppFaceIndexOfHex +//purpose : Return index of the opposite face of the hexahedron +//======================================================================= + +int SMDS_VolumeTool::GetOppFaceIndexOfHex( int faceIndex ) +{ + return Hexa_oppF[ faceIndex ]; +} + //======================================================================= //function : IsLinked //purpose : return true if theNode1 is linked with theNode2 +// If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well //======================================================================= bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1, - const SMDS_MeshNode* theNode2) const + const SMDS_MeshNode* theNode2, + const bool theIgnoreMediumNodes) const { if ( !myVolume ) return false; @@ -777,12 +1205,13 @@ bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1, } // find nodes indices - int i1 = -1, i2 = -1; - for ( int i = 0; i < myVolumeNbNodes; i++ ) { + int i1 = -1, i2 = -1, nbFound = 0; + for ( int i = 0; i < myVolumeNbNodes && nbFound < 2; i++ ) + { if ( myVolumeNodes[ i ] == theNode1 ) - i1 = i; + i1 = i, ++nbFound; else if ( myVolumeNodes[ i ] == theNode2 ) - i2 = i; + i2 = i, ++nbFound; } return IsLinked( i1, i2 ); } @@ -791,25 +1220,50 @@ bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1, //function : IsLinked //purpose : return true if the node with theNode1Index is linked // with the node with theNode2Index +// If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well //======================================================================= bool SMDS_VolumeTool::IsLinked (const int theNode1Index, - const int theNode2Index) const + const int theNode2Index, + bool theIgnoreMediumNodes) const { - if (myVolume->IsPoly()) { + if ( myVolume->IsPoly() ) { return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]); } - int minInd = theNode1Index < theNode2Index ? theNode1Index : theNode2Index; - int maxInd = theNode1Index < theNode2Index ? theNode2Index : theNode1Index; + int minInd = min( theNode1Index, theNode2Index ); + int maxInd = max( theNode1Index, theNode2Index ); if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd ) return false; - switch ( myVolumeNbNodes ) { - case 4: + VolumeType type = GetVolumeType(); + if ( myVolume->IsQuadratic() ) + { + int firstMediumInd = myVolume->NbCornerNodes(); + if ( minInd >= firstMediumInd ) + return false; // both nodes are medium - not linked + if ( maxInd < firstMediumInd ) // both nodes are corners + { + if ( theIgnoreMediumNodes ) + type = quadToLinear(type); // to check linkage of corner nodes only + else + return false; // corner nodes are not linked directly in a quadratic cell + } + } + + switch ( type ) { + case TETRA: return true; - case 5: + case HEXA: + switch ( maxInd - minInd ) { + case 1: return minInd != 3; + case 3: return minInd == 0 || minInd == 4; + case 4: return true; + default:; + } + break; + case PYRAM: if ( maxInd == 4 ) return true; switch ( maxInd - minInd ) { @@ -818,7 +1272,7 @@ bool SMDS_VolumeTool::IsLinked (const int theNode1Index, default:; } break; - case 6: + case PENTA: switch ( maxInd - minInd ) { case 1: return minInd != 2; case 2: return minInd == 0 || minInd == 3; @@ -826,14 +1280,64 @@ bool SMDS_VolumeTool::IsLinked (const int theNode1Index, default:; } break; - case 8: - switch ( maxInd - minInd ) { - case 1: return minInd != 3; - case 3: return minInd == 0 || minInd == 4; - case 4: return true; - default:; + case QUAD_TETRA: + { + switch ( minInd ) { + case 0: if( maxInd==4 || maxInd==6 || maxInd==7 ) return true; + case 1: if( maxInd==4 || maxInd==5 || maxInd==8 ) return true; + case 2: if( maxInd==5 || maxInd==6 || maxInd==9 ) return true; + case 3: if( maxInd==7 || maxInd==8 || maxInd==9 ) return true; + default:; + } + break; + } + case QUAD_HEXA: + { + switch ( minInd ) { + case 0: if( maxInd==8 || maxInd==11 || maxInd==16 ) return true; + case 1: if( maxInd==8 || maxInd==9 || maxInd==17 ) return true; + case 2: if( maxInd==9 || maxInd==10 || maxInd==18 ) return true; + case 3: if( maxInd==10 || maxInd==11 || maxInd==19 ) return true; + case 4: if( maxInd==12 || maxInd==15 || maxInd==16 ) return true; + case 5: if( maxInd==12 || maxInd==13 || maxInd==17 ) return true; + case 6: if( maxInd==13 || maxInd==14 || maxInd==18 ) return true; + case 7: if( maxInd==14 || maxInd==15 || maxInd==19 ) return true; + default:; + } + break; + } + case QUAD_PYRAM: + { + switch ( minInd ) { + case 0: if( maxInd==5 || maxInd==8 || maxInd==9 ) return true; + case 1: if( maxInd==5 || maxInd==6 || maxInd==10 ) return true; + case 2: if( maxInd==6 || maxInd==7 || maxInd==11 ) return true; + case 3: if( maxInd==7 || maxInd==8 || maxInd==12 ) return true; + case 4: if( maxInd==9 || maxInd==10 || maxInd==11 || maxInd==12 ) return true; + default:; + } + break; + } + case QUAD_PENTA: + { + switch ( minInd ) { + case 0: if( maxInd==6 || maxInd==8 || maxInd==12 ) return true; + case 1: if( maxInd==6 || maxInd==7 || maxInd==13 ) return true; + case 2: if( maxInd==7 || maxInd==8 || maxInd==14 ) return true; + case 3: if( maxInd==9 || maxInd==11 || maxInd==12 ) return true; + case 4: if( maxInd==9 || maxInd==10 || maxInd==13 ) return true; + case 5: if( maxInd==10 || maxInd==11 || maxInd==14 ) return true; + default:; + } + break; + } + case HEX_PRISM: + { + const int diff = maxInd-minInd; + if ( diff > 6 ) return false;// not linked top and bottom + if ( diff == 6 ) return true; // linked top and bottom + return diff == 1 || diff == 7; } - break; default:; } return false; @@ -855,46 +1359,221 @@ int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const return -1; } -//======================================================================= -//function : IsFreeFace -//purpose : check that only one volume is build on the face nodes -//======================================================================= +//================================================================================ +/*! + * \brief Fill vector with boundary faces existing in the mesh + * \param faces - vector of found nodes + * \retval int - nb of found faces + */ +//================================================================================ + +int SMDS_VolumeTool::GetAllExistingFaces(vector & faces) const +{ + faces.clear(); + for ( int iF = 0; iF < NbFaces(); ++iF ) { + const SMDS_MeshFace* face = 0; + const SMDS_MeshNode** nodes = GetFaceNodes( iF ); + switch ( NbFaceNodes( iF )) { + case 3: + face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break; + case 4: + face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break; + case 6: + face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], + nodes[3], nodes[4], nodes[5]); break; + case 8: + face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3], + nodes[4], nodes[5], nodes[6], nodes[7]); break; + } + if ( face ) + faces.push_back( face ); + } + return faces.size(); +} + + +//================================================================================ +/*! + * \brief Fill vector with boundary edges existing in the mesh + * \param edges - vector of found edges + * \retval int - nb of found faces + */ +//================================================================================ + +int SMDS_VolumeTool::GetAllExistingEdges(vector & edges) const +{ + edges.clear(); + edges.reserve( myVolumeNbNodes * 2 ); + for ( int i = 0; i < myVolumeNbNodes-1; ++i ) { + for ( int j = i + 1; j < myVolumeNbNodes; ++j ) { + if ( IsLinked( i, j )) { + const SMDS_MeshElement* edge = + SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] ); + if ( edge ) + edges.push_back( edge ); + } + } + } + return edges.size(); +} + +//================================================================================ +/*! + * \brief Return minimal square distance between connected corner nodes + */ +//================================================================================ + +double SMDS_VolumeTool::MinLinearSize2() const +{ + double minSize = 1e+100; + int iQ = myVolume->IsQuadratic() ? 2 : 1; + + // store current face data + int curFace = myCurFace, nbN = myFaceNbNodes; + int* ind = myFaceNodeIndices; + myFaceNodeIndices = NULL; + const SMDS_MeshNode** nodes = myFaceNodes; + myFaceNodes = NULL; + + // it seems that compute distance twice is faster than organization of a sole computing + myCurFace = -1; + for ( int iF = 0; iF < myNbFaces; ++iF ) + { + setFace( iF ); + for ( int iN = 0; iN < myFaceNbNodes; iN += iQ ) + { + XYZ n1( myFaceNodes[ iN ]); + XYZ n2( myFaceNodes[(iN + iQ) % myFaceNbNodes]); + minSize = std::min( minSize, (n1 - n2).SquareMagnitude()); + } + } + // restore current face data + myCurFace = curFace; + myFaceNbNodes = nbN; + myFaceNodeIndices = ind; + delete [] myFaceNodes; myFaceNodes = nodes; + + return minSize; +} + +//================================================================================ +/*! + * \brief Return maximal square distance between connected corner nodes + */ +//================================================================================ + +double SMDS_VolumeTool::MaxLinearSize2() const +{ + double maxSize = -1e+100; + int iQ = myVolume->IsQuadratic() ? 2 : 1; + + // store current face data + int curFace = myCurFace, nbN = myFaceNbNodes; + int* ind = myFaceNodeIndices; + myFaceNodeIndices = NULL; + const SMDS_MeshNode** nodes = myFaceNodes; + myFaceNodes = NULL; + + // it seems that compute distance twice is faster than organization of a sole computing + myCurFace = -1; + for ( int iF = 0; iF < myNbFaces; ++iF ) + { + setFace( iF ); + for ( int iN = 0; iN < myFaceNbNodes; iN += iQ ) + { + XYZ n1( myFaceNodes[ iN ]); + XYZ n2( myFaceNodes[(iN + iQ) % myFaceNbNodes]); + maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude()); + } + } + // restore current face data + myCurFace = curFace; + myFaceNbNodes = nbN; + myFaceNodeIndices = ind; + delete [] myFaceNodes; myFaceNodes = nodes; + + return maxSize; +} + +//================================================================================ +/*! + * \brief fast check that only one volume is build on the face nodes + */ +//================================================================================ + +bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const +{ + const bool isFree = true; + + if (!setFace( faceIndex )) + return !isFree; + + const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex ); + + // a set of facet nodes w/o medium ones and w/o nodes[0] + set< const SMDS_MeshElement* > nodeSet; + const int di = myVolume->IsQuadratic() ? 2 : 1; + for ( int i = di; i < myFaceNbNodes; i += di ) + nodeSet.insert( nodes[i] ); + + SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume ); + SMDS_ElemIteratorPtr nIt; + while ( eIt->more() ) { + const SMDS_MeshElement* vol = eIt->next(); + if ( vol != myVolume ) { + size_t nbShared = 0; + if ( const SMDS_VtkVolume* v = dynamic_cast< const SMDS_VtkVolume* >( vol )) + nIt = v->uniqueNodesIterator(); + else + nIt = vol->nodesIterator(); + while ( nIt->more() ) + if (( nbShared += nodeSet.count( nIt->next() )) == nodeSet.size() ) + { + if ( otherVol ) *otherVol = vol; + return !isFree; + } + } + } + if ( otherVol ) *otherVol = 0; + return isFree; +} + +//================================================================================ +/*! + * \brief Thorough check that only one volume is build on the face nodes + */ +//================================================================================ -bool SMDS_VolumeTool::IsFreeFace( int faceIndex ) +bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const { - const int free = true; + const bool isFree = true; if (!setFace( faceIndex )) - return !free; + return !isFree; const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex ); - int nbFaceNodes = myFaceNbNodes; + const int nbFaceNodes = myFaceNbNodes; - // evaluate nb of face nodes shared by other volume + // evaluate nb of face nodes shared by other volumes int maxNbShared = -1; typedef map< const SMDS_MeshElement*, int > TElemIntMap; TElemIntMap volNbShared; TElemIntMap::iterator vNbIt; - for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) - { + for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) { const SMDS_MeshNode* n = nodes[ iNode ]; - SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator(); + SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume ); while ( eIt->more() ) { const SMDS_MeshElement* elem = eIt->next(); - if ( elem != myVolume && elem->GetType() == SMDSAbs_Volume ) { - int nbShared = 1; - vNbIt = volNbShared.find( elem ); - if ( vNbIt == volNbShared.end() ) - volNbShared.insert ( TElemIntMap::value_type( elem, nbShared )); - else - nbShared = ++(*vNbIt).second; - if ( nbShared > maxNbShared ) - maxNbShared = nbShared; + if ( elem != myVolume ) { + vNbIt = volNbShared.insert( make_pair( elem, 0 )).first; + (*vNbIt).second++; + if ( vNbIt->second > maxNbShared ) + maxNbShared = vNbIt->second; } } } if ( maxNbShared < 3 ) - return free; // is free + return isFree; // is free // find volumes laying on the opposite side of the face // and sharing all nodes @@ -903,56 +1582,81 @@ bool SMDS_VolumeTool::IsFreeFace( int faceIndex ) if ( IsFaceExternal( faceIndex )) intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z ); XYZ p0 ( nodes[0] ), baryCenter; - for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) - { - int nbShared = (*vNbIt).second; + for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); ) { + const int& nbShared = (*vNbIt).second; if ( nbShared >= 3 ) { SMDS_VolumeTool volume( (*vNbIt).first ); volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z ); XYZ intNormal2( baryCenter - p0 ); - if ( intNormal.Dot( intNormal2 ) < 0 ) - continue; // opposite side + if ( intNormal.Dot( intNormal2 ) < 0 ) { + // opposite side + if ( nbShared >= nbFaceNodes ) + { + // a volume shares the whole facet + if ( otherVol ) *otherVol = vNbIt->first; + return !isFree; + } + ++vNbIt; + continue; + } } // remove a volume from volNbShared map - volNbShared.erase( vNbIt ); + volNbShared.erase( vNbIt++ ); + } + + // here volNbShared contains only volumes laying on the opposite side of + // the face and sharing 3 or more but not all face nodes with myVolume + if ( volNbShared.size() < 2 ) { + return isFree; // is free } - // here volNbShared contains only volumes laying on the - // opposite side of the face - if ( volNbShared.empty() ) - return free; // is free // check if the whole area of a face is shared - bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle - for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) + for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) { - SMDS_VolumeTool volume( (*vNbIt).first ); - bool prevLinkShared = false; - int nbSharedLinks = 0; - for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) - { - bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] ); - if ( linkShared ) - nbSharedLinks++; - if ( linkShared && prevLinkShared && - volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] )) - isShared[ iNode ] = true; - prevLinkShared = linkShared; - } - if ( nbSharedLinks == nbFaceNodes ) - return !free; // is not free - if ( nbFaceNodes == 4 ) { - // check traingle parts 1 & 3 - if ( isShared[1] && isShared[3] ) - return !free; // is not free - // check triangle parts 0 & 2; - // 0 part could not be checked in the loop; check it here - if ( isShared[2] && prevLinkShared && - volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) && - volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) ) - return !free; // is not free - } + const SMDS_MeshNode* n = nodes[ iNode ]; + // check if n is shared by one of volumes of volNbShared + bool isShared = false; + SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume ); + while ( eIt->more() && !isShared ) + isShared = volNbShared.count( eIt->next() ); + if ( !isShared ) + return isFree; } - return free; + if ( otherVol ) *otherVol = volNbShared.begin()->first; + return !isFree; + +// if ( !myVolume->IsPoly() ) +// { +// bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle +// for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) { +// SMDS_VolumeTool volume( (*vNbIt).first ); +// bool prevLinkShared = false; +// int nbSharedLinks = 0; +// for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) { +// bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] ); +// if ( linkShared ) +// nbSharedLinks++; +// if ( linkShared && prevLinkShared && +// volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] )) +// isShared[ iNode ] = true; +// prevLinkShared = linkShared; +// } +// if ( nbSharedLinks == nbFaceNodes ) +// return !free; // is not free +// if ( nbFaceNodes == 4 ) { +// // check traingle parts 1 & 3 +// if ( isShared[1] && isShared[3] ) +// return !free; // is not free +// // check triangle parts 0 & 2; +// // 0 part could not be checked in the loop; check it here +// if ( isShared[2] && prevLinkShared && +// volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) && +// volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) ) +// return !free; // is not free +// } +// } +// } +// return free; } //======================================================================= @@ -960,16 +1664,40 @@ bool SMDS_VolumeTool::IsFreeFace( int faceIndex ) //purpose : Return index of a face formed by theFaceNodes //======================================================================= -int SMDS_VolumeTool::GetFaceIndex( const set& theFaceNodes ) +int SMDS_VolumeTool::GetFaceIndex( const set& theFaceNodes, + const int theFaceIndexHint ) const { - for ( int iFace = 0; iFace < myNbFaces; iFace++ ) { - const SMDS_MeshNode** nodes = GetFaceNodes( iFace ); - int nbFaceNodes = NbFaceNodes( iFace ); - set nodeSet; - for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) - nodeSet.insert( nodes[ iNode ] ); - if ( theFaceNodes == nodeSet ) - return iFace; + if ( theFaceIndexHint >= 0 ) + { + int nbNodes = NbFaceNodes( theFaceIndexHint ); + if ( nbNodes == (int) theFaceNodes.size() ) + { + const SMDS_MeshNode** nodes = GetFaceNodes( theFaceIndexHint ); + while ( nbNodes ) + if ( theFaceNodes.count( nodes[ nbNodes-1 ])) + --nbNodes; + else + break; + if ( nbNodes == 0 ) + return theFaceIndexHint; + } + } + for ( int iFace = 0; iFace < myNbFaces; iFace++ ) + { + if ( iFace == theFaceIndexHint ) + continue; + int nbNodes = NbFaceNodes( iFace ); + if ( nbNodes == (int) theFaceNodes.size() ) + { + const SMDS_MeshNode** nodes = GetFaceNodes( iFace ); + while ( nbNodes ) + if ( theFaceNodes.count( nodes[ nbNodes-1 ])) + --nbNodes; + else + break; + if ( nbNodes == 0 ) + return iFace; + } } return -1; } @@ -998,7 +1726,7 @@ int SMDS_VolumeTool::GetFaceIndex( const set& theFaceNodes //purpose : //======================================================================= -bool SMDS_VolumeTool::setFace( int faceIndex ) +bool SMDS_VolumeTool::setFace( int faceIndex ) const { if ( !myVolume ) return false; @@ -1016,64 +1744,120 @@ bool SMDS_VolumeTool::setFace( int faceIndex ) myFaceNodes = NULL; } - if (myVolume->IsPoly()) { + if (myVolume->IsPoly()) + { if (!myPolyedre) { MESSAGE("Warning: bad volumic element"); return false; } - // check orientation - bool isGoodOri = true; - if (myExternalFaces) - isGoodOri = IsFaceExternal( faceIndex ); - // set face nodes int iNode; myFaceNbNodes = myPolyedre->NbFaceNodes(faceIndex + 1); myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1]; - if (isGoodOri) { - for ( iNode = 0; iNode < myFaceNbNodes; iNode++ ) - myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1); - } else { - for ( iNode = 0; iNode < myFaceNbNodes; iNode++ ) - myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, myFaceNbNodes - iNode); - } + for ( iNode = 0; iNode < myFaceNbNodes; iNode++ ) + myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1); myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first - } else { - // choose face node indices - switch ( myVolumeNbNodes ) { - case 4: - myFaceNbNodes = Tetra_nbN[ faceIndex ]; - if ( myExternalFaces ) - myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_RE[ faceIndex ]; - else - myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_R[ faceIndex ]; - break; - case 5: - myFaceNbNodes = Pyramid_nbN[ faceIndex ]; - if ( myExternalFaces ) - myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_RE[ faceIndex ]; - else - myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_R[ faceIndex ]; - break; - case 6: - myFaceNbNodes = Penta_nbN[ faceIndex ]; - if ( myExternalFaces ) - myFaceNodeIndices = myVolForward ? Penta_FE[ faceIndex ] : Penta_RE[ faceIndex ]; - else - myFaceNodeIndices = myVolForward ? Penta_F[ faceIndex ] : Penta_R[ faceIndex ]; - break; - case 8: - myFaceNbNodes = Hexa_nbN[ faceIndex ]; - if ( myExternalFaces ) - myFaceNodeIndices = myVolForward ? Hexa_FE[ faceIndex ] : Hexa_RE[ faceIndex ]; - else - myFaceNodeIndices = Hexa_F[ faceIndex ]; - break; - default: - return false; + // check orientation + if (myExternalFaces) + { + myCurFace = faceIndex; // avoid infinite recursion in IsFaceExternal() + myExternalFaces = false; // force normal computation by IsFaceExternal() + if ( !IsFaceExternal( faceIndex )) + for ( int i = 0, j = myFaceNbNodes; i < j; ++i, --j ) + std::swap( myFaceNodes[i], myFaceNodes[j] ); + myExternalFaces = true; + } + } + else + { + if ( !myAllFacesNodeIndices_F ) + { + // choose data for an element type + switch ( myVolumeNbNodes ) { + case 4: + myAllFacesNodeIndices_F = &Tetra_F [0][0]; + //myAllFacesNodeIndices_FE = &Tetra_F [0][0]; + myAllFacesNodeIndices_RE = &Tetra_RE[0][0]; + myAllFacesNbNodes = Tetra_nbN; + myMaxFaceNbNodes = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]); + break; + case 5: + myAllFacesNodeIndices_F = &Pyramid_F [0][0]; + //myAllFacesNodeIndices_FE = &Pyramid_F [0][0]; + myAllFacesNodeIndices_RE = &Pyramid_RE[0][0]; + myAllFacesNbNodes = Pyramid_nbN; + myMaxFaceNbNodes = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]); + break; + case 6: + myAllFacesNodeIndices_F = &Penta_F [0][0]; + //myAllFacesNodeIndices_FE = &Penta_FE[0][0]; + myAllFacesNodeIndices_RE = &Penta_RE[0][0]; + myAllFacesNbNodes = Penta_nbN; + myMaxFaceNbNodes = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]); + break; + case 8: + myAllFacesNodeIndices_F = &Hexa_F [0][0]; + ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0]; + myAllFacesNodeIndices_RE = &Hexa_RE[0][0]; + myAllFacesNbNodes = Hexa_nbN; + myMaxFaceNbNodes = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]); + break; + case 10: + myAllFacesNodeIndices_F = &QuadTetra_F [0][0]; + //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0]; + myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0]; + myAllFacesNbNodes = QuadTetra_nbN; + myMaxFaceNbNodes = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]); + break; + case 13: + myAllFacesNodeIndices_F = &QuadPyram_F [0][0]; + //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0]; + myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0]; + myAllFacesNbNodes = QuadPyram_nbN; + myMaxFaceNbNodes = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]); + break; + case 15: + myAllFacesNodeIndices_F = &QuadPenta_F [0][0]; + //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0]; + myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0]; + myAllFacesNbNodes = QuadPenta_nbN; + myMaxFaceNbNodes = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]); + break; + case 20: + case 27: + myAllFacesNodeIndices_F = &QuadHexa_F [0][0]; + //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0]; + myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0]; + myAllFacesNbNodes = QuadHexa_nbN; + myMaxFaceNbNodes = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]); + if ( !myIgnoreCentralNodes && myVolumeNbNodes == 27 ) + { + myAllFacesNodeIndices_F = &TriQuadHexa_F [0][0]; + //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0]; + myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0]; + myAllFacesNbNodes = TriQuadHexa_nbN; + myMaxFaceNbNodes = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]); + } + break; + case 12: + myAllFacesNodeIndices_F = &HexPrism_F [0][0]; + //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0]; + myAllFacesNodeIndices_RE = &HexPrism_RE[0][0]; + myAllFacesNbNodes = HexPrism_nbN; + myMaxFaceNbNodes = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]); + break; + default: + return false; + } } + myFaceNbNodes = myAllFacesNbNodes[ faceIndex ]; + // if ( myExternalFaces ) + // myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes ); + // else + // myFaceNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes ); + myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes ); // set face nodes myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1]; @@ -1099,6 +1883,12 @@ SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes) case 5: return PYRAM; case 6: return PENTA; case 8: return HEXA; + case 10: return QUAD_TETRA; + case 13: return QUAD_PYRAM; + case 15: return QUAD_PENTA; + case 20: + case 27: return QUAD_HEXA; + case 12: return HEX_PRISM; default:return UNKNOWN; } } @@ -1111,13 +1901,44 @@ SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes) int SMDS_VolumeTool::NbFaces( VolumeType type ) { switch ( type ) { - case TETRA: return 4; - case PYRAM: return 5; - case PENTA: return 5; - case HEXA : return 6; - default: return 0; + case TETRA : + case QUAD_TETRA: return 4; + case PYRAM : + case QUAD_PYRAM: return 5; + case PENTA : + case QUAD_PENTA: return 5; + case HEXA : + case QUAD_HEXA : return 6; + case HEX_PRISM : return 8; + default: return 0; + } +} + +//================================================================================ +/*! + * \brief Useful to know nb of corner nodes of a quadratic volume + * \param type - volume type + * \retval int - nb of corner nodes + */ +//================================================================================ + +int SMDS_VolumeTool::NbCornerNodes(VolumeType type) +{ + switch ( type ) { + case TETRA : + case QUAD_TETRA: return 4; + case PYRAM : + case QUAD_PYRAM: return 5; + case PENTA : + case QUAD_PENTA: return 6; + case HEXA : + case QUAD_HEXA : return 8; + case HEX_PRISM : return 12; + default: return 0; } + return 0; } + // //======================================================================= //function : GetFaceNodesIndices @@ -1134,8 +1955,14 @@ const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type, switch ( type ) { case TETRA: return Tetra_F[ faceIndex ]; case PYRAM: return Pyramid_F[ faceIndex ]; - case PENTA: return external ? Penta_FE[ faceIndex ] : Penta_F[ faceIndex ]; - case HEXA: return external ? Hexa_FE[ faceIndex ] : Hexa_F[ faceIndex ]; + case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ]; + case HEXA: return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ]; + case QUAD_TETRA: return QuadTetra_F[ faceIndex ]; + case QUAD_PYRAM: return QuadPyram_F[ faceIndex ]; + case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ]; + // what about SMDSEntity_TriQuad_Hexa? + case QUAD_HEXA: return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ]; + case HEX_PRISM: return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ]; default:; } return 0; @@ -1154,8 +1981,33 @@ int SMDS_VolumeTool::NbFaceNodes(VolumeType type, case PYRAM: return Pyramid_nbN[ faceIndex ]; case PENTA: return Penta_nbN[ faceIndex ]; case HEXA: return Hexa_nbN[ faceIndex ]; + case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ]; + case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ]; + case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ]; + // what about SMDSEntity_TriQuad_Hexa? + case QUAD_HEXA: return QuadHexa_nbN[ faceIndex ]; + case HEX_PRISM: return HexPrism_nbN[ faceIndex ]; default:; } return 0; } +//======================================================================= +//function : Element +//purpose : return element +//======================================================================= + +const SMDS_MeshVolume* SMDS_VolumeTool::Element() const +{ + return static_cast( myVolume ); +} + +//======================================================================= +//function : ID +//purpose : return element ID +//======================================================================= + +int SMDS_VolumeTool::ID() const +{ + return myVolume ? myVolume->GetID() : 0; +}