X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMDS%2FSMDS_VolumeTool.cxx;h=ef7c8422265399d952b3ececf608d44585a31cf8;hp=fe3c56adc05b583ad86211243b1c73d0c81cdc42;hb=4b5fd249abfca52265832f10435f18f5bb9c69ae;hpb=5a053fbfd04d90fea6aae6d14cb15260d5aaf5a7 diff --git a/src/SMDS/SMDS_VolumeTool.cxx b/src/SMDS/SMDS_VolumeTool.cxx index fe3c56adc..ef7c84222 100644 --- a/src/SMDS/SMDS_VolumeTool.cxx +++ b/src/SMDS/SMDS_VolumeTool.cxx @@ -1,8 +1,29 @@ +// Copyright (C) 2007-2011 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,7 +32,8 @@ #include "SMDS_MeshElement.hxx" #include "SMDS_MeshNode.hxx" -#include "SMDS_PolyhedralVolumeOfNodes.hxx" +#include "SMDS_VtkVolume.hxx" +#include "SMDS_Mesh.hxx" #include "utilities.h" @@ -166,9 +188,189 @@ static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL { 1, 2, 6, 5, 1 }}; static int Hexa_nbN [] = { 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 == EXTERNAL + { 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_R [4][7] = { // REVERSED + { 0, 4, 1, 5, 2, 6, 0 }, // All faces but a bottom 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_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 == EXTERNAL + { 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_R [5][9] = { // REVERSED + { 0, 5, 1, 6, 2, 7, 3, 8, 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_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 }, // Top face has an internal normal, other - external + { 3, 9, 4, 10,5, 11,3, 3, 3 }, // 0 is bottom, 1 is top face + { 0, 8, 2, 14,5, 11,3, 12,0 }, + { 1, 13,4, 10,5, 14,2, 7, 1 }, + { 0, 12,3, 9, 4, 13,1, 6, 0 }}; +static int QuadPenta_R [5][9] = { // REVERSED + { 0, 6, 1, 7, 2, 8, 0, 0, 0 }, // Bottom face has an internal normal, other - external + { 3, 9, 4, 10,5, 11,3, 3, 3 }, // 0 is bottom, 1 is top face + { 0, 12,3, 11,5, 14,2, 8, 0 }, + { 1, 7, 2, 14,5, 10,4, 13,1 }, + { 0, 6, 1, 13,4, 9, 3, 12,0 }}; +static int QuadPenta_FE [5][9] = { // FORWARD -> EXTERNAL + { 0, 6, 1, 7, 2, 8, 0, 0, 0 }, + { 3,11, 5, 10,4, 9, 3, 3, 3 }, + { 0, 8, 2, 14,5, 11,3, 12,0 }, + { 1, 13,4, 10,5, 14,2, 7, 1 }, + { 0, 12,3, 9, 4, 13,1, 6, 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, 12,3, 11,5, 14,2, 8, 0 }, + { 1, 7, 2, 14,5, 10,4, 13,1 }, + { 0, 6, 1, 13,4, 9, 3, 12,0 }}; +static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 }; + +/* +// 13 +// N5+-----+-----+N6 +// /| /| +// 12+ | 14+ | +// / | / | +// N4+-----+-----+N7 | QUADRATIC +// | | 15 | | HEXAHEDRON +// | | | | +// | 17+ | +18 +// | | | | +// | | | | +// | | | | +// 16+ | +19 | +// | | | | +// | | 9 | | +// | N1+-----+-|---+N2 +// | / | / +// | +8 | +10 +// |/ |/ +// N0+-----+-----+N3 +// 11 +*/ +static int QuadHexa_F [6][9] = { // FORWARD + { 0, 8, 1, 9, 2, 10,3, 11,0 }, // opposite faces are neighbouring, + { 4, 12,5, 13,6, 14,7, 15,4 }, // odd face(1,3,5) normal is internal, even(0,2,4) - external + { 1, 8, 0, 16,4, 12,5, 17,1 }, // same index nodes of opposite faces are linked + { 2, 10,3, 19,7, 14,6, 18,2 }, + { 0, 11,3, 19,7, 15,4, 16,0 }, + { 1, 9, 2, 18,6, 13,5, 17,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 QuadHexa_FE [6][9] = { // FORWARD -> EXTERNAL + { 0, 8, 1, 9, 2, 10,3, 11,0 }, // opposite faces are neighbouring, + { 4, 15,7, 14,6, 13,5, 12,4 }, // all face normals are external, + { 0, 16,4, 12,5, 17,1, 8, 0 }, // links in opposite faces: 0-0, 1-3, 2-2, 3-1 + { 3, 10,2, 18,6, 14,7, 19,3 }, + { 0, 11,3, 19,7, 15,4, 16,0 }, + { 1, 17,5, 13,6, 18,2, 9, 1 }}; +static int QuadHexa_RE [6][9] = { // REVERSED -> EXTERNAL + { 0, 11,3, 10,2, 9, 1, 8, 0 }, // opposite faces are neighbouring, + { 4, 12,5, 13,6, 14,7, 15,4 }, // all face normals are external, + { 0, 8, 1, 17,5, 12,4, 16,0 }, // links in opposite faces: 0-0, 1-3, 2-2, 3-1 + { 3, 19,7, 14,6, 18,2, 10,3 }, + { 0, 16,4, 15,7, 19,3, 11,0 }, + { 1, 9, 2, 18,6, 13,5, 17,1 }}; +static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 }; + + // ======================================================== // to perform some calculations without linkage to CASCADE // ======================================================== +namespace +{ struct XYZ { double x; double y; @@ -177,25 +379,26 @@ 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 Crossed( const XYZ& other ); + inline double Dot( const XYZ& other ); + inline double Magnitude(); }; -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::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); } +} //======================================================================= //function : SMDS_VolumeTool @@ -215,6 +418,7 @@ SMDS_VolumeTool::SMDS_VolumeTool () myFaceNodeIndices( NULL ), myFaceNodes( NULL ) { + //MESSAGE("******************************************************** SMDS_VolumeToo"); } //======================================================================= @@ -235,6 +439,7 @@ SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume) myFaceNodeIndices( NULL ), myFaceNodes( NULL ) { + //MESSAGE("******************************************************** SMDS_VolumeToo"); Set( theVolume ); } @@ -245,14 +450,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; } //======================================================================= @@ -299,22 +501,29 @@ bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume) } if (myVolume->IsPoly()) { - myPolyedre = static_cast( myVolume ); + myPolyedre = dynamic_cast( myVolume ); if (!myPolyedre) { MESSAGE("Warning: bad volumic element"); return false; } - } else { + } + else { switch ( myVolumeNbNodes ) { case 4: case 5: case 6: - case 8: { + case 8: + case 10: + case 13: + case 15: + case 20: { // 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 ]; + 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() ); @@ -368,10 +577,109 @@ void SMDS_VolumeTool::Inverse () SWAP_NODES( myVolumeNodes, 1, 3 ); SWAP_NODES( myVolumeNodes, 5, 7 ); 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; default:; } } +//======================================================================= +//function : GetVolumeType +//purpose : +//======================================================================= + +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; break; + case 5: return PYRAM; break; + case 6: return PENTA; break; + case 8: return HEXA; break; + case 10: return QUAD_TETRA; break; + case 13: return QUAD_PYRAM; break; + case 15: return QUAD_PENTA; break; + case 20: return QUAD_HEXA; break; + default: break; + } + } + + return UNKNOWN; +} + +//======================================================================= +//function : getTetraVolume +//purpose : +//======================================================================= + +static double getTetraVolume(const SMDS_MeshNode* n1, + const SMDS_MeshNode* n2, + const SMDS_MeshNode* n3, + const SMDS_MeshNode* n4) +{ + double X1 = n1->X(); + double Y1 = n1->Y(); + double Z1 = n1->Z(); + + double X2 = n2->X(); + double Y2 = n2->Y(); + double Z2 = n2->Z(); + + double X3 = n3->X(); + double Y3 = n3->Y(); + double Z3 = n3->Z(); + + double X4 = n4->X(); + double Y4 = n4->Y(); + double Z4 = n4->Z(); + + double Q1 = -(X1-X2)*(Y3*Z4-Y4*Z3); + double Q2 = (X1-X3)*(Y2*Z4-Y4*Z2); + double R1 = -(X1-X4)*(Y2*Z3-Y3*Z2); + double R2 = -(X2-X3)*(Y1*Z4-Y4*Z1); + double S1 = (X2-X4)*(Y1*Z3-Y3*Z1); + double S2 = -(X3-X4)*(Y1*Z2-Y2*Z1); + + return (Q1+Q2+R1+R2+S1+S2)/6.0; +} + //======================================================================= //function : GetSize //purpose : Return element volume @@ -379,7 +687,138 @@ void SMDS_VolumeTool::Inverse () double SMDS_VolumeTool::GetSize() const { - return 0; + double V = 0.; + if ( !myVolume ) + return 0.; + + if ( myVolume->IsPoly() ) + { + if ( !myPolyedre ) + return 0.; + + SMDS_Mesh *mesh = SMDS_Mesh::_meshList[myPolyedre->getMeshId()]; + // split a polyhedron into tetrahedrons + + SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this ); + XYZ baryCenter; + me->GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z); + SMDS_MeshNode *bcNode = mesh->AddNode( 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 ) + { + double Vn = getTetraVolume( myFaceNodes[ 0 ], + myFaceNodes[ n-1 ], + myFaceNodes[ n ], + bcNode ); +/// cout <<"++++ " << Vn << " nodes " <GetID() << " " <GetID() << " " <GetID() << " < " << V << endl; + V += externalFace ? -Vn : Vn; + } + } + mesh->RemoveNode(bcNode); + } + else + { + const static int ind[] = { + 0, 1, 3, 6, 11, 19, 32, 46, 66}; + const static int vtab[][4] = { + // tetrahedron + { 0, 1, 2, 3 }, + // pyramid + { 0, 1, 3, 4 }, + { 1, 2, 3, 4 }, + // pentahedron + { 0, 1, 2, 3 }, + { 1, 5, 3, 4 }, + { 1, 5, 2, 3 }, + // hexahedron + { 1, 4, 3, 0 }, + { 4, 1, 6, 5 }, + { 1, 3, 6, 2 }, + { 4, 6, 3, 7 }, + { 1, 4, 6, 3 }, + + // 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++) { + V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ], + myVolumeNodes[ vtab[i][1] ], + myVolumeNodes[ vtab[i][2] ], + myVolumeNodes[ vtab[i][3] ]); + } + } + return V; } //======================================================================= @@ -405,6 +844,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) +{ + // 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 @@ -453,12 +918,16 @@ const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) { - if (myVolume->IsPoly()) { - MESSAGE("Warning: attempt to obtain FaceNodesIndices of polyhedral volume"); - return NULL; - } if ( !setFace( faceIndex )) return 0; + + if (myVolume->IsPoly()) + { + myPolyIndices.resize( myFaceNbNodes + 1 ); + myFaceNodeIndices = & myPolyIndices[0]; + for ( int i = 0; i <= myFaceNbNodes; ++i ) + myFaceNodeIndices[i] = myVolume->GetNodeIndex( myFaceNodes[i] ); + } return myFaceNodeIndices; } @@ -504,12 +973,16 @@ bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) 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: // in a forward pentahedron, the top is internal, in a reversed one - bottom return ( myVolForward ? faceIndex != 1 : faceIndex != 0 ); - case 8: { + case 8: + case 20: { // in a forward hexahedron, even face normal is external, odd - internal bool odd = faceIndex % 2; return ( myVolForward ? !odd : odd ); @@ -536,7 +1009,8 @@ bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, doub XYZ aVec13( p3 - p1 ); XYZ cross = aVec12.Crossed( aVec13 ); - if ( myFaceNbNodes == 4 ) { + //if ( myFaceNbNodes == 4 ) { + if ( myFaceNbNodes >3 ) { XYZ p4 ( myFaceNodes[3] ); XYZ aVec14( p4 - p1 ); XYZ cross2 = aVec13.Crossed( aVec14 ); @@ -556,6 +1030,27 @@ 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) +{ + 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 @@ -617,10 +1112,12 @@ int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const //======================================================================= //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; @@ -667,25 +1164,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: + SMDSAbs_EntityType type = myVolume->GetEntityType(); + if ( myVolume->IsQuadratic() ) + { + int firstMediumInd = myVolume->NbCornerNodes(); + if ( minInd >= firstMediumInd ) + return false; // medium nodes are not linked + if ( maxInd < firstMediumInd ) // both nodes are corners + { + if ( theIgnoreMediumNodes ) + type = SMDSAbs_EntityType( int(type)-1 ); // check linkage of corner nodes + else + return false; // corner nodes are not linked directly in a quadratic cell + } + } + + switch ( type ) { + case SMDSEntity_Tetra: return true; - case 5: + case SMDSEntity_Hexa: + switch ( maxInd - minInd ) { + case 1: return minInd != 3; + case 3: return minInd == 0 || minInd == 4; + case 4: return true; + default:; + } + break; + case SMDSEntity_Pyramid: if ( maxInd == 4 ) return true; switch ( maxInd - minInd ) { @@ -694,7 +1216,7 @@ bool SMDS_VolumeTool::IsLinked (const int theNode1Index, default:; } break; - case 6: + case SMDSEntity_Penta: switch ( maxInd - minInd ) { case 1: return minInd != 2; case 2: return minInd == 0 || minInd == 3; @@ -702,14 +1224,57 @@ 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 SMDSEntity_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 SMDSEntity_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 SMDSEntity_Quad_Pyramid: + { + 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 SMDSEntity_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; } - break; default:; } return false; @@ -731,6 +1296,65 @@ int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const return -1; } +//================================================================================ +/*! + * \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) +{ + faces.clear(); + faces.reserve( NbFaces() ); + 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; ++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(); +} + //======================================================================= //function : IsFreeFace //purpose : check that only one volume is build on the face nodes @@ -751,8 +1375,7 @@ bool SMDS_VolumeTool::IsFreeFace( int faceIndex ) 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(); while ( eIt->more() ) { @@ -760,10 +1383,12 @@ bool SMDS_VolumeTool::IsFreeFace( int faceIndex ) if ( elem != myVolume && elem->GetType() == SMDSAbs_Volume ) { int nbShared = 1; vNbIt = volNbShared.find( elem ); - if ( vNbIt == volNbShared.end() ) + if ( vNbIt == volNbShared.end() ) { volNbShared.insert ( TElemIntMap::value_type( elem, nbShared )); - else + } + else { nbShared = ++(*vNbIt).second; + } if ( nbShared > maxNbShared ) maxNbShared = nbShared; } @@ -779,8 +1404,7 @@ 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++ ) - { + for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) { int nbShared = (*vNbIt).second; if ( nbShared >= 3 ) { SMDS_VolumeTool volume( (*vNbIt).first ); @@ -790,22 +1414,22 @@ bool SMDS_VolumeTool::IsFreeFace( int faceIndex ) continue; // opposite side } // 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 - if ( volNbShared.empty() ) + 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 ( 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++ ) - { + for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) { bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] ); if ( linkShared ) nbSharedLinks++; @@ -900,16 +1524,8 @@ bool SMDS_VolumeTool::setFace( int faceIndex ) // check orientation bool isGoodOri = true; - if (myExternalFaces) { - // get natural orientation - XYZ aNormal, baryCenter, p0 (myPolyedre->GetFaceNode(faceIndex + 1, 1)); - SMDS_VolumeTool vTool (myPolyedre); - vTool.GetFaceNormal(faceIndex, aNormal.x, aNormal.y, aNormal.z); - vTool.GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z); - XYZ insideVec (baryCenter - p0); - if (insideVec.Dot(aNormal) > 0) - isGoodOri = false; - } + if (myExternalFaces) + isGoodOri = IsFaceExternal( faceIndex ); // set face nodes int iNode; @@ -924,7 +1540,8 @@ bool SMDS_VolumeTool::setFace( int faceIndex ) } myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first - } else { + } + else { // choose face node indices switch ( myVolumeNbNodes ) { case 4: @@ -955,6 +1572,34 @@ bool SMDS_VolumeTool::setFace( int faceIndex ) else myFaceNodeIndices = Hexa_F[ faceIndex ]; break; + case 10: + myFaceNbNodes = QuadTetra_nbN[ faceIndex ]; + if ( myExternalFaces ) + myFaceNodeIndices = myVolForward ? QuadTetra_F[ faceIndex ] : QuadTetra_RE[ faceIndex ]; + else + myFaceNodeIndices = myVolForward ? QuadTetra_F[ faceIndex ] : QuadTetra_R[ faceIndex ]; + break; + case 13: + myFaceNbNodes = QuadPyram_nbN[ faceIndex ]; + if ( myExternalFaces ) + myFaceNodeIndices = myVolForward ? QuadPyram_F[ faceIndex ] : QuadPyram_RE[ faceIndex ]; + else + myFaceNodeIndices = myVolForward ? QuadPyram_F[ faceIndex ] : QuadPyram_R[ faceIndex ]; + break; + case 15: + myFaceNbNodes = QuadPenta_nbN[ faceIndex ]; + if ( myExternalFaces ) + myFaceNodeIndices = myVolForward ? QuadPenta_FE[ faceIndex ] : QuadPenta_RE[ faceIndex ]; + else + myFaceNodeIndices = myVolForward ? QuadPenta_F[ faceIndex ] : QuadPenta_R[ faceIndex ]; + break; + case 20: + myFaceNbNodes = QuadHexa_nbN[ faceIndex ]; + if ( myExternalFaces ) + myFaceNodeIndices = myVolForward ? QuadHexa_FE[ faceIndex ] : QuadHexa_RE[ faceIndex ]; + else + myFaceNodeIndices = QuadHexa_F[ faceIndex ]; + break; default: return false; } @@ -983,6 +1628,10 @@ 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: return QUAD_HEXA; default:return UNKNOWN; } } @@ -995,14 +1644,43 @@ 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; + 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; 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; + default: return 0; + } + return 0; +} + // + //======================================================================= //function : GetFaceNodesIndices //purpose : Return the array of face nodes indices @@ -1020,6 +1698,10 @@ const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type, 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 QUAD_TETRA: return QuadTetra_F[ faceIndex ]; + case QUAD_PYRAM: return QuadPyram_F[ faceIndex ]; + case QUAD_PENTA: return external ? QuadPenta_FE[ faceIndex ] : QuadPenta_F[ faceIndex ]; + case QUAD_HEXA: return external ? QuadHexa_FE[ faceIndex ] : QuadHexa_F[ faceIndex ]; default:; } return 0; @@ -1038,8 +1720,31 @@ 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 ]; + case QUAD_HEXA: return QuadHexa_nbN[ faceIndex ]; default:; } return 0; } +//======================================================================= +//function : Get +//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; +}