X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MeshEditor.cxx;h=3ab795264d54f260c7c2bc43e15165ed8aeeede2;hp=71a4d3c1136e7bb1d70a09471d8f6d99c137ee71;hb=05ae10badd80e02346c06faf21f68062400f3f20;hpb=e4737e85f0da6d3f90fd08f6be1c2825195fe16f diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx index 71a4d3c11..3ab795264 100644 --- a/src/SMESH/SMESH_MeshEditor.cxx +++ b/src/SMESH/SMESH_MeshEditor.cxx @@ -1,23 +1,23 @@ // SMESH SMESH : idl implementation based on 'SMESH' unit's classes // // Copyright (C) 2003 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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org +// 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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org // // // @@ -30,33 +30,56 @@ #include "SMDS_FaceOfNodes.hxx" #include "SMDS_VolumeTool.hxx" +#include "SMDS_EdgePosition.hxx" +#include "SMDS_PolyhedralVolumeOfNodes.hxx" +#include "SMDS_FacePosition.hxx" +#include "SMDS_SpacePosition.hxx" + #include "SMESHDS_Group.hxx" #include "SMESHDS_Mesh.hxx" + #include "SMESH_subMesh.hxx" +#include "SMESH_ControlsDef.hxx" #include "utilities.h" -#include #include #include +#include +#include #include #include #include #include +#include +#include #include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include using namespace std; +using namespace SMESH::Controls; -typedef map TNodeNodeMap; -typedef map > TNodeOfNodeListMap; -typedef map > TElemOfNodeListMap; +typedef map TNodeNodeMap; +typedef map > TElemOfNodeListMap; +typedef map > TElemOfElemListMap; +typedef map > TNodeOfNodeListMap; +typedef TNodeOfNodeListMap::iterator TNodeOfNodeListMapItr; +typedef map > TElemOfVecOfNnlmiMap; //======================================================================= //function : SMESH_MeshEditor -//purpose : +//purpose : //======================================================================= SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh ): @@ -76,7 +99,7 @@ bool SMESH_MeshEditor::Remove (const list< int >& theIDs, SMESHDS_Mesh* aMesh = GetMeshDS(); set< SMESH_subMesh *> smmap; - + list::const_iterator it = theIDs.begin(); for ( ; it != theIDs.end(); it++ ) { @@ -202,8 +225,8 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1, if (!F2) return false; // 1 +--+ A theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A - // | /| theTria2: ( B A 2 ) B->1 ( 1 A 2 ) |\ | - // |/ | | \| + // | /| theTria2: ( B A 2 ) B->1 ( 1 A 2 ) |\ | + // |/ | | \| // B +--+ 2 B +--+ 2 // put nodes in array and find out indices of the same ones @@ -323,8 +346,8 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1, if (!F2) return false; // 1 +--+ A tr1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A - // | /| tr2: ( B A 2 ) B->1 ( 1 A 2 ) |\ | - // |/ | | \| + // | /| tr2: ( B A 2 ) B->1 ( 1 A 2 ) |\ | + // |/ | | \| // B +--+ 2 B +--+ 2 // put nodes in array @@ -366,7 +389,7 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1, //MESSAGE( tr1 << tr2 ); return true; - + } //======================================================================= @@ -455,50 +478,84 @@ bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1, //======================================================================= //function : Reorient -//purpose : Reverse the normal of theFace -// Return false if theFace is null +//purpose : Reverse theElement orientation //======================================================================= -bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theFace) +bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem) { - if (!theFace) return false; - const SMDS_FaceOfNodes* F = dynamic_cast( theFace ); - if (!F) return false; + if (!theElem) + return false; + SMDS_ElemIteratorPtr it = theElem->nodesIterator(); + if ( !it || !it->more() ) + return false; - const SMDS_MeshNode* aNodes [4], *tmpNode; - int i = 0; - SMDS_ElemIteratorPtr it = theFace->nodesIterator(); - while ( it->more() ) - aNodes[ i++ ] = static_cast( it->next() ); + switch ( theElem->GetType() ) { - // exchange nodes with indeces 0 and 2 - tmpNode = aNodes[ 0 ]; - aNodes[ 0 ] = aNodes[ 2 ]; - aNodes[ 2 ] = tmpNode; + case SMDSAbs_Edge: + case SMDSAbs_Face: + { + int i = theElem->NbNodes(); + vector aNodes( i ); + while ( it->more() ) + aNodes[ --i ]= static_cast( it->next() ); + return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], theElem->NbNodes() ); + } + case SMDSAbs_Volume: + { + if (theElem->IsPoly()) { + const SMDS_PolyhedralVolumeOfNodes* aPolyedre = + static_cast( theElem ); + if (!aPolyedre) { + MESSAGE("Warning: bad volumic element"); + return false; + } - //MESSAGE( theFace ); + int nbFaces = aPolyedre->NbFaces(); + vector poly_nodes; + vector quantities (nbFaces); - GetMeshDS()->ChangeElementNodes( theFace, aNodes, theFace->NbNodes() ); + // reverse each face of the polyedre + for (int iface = 1; iface <= nbFaces; iface++) { + int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface); + quantities[iface - 1] = nbFaceNodes; - //MESSAGE( theFace ); + for (inode = nbFaceNodes; inode >= 1; inode--) { + const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode); + poly_nodes.push_back(curNode); + } + } - return true; + return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities ); + + } else { + SMDS_VolumeTool vTool; + if ( !vTool.Set( theElem )) + return false; + vTool.Inverse(); + return GetMeshDS()->ChangeElementNodes( theElem, vTool.GetNodes(), vTool.NbNodes() ); + } + } + default:; + } + + return false; } //======================================================================= //function : getBadRate -//purpose : +//purpose : //======================================================================= static double getBadRate (const SMDS_MeshElement* theElem, SMESH::Controls::NumericalFunctorPtr& theCrit) { - TColgp_SequenceOfXYZ P; + SMESH::Controls::TSequenceOfXYZ P; if ( !theElem || !theCrit->GetPoints( theElem, P )) return 1e100; return theCrit->GetBadRate( theCrit->GetValue( P ), theElem->NbNodes() ); + //return theCrit->GetBadRate( theCrit->GetValue( theElem->GetID() ), theElem->NbNodes() ); } - + //======================================================================= //function : QuadToTri //purpose : Cut quadrangles into triangles. @@ -534,7 +591,7 @@ bool SMESH_MeshEditor::QuadToTri (set & theElems, SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] ); SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] ); aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit ); - + SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] ); SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] ); aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit ); @@ -542,7 +599,7 @@ bool SMESH_MeshEditor::QuadToTri (set & theElems, int aShapeId = FindShape( elem ); //MESSAGE( "aBadRate1 = " << aBadRate1 << "; aBadRate2 = " << aBadRate2 // << " ShapeID = " << aShapeId << endl << elem ); - + if ( aBadRate1 <= aBadRate2 ) { // tr1 + tr2 is better aMesh->ChangeElementNodes( elem, aNodes, 3 ); @@ -568,13 +625,49 @@ bool SMESH_MeshEditor::QuadToTri (set & theElems, } //======================================================================= -//function : addToSameGroups +//function : BestSplit +//purpose : Find better diagonal for cutting. +//======================================================================= +int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement* theQuad, + SMESH::Controls::NumericalFunctorPtr theCrit) +{ + if (!theCrit.get()) + return -1; + + if (!theQuad || theQuad->GetType() != SMDSAbs_Face || theQuad->NbNodes() != 4) + return -1; + + // retrieve element nodes + const SMDS_MeshNode* aNodes [4]; + SMDS_ElemIteratorPtr itN = theQuad->nodesIterator(); + int i = 0; + while (itN->more()) + aNodes[ i++ ] = static_cast( itN->next() ); + + // compare two sets of possible triangles + double aBadRate1, aBadRate2; // to what extent a set is bad + SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] ); + SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] ); + aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit ); + + SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] ); + SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] ); + aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit ); + + if (aBadRate1 <= aBadRate2) // tr1 + tr2 is better + return 1; // diagonal 1-3 + + return 2; // diagonal 2-4 +} + +//======================================================================= +//function : AddToSameGroups //purpose : add elemToAdd to the groups the elemInGroups belongs to //======================================================================= -static void addToSameGroups (const SMDS_MeshElement* elemToAdd, - const SMDS_MeshElement* elemInGroups, - SMESHDS_Mesh * aMesh) +void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd, + const SMDS_MeshElement* elemInGroups, + SMESHDS_Mesh * aMesh) { const set& groups = aMesh->GetGroups(); set::const_iterator grIt = groups.begin(); @@ -630,7 +723,7 @@ bool SMESH_MeshEditor::QuadToTri (std::set & theElems, if ( aShapeId ) aMesh->SetMeshElementOnShape( newElem, aShapeId ); - addToSameGroups( newElem, elem, aMesh ); + AddToSameGroups( newElem, elem, aMesh ); } return true; @@ -638,7 +731,7 @@ bool SMESH_MeshEditor::QuadToTri (std::set & theElems, //======================================================================= //function : getAngle -//purpose : +//purpose : //======================================================================= double getAngle(const SMDS_MeshElement * tr1, @@ -649,7 +742,7 @@ double getAngle(const SMDS_MeshElement * tr1, double angle = 2*PI; // bad angle // get normals - TColgp_SequenceOfXYZ P1, P2; + SMESH::Controls::TSequenceOfXYZ P1, P2; if ( !SMESH::Controls::NumericalFunctor::GetPoints( tr1, P1 ) || !SMESH::Controls::NumericalFunctor::GetPoints( tr2, P2 )) return angle; @@ -659,7 +752,7 @@ double getAngle(const SMDS_MeshElement * tr1, gp_Vec N2 = gp_Vec( P2(2) - P2(1) ) ^ gp_Vec( P2(3) - P2(1) ); if ( N2.SquareMagnitude() <= gp::Resolution() ) return angle; - + // find the first diagonal node n1 in the triangles: // take in account a diagonal link orientation const SMDS_MeshElement *nFirst[2], *tr[] = { tr1, tr2 }; @@ -779,7 +872,7 @@ bool SMESH_MeshEditor::TriToQuad (set & theElems, itLE = mapLi_listEl.find( linkID ); if ( itLE != mapLi_listEl.end() ) { - if ((*itLE).second.size() > 1 ) // consider only 2 elems adjacent by a link + if ((*itLE).second.size() > 1 ) // consider only 2 elems adjacent by a link continue; const SMDS_MeshElement* elem2 = (*itLE).second.front(); // if ( FindShape( elem ) != FindShape( elem2 )) @@ -809,7 +902,7 @@ bool SMESH_MeshEditor::TriToQuad (set & theElems, } // Algo: fuse triangles into quadrangles - + while ( ! mapEl_setLi.empty() ) { // Look for the start element: @@ -951,17 +1044,17 @@ bool SMESH_MeshEditor::TriToQuad (set & theElems, } // if ( startElem ) } // while ( startElem || !startLinks.empty() ) } // while ( ! mapEl_setLi.empty() ) - + return true; } -#define DUMPSO(txt) \ +/*#define DUMPSO(txt) \ // cout << txt << endl; //============================================================================= -/*! - * - */ +// +// +// //============================================================================= static void swap( int i1, int i2, int idNodes[], gp_Pnt P[] ) { @@ -1041,7 +1134,7 @@ bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh, } DUMPSO( "========================================"); - + set faceNodes; // ids of bottom face nodes, to be found set checkedId1; // ids of tried 2-nd nodes Standard_Real leastDist = DBL_MAX; // dist of the 4-th node from 123 plane @@ -1060,7 +1153,7 @@ bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh, checkedId1.insert ( id1 ); break; } - + // Find the 3-d node so that 1-2-3 triangle to be on a hexa face, // ie that all but meybe one (id3 which is on the same face) nodes // lay on the same side from the triangle plane. @@ -1161,7 +1254,7 @@ bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh, } } - + // Set nodes of the found bottom face in good order DUMPSO( " Found bottom face: "); i = SortQuadNodes( theMesh, idNodes ); @@ -1185,7 +1278,7 @@ bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh, Standard_Real upDirSize = upDir.Magnitude(); if ( upDirSize <= gp::Resolution() ) return false; upDir / upDirSize; - + // Assure that the bottom face normal points up gp_Vec Nb = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) ); Nb += gp_Vec (P[0], P[2]).Crossed( gp_Vec (P[0], P[3]) ); @@ -1233,7 +1326,7 @@ bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh, // } return true; -} +}*/ //======================================================================= //function : laplacianSmooth @@ -1241,52 +1334,76 @@ bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh, // connected to that node along an element edge //======================================================================= -void laplacianSmooth(SMESHDS_Mesh * theMesh, - const SMDS_MeshNode* theNode, - const set & theElems, - const set & theFixedNodes) +void laplacianSmooth(const SMDS_MeshNode* theNode, + const Handle(Geom_Surface)& theSurface, + map< const SMDS_MeshNode*, gp_XY* >& theUVMap) { // find surrounding nodes + set< const SMDS_MeshNode* > nodeSet; SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(); while ( elemIt->more() ) { const SMDS_MeshElement* elem = elemIt->next(); - if ( theElems.find( elem ) == theElems.end() ) + if ( elem->GetType() != SMDSAbs_Face ) continue; - int i = 0, iNode = 0; - const SMDS_MeshNode* aNodes [4]; + // put all nodes in array + int nbNodes = 0, iNode = 0; + vector< const SMDS_MeshNode*> aNodes( elem->NbNodes() ); SMDS_ElemIteratorPtr itN = elem->nodesIterator(); while ( itN->more() ) { - aNodes[ i ] = static_cast( itN->next() ); - if ( aNodes[ i ] == theNode ) - iNode = i; - else - nodeSet.insert( aNodes[ i ] ); - i++; - } - if ( elem->NbNodes() == 4 ) { // remove an opposite node - iNode += ( iNode < 2 ) ? 2 : -2; - nodeSet.erase( aNodes[ iNode ]); + aNodes[ nbNodes ] = static_cast( itN->next() ); + if ( aNodes[ nbNodes ] == theNode ) + iNode = nbNodes; // index of theNode within aNodes + nbNodes++; } + // add linked nodes + int iAfter = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1; + nodeSet.insert( aNodes[ iAfter ]); + int iBefore = ( iNode == 0 ) ? nbNodes - 1 : iNode - 1; + nodeSet.insert( aNodes[ iBefore ]); } // compute new coodrs + double coord[] = { 0., 0., 0. }; set< const SMDS_MeshNode* >::iterator nodeSetIt = nodeSet.begin(); for ( ; nodeSetIt != nodeSet.end(); nodeSetIt++ ) { const SMDS_MeshNode* node = (*nodeSetIt); - coord[0] += node->X(); - coord[1] += node->Y(); - coord[2] += node->Z(); - } - double nbNodes = nodeSet.size(); - theMesh->MoveNode (theNode, - coord[0]/nbNodes, - coord[1]/nbNodes, - coord[2]/nbNodes); + if ( theSurface.IsNull() ) { // smooth in 3D + coord[0] += node->X(); + coord[1] += node->Y(); + coord[2] += node->Z(); + } + else { // smooth in 2D + ASSERT( theUVMap.find( node ) != theUVMap.end() ); + gp_XY* uv = theUVMap[ node ]; + coord[0] += uv->X(); + coord[1] += uv->Y(); + } + } + int nbNodes = nodeSet.size(); + if ( !nbNodes ) + return; + coord[0] /= nbNodes; + coord[1] /= nbNodes; + + if ( !theSurface.IsNull() ) { + ASSERT( theUVMap.find( theNode ) != theUVMap.end() ); + theUVMap[ theNode ]->SetCoord( coord[0], coord[1] ); + gp_Pnt p3d = theSurface->Value( coord[0], coord[1] ); + coord[0] = p3d.X(); + coord[1] = p3d.Y(); + coord[2] = p3d.Z(); + } + else + coord[2] /= nbNodes; + + // move node + + const_cast< SMDS_MeshNode* >( theNode )->setXYZ(coord[0],coord[1],coord[2]); } //======================================================================= @@ -1295,33 +1412,38 @@ void laplacianSmooth(SMESHDS_Mesh * theMesh, // surrounding elements //======================================================================= -void centroidalSmooth(SMESHDS_Mesh * theMesh, - const SMDS_MeshNode* theNode, - const set & theElems, - const set & theFixedNodes) +void centroidalSmooth(const SMDS_MeshNode* theNode, + const Handle(Geom_Surface)& theSurface, + map< const SMDS_MeshNode*, gp_XY* >& theUVMap) { gp_XYZ aNewXYZ(0.,0.,0.); SMESH::Controls::Area anAreaFunc; double totalArea = 0.; int nbElems = 0; + // compute new XYZ + SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(); while ( elemIt->more() ) { const SMDS_MeshElement* elem = elemIt->next(); - if ( theElems.find( elem ) == theElems.end() ) + if ( elem->GetType() != SMDSAbs_Face ) continue; - nbElems++; gp_XYZ elemCenter(0.,0.,0.); - TColgp_SequenceOfXYZ aNodePoints; + SMESH::Controls::TSequenceOfXYZ aNodePoints; SMDS_ElemIteratorPtr itN = elem->nodesIterator(); while ( itN->more() ) { const SMDS_MeshNode* aNode = static_cast( itN->next() ); gp_XYZ aP( aNode->X(), aNode->Y(), aNode->Z() ); - aNodePoints.Append( aP ); + aNodePoints.push_back( aP ); + if ( !theSurface.IsNull() ) { // smooth in 2D + ASSERT( theUVMap.find( aNode ) != theUVMap.end() ); + gp_XY* uv = theUVMap[ aNode ]; + aP.SetCoord( uv->X(), uv->Y(), 0. ); + } elemCenter += aP; } double elemArea = anAreaFunc.GetValue( aNodePoints ); @@ -1330,10 +1452,38 @@ void centroidalSmooth(SMESHDS_Mesh * theMesh, aNewXYZ += elemCenter * elemArea; } aNewXYZ /= totalArea; - theMesh->MoveNode (theNode, - aNewXYZ.X(), - aNewXYZ.Y(), - aNewXYZ.Z()); + if ( !theSurface.IsNull() ) { + ASSERT( theUVMap.find( theNode ) != theUVMap.end() ); + theUVMap[ theNode ]->SetCoord( aNewXYZ.X(), aNewXYZ.Y() ); + aNewXYZ = theSurface->Value( aNewXYZ.X(), aNewXYZ.Y() ).XYZ(); + } + + // move node + + const_cast< SMDS_MeshNode* >( theNode )->setXYZ(aNewXYZ.X(),aNewXYZ.Y(),aNewXYZ.Z()); +} + +//======================================================================= +//function : getClosestUV +//purpose : return UV of closest projection +//======================================================================= + +static bool getClosestUV (Extrema_GenExtPS& projector, + const gp_Pnt& point, + gp_XY & result) +{ + projector.Perform( point ); + if ( projector.IsDone() ) { + double u, v, minVal = DBL_MAX; + for ( int i = projector.NbExt(); i > 0; i-- ) + if ( projector.Value( i ) < minVal ) { + minVal = projector.Value( i ); + projector.Point( i ).Parameter( u, v ); + } + result.SetCoord( u, v ); + return true; + } + return false; } //======================================================================= @@ -1350,129 +1500,472 @@ void SMESH_MeshEditor::Smooth (set & theElems, set & theFixedNodes, const SmoothMethod theSmoothMethod, const int theNbIterations, - double theTgtAspectRatio) + double theTgtAspectRatio, + const bool the2D) { MESSAGE((theSmoothMethod==LAPLACIAN ? "LAPLACIAN" : "CENTROIDAL") << "--::Smooth()"); + if ( theTgtAspectRatio < 1.0 ) + theTgtAspectRatio = 1.0; + + SMESH::Controls::AspectRatio aQualityFunc; + SMESHDS_Mesh* aMesh = GetMeshDS(); + if ( theElems.empty() ) { - // add all faces + // add all faces to theElems SMDS_FaceIteratorPtr fIt = aMesh->facesIterator(); while ( fIt->more() ) theElems.insert( fIt->next() ); } + // get all face ids theElems are on + set< int > faceIdSet; + set< const SMDS_MeshElement* >::iterator itElem; + if ( the2D ) + for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) { + int fId = FindShape( *itElem ); + // check that corresponding submesh exists and a shape is face + if (fId && + faceIdSet.find( fId ) == faceIdSet.end() && + aMesh->MeshElements( fId )) { + TopoDS_Shape F = aMesh->IndexToShape( fId ); + if ( !F.IsNull() && F.ShapeType() == TopAbs_FACE ) + faceIdSet.insert( fId ); + } + } + faceIdSet.insert( 0 ); // to smooth elements that are not on any TopoDS_Face - set setMovableNodes; - - // Fill setMovableNodes + // =============================================== + // smooth elements on each TopoDS_Face separately + // =============================================== - map< const SMDS_MeshNode*, int > mapNodeNbFaces; - set< const SMDS_MeshElement* >::iterator itElem; - for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) + set< int >::reverse_iterator fId = faceIdSet.rbegin(); // treate 0 fId at the end + for ( ; fId != faceIdSet.rend(); ++fId ) { - const SMDS_MeshElement* elem = (*itElem); - if ( !elem || elem->GetType() != SMDSAbs_Face ) - continue; - - SMDS_ElemIteratorPtr itN = elem->nodesIterator(); - while ( itN->more() ) { - const SMDS_MeshNode* node = - static_cast( itN->next() ); + // get face surface and submesh + Handle(Geom_Surface) surface; + SMESHDS_SubMesh* faceSubMesh = 0; + TopoDS_Face face; + double fToler2 = 0, vPeriod = 0., uPeriod = 0., f,l; + double u1 = 0, u2 = 0, v1 = 0, v2 = 0; + bool isUPeriodic = false, isVPeriodic = false; + if ( *fId ) { + face = TopoDS::Face( aMesh->IndexToShape( *fId )); + surface = BRep_Tool::Surface( face ); + faceSubMesh = aMesh->MeshElements( *fId ); + fToler2 = BRep_Tool::Tolerance( face ); + fToler2 *= fToler2 * 10.; + isUPeriodic = surface->IsUPeriodic(); + if ( isUPeriodic ) + vPeriod = surface->UPeriod(); + isVPeriodic = surface->IsVPeriodic(); + if ( isVPeriodic ) + uPeriod = surface->VPeriod(); + surface->Bounds( u1, u2, v1, v2 ); + } + // --------------------------------------------------------- + // for elements on a face, find movable and fixed nodes and + // compute UV for them + // --------------------------------------------------------- + bool checkBoundaryNodes = false; + set setMovableNodes; + map< const SMDS_MeshNode*, gp_XY* > uvMap, uvMap2; + list< gp_XY > listUV; // uvs the 2 uvMaps refer to + list< const SMDS_MeshElement* > elemsOnFace; + + Extrema_GenExtPS projector; + GeomAdaptor_Surface surfAdaptor; + if ( !surface.IsNull() ) { + surfAdaptor.Load( surface ); + projector.Initialize( surfAdaptor, 20,20, 1e-5,1e-5 ); + } + int nbElemOnFace = 0; + itElem = theElems.begin(); + // loop on not yet smoothed elements: look for elems on a face + while ( itElem != theElems.end() ) + { + if ( faceSubMesh && nbElemOnFace == faceSubMesh->NbElements() ) + break; // all elements found - if ( theFixedNodes.find( node ) != theFixedNodes.end() ) + const SMDS_MeshElement* elem = (*itElem); + if ( !elem || elem->GetType() != SMDSAbs_Face || elem->NbNodes() < 3 || + ( faceSubMesh && !faceSubMesh->Contains( elem ))) { + ++itElem; continue; + } + elemsOnFace.push_back( elem ); + theElems.erase( itElem++ ); + nbElemOnFace++; + + // get movable nodes of elem + const SMDS_MeshNode* node; + SMDS_TypeOfPosition posType; + SMDS_ElemIteratorPtr itN = elem->nodesIterator(); + while ( itN->more() ) { + node = static_cast( itN->next() ); + const SMDS_PositionPtr& pos = node->GetPosition(); + posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE; + if (posType != SMDS_TOP_EDGE && + posType != SMDS_TOP_VERTEX && + theFixedNodes.find( node ) == theFixedNodes.end()) + { + // check if all faces around the node are on faceSubMesh + // because a node on edge may be bound to face + SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(); + bool all = true; + if ( faceSubMesh ) { + while ( eIt->more() && all ) { + const SMDS_MeshElement* e = eIt->next(); + if ( e->GetType() == SMDSAbs_Face ) + all = faceSubMesh->Contains( e ); + } + } + if ( all ) + setMovableNodes.insert( node ); + else + checkBoundaryNodes = true; + } + if ( posType == SMDS_TOP_3DSPACE ) + checkBoundaryNodes = true; + } - // if node is on edge => it is fixed - SMDS_PositionPtr aPositionPtr = node->GetPosition(); - if ( aPositionPtr.get() && - (aPositionPtr->GetTypeOfPosition() == SMDS_TOP_EDGE || - aPositionPtr->GetTypeOfPosition() == SMDS_TOP_VERTEX)) { - theFixedNodes.insert( node ); + if ( surface.IsNull() ) continue; + + // get nodes to check UV + list< const SMDS_MeshNode* > uvCheckNodes; + itN = elem->nodesIterator(); + while ( itN->more() ) { + node = static_cast( itN->next() ); + if ( uvMap.find( node ) == uvMap.end() ) + uvCheckNodes.push_back( node ); + // add nodes of elems sharing node +// SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(); +// while ( eIt->more() ) { +// const SMDS_MeshElement* e = eIt->next(); +// if ( e != elem && e->GetType() == SMDSAbs_Face ) { +// SMDS_ElemIteratorPtr nIt = e->nodesIterator(); +// while ( nIt->more() ) { +// const SMDS_MeshNode* n = +// static_cast( nIt->next() ); +// if ( uvMap.find( n ) == uvMap.end() ) +// uvCheckNodes.push_back( n ); +// } +// } +// } + } + // check UV on face + list< const SMDS_MeshNode* >::iterator n = uvCheckNodes.begin(); + for ( ; n != uvCheckNodes.end(); ++n ) + { + node = *n; + gp_XY uv( 0, 0 ); + const SMDS_PositionPtr& pos = node->GetPosition(); + posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE; + // get existing UV + switch ( posType ) { + case SMDS_TOP_FACE: { + SMDS_FacePosition* fPos = ( SMDS_FacePosition* ) pos.get(); + uv.SetCoord( fPos->GetUParameter(), fPos->GetVParameter() ); + break; + } + case SMDS_TOP_EDGE: { + TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() ); + Handle(Geom2d_Curve) pcurve; + if ( !S.IsNull() && S.ShapeType() == TopAbs_EDGE ) + pcurve = BRep_Tool::CurveOnSurface( TopoDS::Edge( S ), face, f,l ); + if ( !pcurve.IsNull() ) { + double u = (( SMDS_EdgePosition* ) pos.get() )->GetUParameter(); + uv = pcurve->Value( u ).XY(); + } + break; + } + case SMDS_TOP_VERTEX: { + TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() ); + if ( !S.IsNull() && S.ShapeType() == TopAbs_VERTEX ) + uv = BRep_Tool::Parameters( TopoDS::Vertex( S ), face ).XY(); + break; + } + default:; + } + // check existing UV + bool project = true; + gp_Pnt pNode ( node->X(), node->Y(), node->Z() ); + double dist1 = DBL_MAX, dist2 = 0; + if ( posType != SMDS_TOP_3DSPACE ) { + dist1 = pNode.SquareDistance( surface->Value( uv.X(), uv.Y() )); + project = dist1 > fToler2; + } + if ( project ) { // compute new UV + gp_XY newUV; + if ( !getClosestUV( projector, pNode, newUV )) { + MESSAGE("Node Projection Failed " << node); + } + else { + if ( isUPeriodic ) + newUV.SetX( ElCLib::InPeriod( newUV.X(), u1, u2 )); + if ( isVPeriodic ) + newUV.SetY( ElCLib::InPeriod( newUV.Y(), v1, v2 )); + // check new UV + if ( posType != SMDS_TOP_3DSPACE ) + dist2 = pNode.SquareDistance( surface->Value( newUV.X(), newUV.Y() )); + if ( dist2 < dist1 ) + uv = newUV; + } + } + // store UV in the map + listUV.push_back( uv ); + uvMap.insert( make_pair( node, &listUV.back() )); + } + } // loop on not yet smoothed elements + + if ( !faceSubMesh || nbElemOnFace != faceSubMesh->NbElements() ) + checkBoundaryNodes = true; + + // fix nodes on mesh boundary + + if ( checkBoundaryNodes ) + { + typedef pair TLink; + map< TLink, int > linkNbMap; // how many times a link encounters in elemsOnFace + map< TLink, int >::iterator link_nb; + // put all elements links to linkNbMap + list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin(); + for ( ; elemIt != elemsOnFace.end(); ++elemIt ) + { + // put elem nodes in array + vector< const SMDS_MeshNode* > nodes; + nodes.reserve( (*elemIt)->NbNodes() + 1 ); + SMDS_ElemIteratorPtr itN = (*elemIt)->nodesIterator(); + while ( itN->more() ) + nodes.push_back( static_cast( itN->next() )); + nodes.push_back( nodes.front() ); + // loop on elem links: insert them in linkNbMap + for ( int iN = 1; iN < nodes.size(); ++iN ) { + TLink link; + if ( nodes[ iN-1 ]->GetID() < nodes[ iN ]->GetID() ) + link = make_pair( nodes[ iN-1 ], nodes[ iN ] ); + else + link = make_pair( nodes[ iN ], nodes[ iN-1 ] ); + link_nb = linkNbMap.find( link ); + if ( link_nb == linkNbMap.end() ) + linkNbMap.insert( make_pair ( link, 1 )); + else + link_nb->second++; + } + } + // remove nodes that are in links encountered only once from setMovableNodes + for ( link_nb = linkNbMap.begin(); link_nb != linkNbMap.end(); ++link_nb ) { + if ( link_nb->second == 1 ) { + setMovableNodes.erase( link_nb->first.first ); + setMovableNodes.erase( link_nb->first.second ); + } } - // fill mapNodeNbFaces in order to detect fixed boundary nodes - map::iterator nodeNbFacesIt = - mapNodeNbFaces.find ( node ); - if ( nodeNbFacesIt == mapNodeNbFaces.end() ) - mapNodeNbFaces.insert( map::value_type( node, 1 )); - else - (*nodeNbFacesIt).second++; } - } - // put not fixed nodes in setMovableNodes - map::iterator nodeNbFacesIt = - mapNodeNbFaces.begin(); - for ( ; nodeNbFacesIt != mapNodeNbFaces.end(); nodeNbFacesIt++ ) { - const SMDS_MeshNode* node = (*nodeNbFacesIt).first; - // a node is on free boundary if it is shared by 1-2 faces - if ( (*nodeNbFacesIt).second > 2 ) - setMovableNodes.insert( node ); - else - theFixedNodes.insert( node ); - } - // SMOOTHING // + // ----------------------------------------------------- + // for nodes on seam edge, compute one more UV ( uvMap2 ); + // find movable nodes linked to nodes on seam and which + // are to be smoothed using the second UV ( uvMap2 ) + // ----------------------------------------------------- - if ( theTgtAspectRatio < 1.0 ) - theTgtAspectRatio = 1.0; + set nodesNearSeam; // to smooth using uvMap2 + if ( !surface.IsNull() ) + { + TopExp_Explorer eExp( face, TopAbs_EDGE ); + for ( ; eExp.More(); eExp.Next() ) + { + TopoDS_Edge edge = TopoDS::Edge( eExp.Current() ); + if ( !BRep_Tool::IsClosed( edge, face )) + continue; + SMESHDS_SubMesh* sm = aMesh->MeshElements( edge ); + if ( !sm ) continue; + // find out which parameter varies for a node on seam + double f,l; + gp_Pnt2d uv1, uv2; + Handle(Geom2d_Curve) pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l ); + if ( pcurve.IsNull() ) continue; + uv1 = pcurve->Value( f ); + edge.Reverse(); + pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l ); + if ( pcurve.IsNull() ) continue; + uv2 = pcurve->Value( f ); + int iPar = Abs( uv1.X() - uv2.X() ) > Abs( uv1.Y() - uv2.Y() ) ? 1 : 2; + // assure uv1 < uv2 + if ( uv1.Coord( iPar ) > uv2.Coord( iPar )) { + gp_Pnt2d tmp = uv1; uv1 = uv2; uv2 = tmp; + } + // get nodes on seam and its vertices + list< const SMDS_MeshNode* > seamNodes; + SMDS_NodeIteratorPtr nSeamIt = sm->GetNodes(); + while ( nSeamIt->more() ) + seamNodes.push_back( nSeamIt->next() ); + TopExp_Explorer vExp( edge, TopAbs_VERTEX ); + for ( ; vExp.More(); vExp.Next() ) { + sm = aMesh->MeshElements( vExp.Current() ); + if ( sm ) { + nSeamIt = sm->GetNodes(); + while ( nSeamIt->more() ) + seamNodes.push_back( nSeamIt->next() ); + } + } + // loop on nodes on seam + list< const SMDS_MeshNode* >::iterator noSeIt = seamNodes.begin(); + for ( ; noSeIt != seamNodes.end(); ++noSeIt ) + { + const SMDS_MeshNode* nSeam = *noSeIt; + map< const SMDS_MeshNode*, gp_XY* >::iterator n_uv = uvMap.find( nSeam ); + if ( n_uv == uvMap.end() ) + continue; + // set the first UV + n_uv->second->SetCoord( iPar, uv1.Coord( iPar )); + // set the second UV + listUV.push_back( *n_uv->second ); + listUV.back().SetCoord( iPar, uv2.Coord( iPar )); + if ( uvMap2.empty() ) + uvMap2 = uvMap; // copy the uvMap contents + uvMap2[ nSeam ] = &listUV.back(); + + // collect movable nodes linked to ones on seam in nodesNearSeam + SMDS_ElemIteratorPtr eIt = nSeam->GetInverseElementIterator(); + while ( eIt->more() ) + { + const SMDS_MeshElement* e = eIt->next(); + if ( e->GetType() != SMDSAbs_Face ) + continue; + int nbUseMap1 = 0, nbUseMap2 = 0; + SMDS_ElemIteratorPtr nIt = e->nodesIterator(); + while ( nIt->more() ) + { + const SMDS_MeshNode* n = + static_cast( nIt->next() ); + if (n == nSeam || + setMovableNodes.find( n ) == setMovableNodes.end() ) + continue; + // add only nodes being closer to uv2 than to uv1 + gp_Pnt pMid (0.5 * ( n->X() + nSeam->X() ), + 0.5 * ( n->Y() + nSeam->Y() ), + 0.5 * ( n->Z() + nSeam->Z() )); + gp_XY uv; + getClosestUV( projector, pMid, uv ); + if ( uv.Coord( iPar ) > uvMap[ n ]->Coord( iPar ) ) { + nodesNearSeam.insert( n ); + nbUseMap2++; + } + else + nbUseMap1++; + } + // for centroidalSmooth all element nodes must + // be on one side of a seam + if ( theSmoothMethod == CENTROIDAL && nbUseMap1 && nbUseMap2 ) + { + SMDS_ElemIteratorPtr nIt = e->nodesIterator(); + while ( nIt->more() ) { + const SMDS_MeshNode* n = + static_cast( nIt->next() ); + setMovableNodes.erase( n ); + } + } + } + } // loop on nodes on seam + } // loop on edge of a face + } // if ( !face.IsNull() ) - SMESH::Controls::AspectRatio aQualityFunc; + if ( setMovableNodes.empty() ) { + MESSAGE( "Face id : " << *fId << " - NO SMOOTHING: no nodes to move!!!"); + continue; // goto next face + } - for ( int it = 0; it < theNbIterations; it++ ) - { - Standard_Real maxDisplacement = 0.; - set::iterator movableNodesIt - = setMovableNodes.begin(); - for ( ; movableNodesIt != setMovableNodes.end(); movableNodesIt++ ) + // ------------- + // SMOOTHING // + // ------------- + + int it = -1; + double maxRatio = -1., maxDisplacement = -1.; + set::iterator nodeToMove; + for ( it = 0; it < theNbIterations; it++ ) { - const SMDS_MeshNode* node = (*movableNodesIt); - gp_XYZ aPrevPos ( node->X(), node->Y(), node->Z() ); + maxDisplacement = 0.; + nodeToMove = setMovableNodes.begin(); + for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) + { + const SMDS_MeshNode* node = (*nodeToMove); + gp_XYZ aPrevPos ( node->X(), node->Y(), node->Z() ); - // smooth - if ( theSmoothMethod == LAPLACIAN ) - laplacianSmooth( aMesh, node, theElems, theFixedNodes ); - else - centroidalSmooth( aMesh, node, theElems, theFixedNodes ); + // smooth + bool map2 = ( nodesNearSeam.find( node ) != nodesNearSeam.end() ); + if ( theSmoothMethod == LAPLACIAN ) + laplacianSmooth( node, surface, map2 ? uvMap2 : uvMap ); + else + centroidalSmooth( node, surface, map2 ? uvMap2 : uvMap ); - // displacement - gp_XYZ aNewPos ( node->X(), node->Y(), node->Z() ); - Standard_Real aDispl = (aPrevPos - aNewPos).SquareModulus(); - if ( aDispl > maxDisplacement ) - maxDisplacement = aDispl; - } - // no node movement => exit - if ( maxDisplacement < 1.e-16 ) { - MESSAGE("-- no node movement -- maxDisplacement: " << maxDisplacement << " it "<< it); - break; - } + // node displacement + gp_XYZ aNewPos ( node->X(), node->Y(), node->Z() ); + Standard_Real aDispl = (aPrevPos - aNewPos).SquareModulus(); + if ( aDispl > maxDisplacement ) + maxDisplacement = aDispl; + } + // no node movement => exit + if ( maxDisplacement < 1.e-16 ) { + MESSAGE("-- no node movement --"); + break; + } - // check elements quality - double maxRatio = 0; - for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) + // check elements quality + maxRatio = 0; + list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin(); + for ( ; elemIt != elemsOnFace.end(); ++elemIt ) + { + const SMDS_MeshElement* elem = (*elemIt); + if ( !elem || elem->GetType() != SMDSAbs_Face ) + continue; + SMESH::Controls::TSequenceOfXYZ aPoints; + if ( aQualityFunc.GetPoints( elem, aPoints )) { + double aValue = aQualityFunc.GetValue( aPoints ); + if ( aValue > maxRatio ) + maxRatio = aValue; + } + } + if ( maxRatio <= theTgtAspectRatio ) { + MESSAGE("-- quality achived --"); + break; + } + if (it+1 == theNbIterations) { + MESSAGE("-- Iteration limit exceeded --"); + } + } // smoothing iterations + + MESSAGE(" Face id: " << *fId << + " Nb iterstions: " << it << + " Displacement: " << maxDisplacement << + " Aspect Ratio " << maxRatio); + + // --------------------------------------- + // new nodes positions are computed, + // record movement in DS and set new UV + // --------------------------------------- + + nodeToMove = setMovableNodes.begin(); + for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) { - const SMDS_MeshElement* elem = (*itElem); - if ( !elem || elem->GetType() != SMDSAbs_Face ) - continue; - TColgp_SequenceOfXYZ aPoints; - if ( aQualityFunc.GetPoints( elem, aPoints )) { - double aValue = aQualityFunc.GetValue( aPoints ); - if ( aValue > maxRatio ) - maxRatio = aValue; + SMDS_MeshNode* node = const_cast< SMDS_MeshNode* > (*nodeToMove); + aMesh->MoveNode( node, node->X(), node->Y(), node->Z() ); + map< const SMDS_MeshNode*, gp_XY* >::iterator node_uv = uvMap.find( node ); + if ( node_uv != uvMap.end() ) { + gp_XY* uv = node_uv->second; + node->SetPosition + ( SMDS_PositionPtr( new SMDS_FacePosition( *fId, uv->X(), uv->Y() ))); } } - if ( maxRatio <= theTgtAspectRatio ) { - MESSAGE("-- quality achived -- maxRatio " << maxRatio << " it "<< it); - break; - } - if (it+1 == theNbIterations) { - MESSAGE("-- Iteration limit exceeded --"); - } - } + + } // loop on face ids } //======================================================================= //function : isReverse -//purpose : +//purpose : Return true if normal of prevNodes is not co-directied with +// gp_Vec(prevNodes[iNotSame],nextNodes[iNotSame]). +// iNotSame is where prevNodes and nextNodes are different //======================================================================= static bool isReverse(const SMDS_MeshNode* prevNodes[], @@ -1500,42 +1993,39 @@ static bool isReverse(const SMDS_MeshNode* prevNodes[], //======================================================================= //function : sweepElement -//purpose : +//purpose : //======================================================================= -static void sweepElement(SMESHDS_Mesh* aMesh, - const SMDS_MeshElement* elem, - const TNodeOfNodeListMap& mapNewNodes ) +static void sweepElement(SMESHDS_Mesh* aMesh, + const SMDS_MeshElement* elem, + const vector & newNodesItVec, + list& newElems) { // Loop on elem nodes: // find new nodes and detect same nodes indices - list::const_iterator itNN[ 4 ]; - const SMDS_MeshNode* prevNod[ 4 ], *nextNod[ 4 ]; - int nbSame = 0, iNotSameNode = 0, iSameNode = 0; - - TNodeOfNodeListMap::const_iterator mapIt; - int iNode = 0; - SMDS_ElemIteratorPtr itN = elem->nodesIterator(); - while ( itN->more() ) + int nbNodes = elem->NbNodes(); + list::const_iterator itNN[ nbNodes ]; + const SMDS_MeshNode* prevNod[ nbNodes ], *nextNod[ nbNodes ]; + int iNode, nbSame = 0, iNotSameNode = 0, iSameNode = 0; + + for ( iNode = 0; iNode < nbNodes; iNode++ ) { - const SMDS_MeshNode* node = - static_cast( itN->next() ); - mapIt = mapNewNodes.find( node ); - if ( mapIt == mapNewNodes.end() ) - return; // not duplicated node + TNodeOfNodeListMapItr nnIt = newNodesItVec[ iNode ]; + const SMDS_MeshNode* node = nnIt->first; + const list< const SMDS_MeshNode* > & listNewNodes = nnIt->second; + if ( listNewNodes.empty() ) + return; - itNN[ iNode ] = (*mapIt).second.begin(); + itNN[ iNode ] = listNewNodes.begin(); prevNod[ iNode ] = node; - nextNod[ iNode ] = (*mapIt).second.front(); + nextNod[ iNode ] = listNewNodes.front(); if ( prevNod[ iNode ] != nextNod [ iNode ]) iNotSameNode = iNode; else { iSameNode = iNode; nbSame++; } - iNode++; } - int nbNodes = iNode; if ( nbSame == nbNodes || nbSame > 2) { MESSAGE( " Too many same nodes of element " << elem->GetID() ); return; @@ -1550,8 +2040,8 @@ static void sweepElement(SMESHDS_Mesh* aMesh, // check element orientation int i0 = 0, i2 = 2; - if ( nbNodes > 2 && isReverse( prevNod, nextNod, nbNodes, iNotSameNode )) { -// MESSAGE("Reversed elem " << elem->GetID() ); + if ( nbNodes > 2 && !isReverse( prevNod, nextNod, nbNodes, iNotSameNode )) { + //MESSAGE("Reversed elem " << elem ); i0 = 2; i2 = 0; if ( nbSame > 0 ) { @@ -1562,7 +2052,7 @@ static void sweepElement(SMESHDS_Mesh* aMesh, } // make new elements - int iStep, nbSteps = (*mapIt).second.size(); + int iStep, nbSteps = newNodesItVec[ 0 ]->second.size(); for (iStep = 0; iStep < nbSteps; iStep++ ) { // get next nodes @@ -1570,69 +2060,101 @@ static void sweepElement(SMESHDS_Mesh* aMesh, nextNod[ iNode ] = *itNN[ iNode ]; itNN[ iNode ]++; } + SMDS_MeshElement* aNewElem = 0; switch ( nbNodes ) { + case 0: + return; + case 1: { // NODE + if ( nbSame == 0 ) + aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ] ); + break; + } case 2: { // EDGE if ( nbSame == 0 ) - aMesh->AddFace( prevNod[ 0 ], prevNod[ 1 ], nextNod[ 1 ], nextNod[ 0 ] ); + aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ], + nextNod[ 1 ], nextNod[ 0 ] ); else - aMesh->AddFace( prevNod[ 0 ], prevNod[ 1 ], nextNod[ iNotSameNode ] ); + aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ], + nextNod[ iNotSameNode ] ); break; } case 3: { // TRIANGLE - if ( nbSame == 0 ) // --- 1 pentahedron - { - aMesh->AddVolume (prevNod[ i2 ], prevNod[ 1 ], prevNod[ i0 ], - nextNod[ i2 ], nextNod[ 1 ], nextNod[ i0 ] ); - } - else if ( nbSame == 1 ) // --- 2 tetrahedrons - { - aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], - nextNod[ iBeforeSame ]); - aMesh->AddVolume (nextNod[ i2 ], nextNod[ 1 ], nextNod[ i0 ], - prevNod[ iAfterSame ]); - } - else // 2 same nodes: --- 1 tetrahedron - { - aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], - nextNod[ iNotSameNode ]); - } + if ( nbSame == 0 ) // --- pentahedron + aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], + nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ] ); + + else if ( nbSame == 1 ) // --- pyramid + aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ], + nextNod[ iAfterSame ], nextNod[ iBeforeSame ], + nextNod[ iSameNode ]); + + else // 2 same nodes: --- tetrahedron + aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], + nextNod[ iNotSameNode ]); break; } case 4: { // QUADRANGLE - if ( nbSame == 0 ) // --- 1 hexahedron + if ( nbSame == 0 ) // --- hexahedron + aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], prevNod[ 3 ], + nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ], nextNod[ 3 ]); + + else if ( nbSame == 1 ) // --- pyramid + pentahedron { - aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], prevNod[ 3 ], - nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ], nextNod[ 3 ]); + aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ], + nextNod[ iAfterSame ], nextNod[ iBeforeSame ], + nextNod[ iSameNode ]); + newElems.push_back( aNewElem ); + aNewElem = aMesh->AddVolume (prevNod[ iAfterSame ], prevNod[ iOpposSame ], + prevNod[ iBeforeSame ], nextNod[ iAfterSame ], + nextNod[ iOpposSame ], nextNod[ iBeforeSame ] ); } - else if ( nbSame == 1 ) // --- 2 tetrahedrons + 1 pentahedron - { - aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iSameNode ], - prevNod[ iAfterSame ], nextNod[ iBeforeSame ]); - aMesh->AddVolume (nextNod[ iAfterSame ], nextNod[ iSameNode ], - nextNod[ iBeforeSame ], prevNod[ iAfterSame ]); - aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iOpposSame ], prevNod[ iAfterSame ], - nextNod[ iBeforeSame ], nextNod[ iOpposSame ], nextNod[ iAfterSame ] ); - } - else if ( nbSame == 2 ) // 1 pentahedron + else if ( nbSame == 2 ) // pentahedron { if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] ) // iBeforeSame is same too - aMesh->AddVolume (prevNod[ iOpposSame ], prevNod[ iBeforeSame ], nextNod[ iOpposSame ], - prevNod[ iAfterSame ], prevNod[ iSameNode ], nextNod[ iAfterSame ]); + aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iOpposSame ], + nextNod[ iOpposSame ], prevNod[ iSameNode ], + prevNod[ iAfterSame ], nextNod[ iAfterSame ]); else // iAfterSame is same too - aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iSameNode ], nextNod[ iBeforeSame ], - prevNod[ iOpposSame ], prevNod[ iAfterSame ], nextNod[ iOpposSame ]); + aNewElem = aMesh->AddVolume (prevNod[ iSameNode ], prevNod[ iBeforeSame ], + nextNod[ iBeforeSame ], prevNod[ iAfterSame ], + prevNod[ iOpposSame ], nextNod[ iOpposSame ]); } break; } - default: - return; + default: { + // realized for extrusion only + vector polyedre_nodes (nbNodes*2 + 4*nbNodes); + vector quantities (nbNodes + 2); + + quantities[0] = nbNodes; // bottom of prism + for (int inode = 0; inode < nbNodes; inode++) { + polyedre_nodes[inode] = prevNod[inode]; + } + + quantities[1] = nbNodes; // top of prism + for (int inode = 0; inode < nbNodes; inode++) { + polyedre_nodes[nbNodes + inode] = nextNod[inode]; + } + + for (int iface = 0; iface < nbNodes; iface++) { + quantities[iface + 2] = 4; + int inextface = (iface == nbNodes - 1) ? 0 : iface + 1; + polyedre_nodes[2*nbNodes + 4*iface + 0] = prevNod[iface]; + polyedre_nodes[2*nbNodes + 4*iface + 1] = prevNod[inextface]; + polyedre_nodes[2*nbNodes + 4*iface + 2] = nextNod[inextface]; + polyedre_nodes[2*nbNodes + 4*iface + 3] = nextNod[iface]; + } + aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities); } + } + if ( aNewElem ) + newElems.push_back( aNewElem ); // set new prev nodes for ( iNode = 0; iNode < nbNodes; iNode++ ) @@ -1641,9 +2163,186 @@ static void sweepElement(SMESHDS_Mesh* aMesh, } // for steps } +//======================================================================= +//function : makeWalls +//purpose : create 1D and 2D elements around swept elements +//======================================================================= + +static void makeWalls (SMESHDS_Mesh* aMesh, + TNodeOfNodeListMap & mapNewNodes, + TElemOfElemListMap & newElemsMap, + TElemOfVecOfNnlmiMap & elemNewNodesMap, + set& elemSet) +{ + ASSERT( newElemsMap.size() == elemNewNodesMap.size() ); + + // Find nodes belonging to only one initial element - sweep them to get edges. + + TNodeOfNodeListMapItr nList = mapNewNodes.begin(); + for ( ; nList != mapNewNodes.end(); nList++ ) + { + const SMDS_MeshNode* node = + static_cast( nList->first ); + SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(); + int nbInitElems = 0; + while ( eIt->more() && nbInitElems < 2 ) + if ( elemSet.find( eIt->next() ) != elemSet.end() ) + nbInitElems++; + if ( nbInitElems < 2 ) { + vector newNodesItVec( 1, nList ); + list newEdges; + sweepElement( aMesh, node, newNodesItVec, newEdges ); + } + } + + // Make a ceiling for each element ie an equal element of last new nodes. + // Find free links of faces - make edges and sweep them into faces. + + TElemOfElemListMap::iterator itElem = newElemsMap.begin(); + TElemOfVecOfNnlmiMap::iterator itElemNodes = elemNewNodesMap.begin(); + for ( ; itElem != newElemsMap.end(); itElem++, itElemNodes++ ) + { + const SMDS_MeshElement* elem = itElem->first; + vector& vecNewNodes = itElemNodes->second; + + if ( elem->GetType() == SMDSAbs_Edge ) + { + // create a ceiling edge + aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(), + vecNewNodes[ 1 ]->second.back() ); + } + if ( elem->GetType() != SMDSAbs_Face ) + continue; + + bool hasFreeLinks = false; + + set avoidSet; + avoidSet.insert( elem ); + + // loop on a face nodes + set aFaceLastNodes; + int iNode, nbNodes = vecNewNodes.size(); + for ( iNode = 0; iNode < nbNodes; iNode++ ) + { + aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() ); + // look for free links of a face + int iNext = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1; + const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first; + const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first; + // check if a link is free + if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) + { + hasFreeLinks = true; + // make an edge and a ceiling for a new edge + if ( !aMesh->FindEdge( n1, n2 )) + aMesh->AddEdge( n1, n2 ); + n1 = vecNewNodes[ iNode ]->second.back(); + n2 = vecNewNodes[ iNext ]->second.back(); + if ( !aMesh->FindEdge( n1, n2 )) + aMesh->AddEdge( n1, n2 ); + } + } + // sweep free links into faces + + if ( hasFreeLinks ) + { + list & newVolumes = itElem->second; + int iStep, nbSteps = vecNewNodes[0]->second.size(); + int iVol, volNb, nbVolumesByStep = newVolumes.size() / nbSteps; + + set initNodeSet, faceNodeSet; + for ( iNode = 0; iNode < nbNodes; iNode++ ) + initNodeSet.insert( vecNewNodes[ iNode ]->first ); + + for ( volNb = 0; volNb < nbVolumesByStep; volNb++ ) + { + list::iterator v = newVolumes.begin(); + iVol = 0; + while ( iVol++ < volNb ) v++; + // find indices of free faces of a volume + list< int > fInd; + SMDS_VolumeTool vTool( *v ); + int iF, nbF = vTool.NbFaces(); + for ( iF = 0; iF < nbF; iF ++ ) + if (vTool.IsFreeFace( iF ) && + vTool.GetFaceNodes( iF, faceNodeSet ) && + initNodeSet != faceNodeSet) // except an initial face + fInd.push_back( iF ); + if ( fInd.empty() ) + continue; + + // create faces for all steps + for ( iStep = 0; iStep < nbSteps; iStep++ ) + { + vTool.Set( *v ); + vTool.SetExternalNormal(); + list< int >::iterator ind = fInd.begin(); + for ( ; ind != fInd.end(); ind++ ) + { + const SMDS_MeshNode** nodes = vTool.GetFaceNodes( *ind ); + switch ( vTool.NbFaceNodes( *ind ) ) { + case 3: + aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ); break; + case 4: + aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ); break; + default: + { + int nbPolygonNodes = vTool.NbFaceNodes( *ind ); + vector polygon_nodes (nbPolygonNodes); + for (int inode = 0; inode < nbPolygonNodes; inode++) { + polygon_nodes[inode] = nodes[inode]; + } + aMesh->AddPolygonalFace(polygon_nodes); + break; + } + } + } + // go to the next volume + iVol = 0; + while ( iVol++ < nbVolumesByStep ) v++; + } + } + } // sweep free links into faces + + // make a ceiling face with a normal external to a volume + + SMDS_VolumeTool lastVol( itElem->second.back() ); + int iF = lastVol.GetFaceIndex( aFaceLastNodes ); + if ( iF >= 0 ) + { + lastVol.SetExternalNormal(); + const SMDS_MeshNode** nodes = lastVol.GetFaceNodes( iF ); + switch ( lastVol.NbFaceNodes( iF ) ) { + case 3: + if (!hasFreeLinks || + !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ])) + aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ); + break; + case 4: + if (!hasFreeLinks || + !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ])) + aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ); + break; + default: + { + int nbPolygonNodes = lastVol.NbFaceNodes( iF ); + vector polygon_nodes (nbPolygonNodes); + for (int inode = 0; inode < nbPolygonNodes; inode++) { + polygon_nodes[inode] = nodes[inode]; + } + if (!hasFreeLinks || !aMesh->FindFace(polygon_nodes)) + aMesh->AddPolygonalFace(polygon_nodes); + } + break; + } + } + + } // loop on swept elements +} + //======================================================================= //function : RotationSweep -//purpose : +//purpose : //======================================================================= void SMESH_MeshEditor::RotationSweep(set & theElems, @@ -1652,6 +2351,7 @@ void SMESH_MeshEditor::RotationSweep(set & theElems, const int theNbSteps, const double theTol) { + MESSAGE( "RotationSweep()"); gp_Trsf aTrsf; aTrsf.SetRotation( theAxis, theAngle ); @@ -1661,6 +2361,72 @@ void SMESH_MeshEditor::RotationSweep(set & theElems, SMESHDS_Mesh* aMesh = GetMeshDS(); TNodeOfNodeListMap mapNewNodes; + TElemOfVecOfNnlmiMap mapElemNewNodes; + TElemOfElemListMap newElemsMap; + + // loop on theElems + set< const SMDS_MeshElement* >::iterator itElem; + for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) + { + const SMDS_MeshElement* elem = (*itElem); + if ( !elem ) + continue; + vector & newNodesItVec = mapElemNewNodes[ elem ]; + newNodesItVec.reserve( elem->NbNodes() ); + + // loop on elem nodes + SMDS_ElemIteratorPtr itN = elem->nodesIterator(); + while ( itN->more() ) { + + // check if a node has been already sweeped + const SMDS_MeshNode* node = + static_cast( itN->next() ); + TNodeOfNodeListMapItr nIt = mapNewNodes.find( node ); + if ( nIt == mapNewNodes.end() ) + { + nIt = mapNewNodes.insert( make_pair( node, list() )).first; + list& listNewNodes = nIt->second; + + // make new nodes + gp_XYZ aXYZ( node->X(), node->Y(), node->Z() ); + double coord[3]; + aXYZ.Coord( coord[0], coord[1], coord[2] ); + bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol ); + const SMDS_MeshNode * newNode = node; + for ( int i = 0; i < theNbSteps; i++ ) { + if ( !isOnAxis ) { + aTrsf.Transforms( coord[0], coord[1], coord[2] ); + newNode = aMesh->AddNode( coord[0], coord[1], coord[2] ); + } + listNewNodes.push_back( newNode ); + } + } + newNodesItVec.push_back( nIt ); + } + // make new elements + sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem] ); + } + + makeWalls( aMesh, mapNewNodes, newElemsMap, mapElemNewNodes, theElems ); + +} +//======================================================================= +//function : ExtrusionSweep +//purpose : +//======================================================================= + +void SMESH_MeshEditor::ExtrusionSweep(set & theElems, + const gp_Vec& theStep, + const int theNbSteps) +{ + gp_Trsf aTrsf; + aTrsf.SetTranslation( theStep ); + + SMESHDS_Mesh* aMesh = GetMeshDS(); + + TNodeOfNodeListMap mapNewNodes; + TElemOfVecOfNnlmiMap mapElemNewNodes; + TElemOfElemListMap newElemsMap; // loop on theElems set< const SMDS_MeshElement* >::iterator itElem; @@ -1668,11 +2434,12 @@ void SMESH_MeshEditor::RotationSweep(set & theElems, { // check element type const SMDS_MeshElement* elem = (*itElem); - if ( !elem || - (elem->GetType() != SMDSAbs_Face && - elem->GetType() != SMDSAbs_Edge )) + if ( !elem ) continue; + vector & newNodesItVec = mapElemNewNodes[ elem ]; + newNodesItVec.reserve( elem->NbNodes() ); + // loop on elem nodes SMDS_ElemIteratorPtr itN = elem->nodesIterator(); while ( itN->more() ) { @@ -1680,87 +2447,369 @@ void SMESH_MeshEditor::RotationSweep(set & theElems, // check if a node has been already sweeped const SMDS_MeshNode* node = static_cast( itN->next() ); - if (mapNewNodes.find( node ) != mapNewNodes.end() ) - continue; - - list& listNewNodes = mapNewNodes[ node ]; + TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node ); + if ( nIt == mapNewNodes.end() ) + { + nIt = mapNewNodes.insert( make_pair( node, list() )).first; + list& listNewNodes = nIt->second; - // make new nodes - gp_XYZ aXYZ( node->X(), node->Y(), node->Z() ); - double coord[3]; - aXYZ.Coord( coord[0], coord[1], coord[2] ); - bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol ); - const SMDS_MeshNode * newNode = node; - for ( int i = 0; i < theNbSteps; i++ ) { - if ( !isOnAxis ) { + // make new nodes + double coord[] = { node->X(), node->Y(), node->Z() }; + for ( int i = 0; i < theNbSteps; i++ ) { aTrsf.Transforms( coord[0], coord[1], coord[2] ); - newNode = aMesh->AddNode( coord[0], coord[1], coord[2] ); + const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] ); + listNewNodes.push_back( newNode ); } - listNewNodes.push_back( newNode ); } + newNodesItVec.push_back( nIt ); } // make new elements - sweepElement( aMesh, elem, mapNewNodes ); + sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem] ); + + // fill history + TColStd_ListOfInteger ListNewID; + list tmpList = newElemsMap[elem]; + for(list::iterator ite = tmpList.begin(); + ite!=tmpList.end(); ite++) { + ListNewID.Append((*ite)->GetID()); + } + myExtrusionHistory.Bind(elem->GetID(),ListNewID); + // end fill history + } + makeWalls( aMesh, mapNewNodes, newElemsMap, mapElemNewNodes, theElems ); } + //======================================================================= -//function : ExtrusionSweep -//purpose : +//class : SMESH_MeshEditor_PathPoint +//purpose : auxiliary class //======================================================================= +class SMESH_MeshEditor_PathPoint { +public: + SMESH_MeshEditor_PathPoint() { + myPnt.SetCoord(99., 99., 99.); + myTgt.SetCoord(1.,0.,0.); + myAngle=0.; + myPrm=0.; + } + void SetPnt(const gp_Pnt& aP3D){ + myPnt=aP3D; + } + void SetTangent(const gp_Dir& aTgt){ + myTgt=aTgt; + } + void SetAngle(const double& aBeta){ + myAngle=aBeta; + } + void SetParameter(const double& aPrm){ + myPrm=aPrm; + } + const gp_Pnt& Pnt()const{ + return myPnt; + } + const gp_Dir& Tangent()const{ + return myTgt; + } + double Angle()const{ + return myAngle; + } + double Parameter()const{ + return myPrm; + } -void SMESH_MeshEditor::ExtrusionSweep(set & theElems, - const gp_Vec& theStep, - const int theNbSteps) -{ - gp_Trsf aTrsf; - aTrsf.SetTranslation( theStep ); +protected: + gp_Pnt myPnt; + gp_Dir myTgt; + double myAngle; + double myPrm; +}; - SMESHDS_Mesh* aMesh = GetMeshDS(); +//======================================================================= +//function : ExtrusionAlongTrack +//purpose : +//======================================================================= +SMESH_MeshEditor::Extrusion_Error + SMESH_MeshEditor::ExtrusionAlongTrack (std::set & theElements, + SMESH_subMesh* theTrack, + const SMDS_MeshNode* theN1, + const bool theHasAngles, + std::list& theAngles, + const bool theHasRefPoint, + const gp_Pnt& theRefPoint) +{ + MESSAGE("SMESH_MeshEditor::ExtrusionAlongTrack") + int j, aNbTP, aNbE, aNb; + double aT1, aT2, aT, aAngle, aX, aY, aZ; + std::list aPrms; + std::list::iterator aItD; + std::set< const SMDS_MeshElement* >::iterator itElem; + + Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2; + gp_Pnt aP3D, aV0; + gp_Vec aVec; + gp_XYZ aGC; + Handle(Geom_Curve) aC3D; + TopoDS_Edge aTrackEdge; + TopoDS_Vertex aV1, aV2; + + SMDS_ElemIteratorPtr aItE; + SMDS_NodeIteratorPtr aItN; + SMDSAbs_ElementType aTypeE; TNodeOfNodeListMap mapNewNodes; + TElemOfVecOfNnlmiMap mapElemNewNodes; + TElemOfElemListMap newElemsMap; + + aTolVec=1.e-7; + aTolVec2=aTolVec*aTolVec; + + // 1. Check data + aNbE = theElements.size(); + // nothing to do + if ( !aNbE ) + return EXTR_NO_ELEMENTS; + + // 1.1 Track Pattern + ASSERT( theTrack ); + + SMESHDS_SubMesh* pSubMeshDS=theTrack->GetSubMeshDS(); + + aItE = pSubMeshDS->GetElements(); + while ( aItE->more() ) { + const SMDS_MeshElement* pE = aItE->next(); + aTypeE = pE->GetType(); + // Pattern must contain links only + if ( aTypeE != SMDSAbs_Edge ) + return EXTR_PATH_NOT_EDGE; + } - // loop on theElems - set< const SMDS_MeshElement* >::iterator itElem; - for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) - { + const TopoDS_Shape& aS = theTrack->GetSubShape(); + // Sub shape for the Pattern must be an Edge + if ( aS.ShapeType() != TopAbs_EDGE ) + return EXTR_BAD_PATH_SHAPE; + + aTrackEdge = TopoDS::Edge( aS ); + // the Edge must not be degenerated + if ( BRep_Tool::Degenerated( aTrackEdge ) ) + return EXTR_BAD_PATH_SHAPE; + + TopExp::Vertices( aTrackEdge, aV1, aV2 ); + aT1=BRep_Tool::Parameter( aV1, aTrackEdge ); + aT2=BRep_Tool::Parameter( aV2, aTrackEdge ); + + aItN = theTrack->GetFather()->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes(); + const SMDS_MeshNode* aN1 = aItN->next(); + + aItN = theTrack->GetFather()->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes(); + const SMDS_MeshNode* aN2 = aItN->next(); + + // starting node must be aN1 or aN2 + if ( !( aN1 == theN1 || aN2 == theN1 ) ) + return EXTR_BAD_STARTING_NODE; + + aNbTP = pSubMeshDS->NbNodes() + 2; + + // 1.2. Angles + vector aAngles( aNbTP ); + + for ( j=0; j < aNbTP; ++j ) { + aAngles[j] = 0.; + } + + if ( theHasAngles ) { + aItD = theAngles.begin(); + for ( j=1; (aItD != theAngles.end()) && (jGetNodes(); + while ( aItN->more() ) { + const SMDS_MeshNode* pNode = aItN->next(); + const SMDS_EdgePosition* pEPos = + static_cast( pNode->GetPosition().get() ); + aT = pEPos->GetUParameter(); + aPrms.push_back( aT ); + } + + // sort parameters + aPrms.sort(); + if ( aN1 == theN1 ) { + if ( aT1 > aT2 ) { + aPrms.reverse(); + } + } + else { + if ( aT2 > aT1 ) { + aPrms.reverse(); + } + } + + // 3. Path Points + SMESH_MeshEditor_PathPoint aPP; + vector aPPs( aNbTP ); + // + aC3D = BRep_Tool::Curve( aTrackEdge, aTx1, aTx2 ); + // + aItD = aPrms.begin(); + for ( j=0; aItD != aPrms.end(); ++aItD, ++j ) { + aT = *aItD; + aC3D->D1( aT, aP3D, aVec ); + aL2 = aVec.SquareMagnitude(); + if ( aL2 < aTolVec2 ) + return EXTR_CANT_GET_TANGENT; + + gp_Dir aTgt( aVec ); + aAngle = aAngles[j]; + + aPP.SetPnt( aP3D ); + aPP.SetTangent( aTgt ); + aPP.SetAngle( aAngle ); + aPP.SetParameter( aT ); + aPPs[j]=aPP; + } + + // 3. Center of rotation aV0 + aV0 = theRefPoint; + if ( !theHasRefPoint ) { + aNb = 0; + aGC.SetCoord( 0.,0.,0. ); + + itElem = theElements.begin(); + for ( ; itElem != theElements.end(); itElem++ ) { + const SMDS_MeshElement* elem = (*itElem); + + SMDS_ElemIteratorPtr itN = elem->nodesIterator(); + while ( itN->more() ) { + const SMDS_MeshNode* node = static_cast( itN->next() ); + aX = node->X(); + aY = node->Y(); + aZ = node->Z(); + + if ( mapNewNodes.find( node ) == mapNewNodes.end() ) { + list aLNx; + mapNewNodes[node] = aLNx; + // + gp_XYZ aXYZ( aX, aY, aZ ); + aGC += aXYZ; + ++aNb; + } + } + } + aGC /= aNb; + aV0.SetXYZ( aGC ); + } // if (!theHasRefPoint) { + mapNewNodes.clear(); + + // 4. Processing the elements + SMESHDS_Mesh* aMesh = GetMeshDS(); + + for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) { // check element type const SMDS_MeshElement* elem = (*itElem); - if ( !elem || - (elem->GetType() != SMDSAbs_Face && - elem->GetType() != SMDSAbs_Edge)) + aTypeE = elem->GetType(); + if ( !elem || ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge ) ) continue; + vector & newNodesItVec = mapElemNewNodes[ elem ]; + newNodesItVec.reserve( elem->NbNodes() ); + // loop on elem nodes SMDS_ElemIteratorPtr itN = elem->nodesIterator(); while ( itN->more() ) { - // check if a node has been already sweeped + // check if a node has been already processed const SMDS_MeshNode* node = - static_cast( itN->next() ); - if (mapNewNodes.find( node ) != mapNewNodes.end() ) - continue; - - list& listNewNodes = mapNewNodes[ node ]; - - // make new nodes - double coord[3]; - coord[0] = node->X(); - coord[1] = node->Y(); - coord[2] = node->Z(); - for ( int i = 0; i < theNbSteps; i++ ) { - aTrsf.Transforms( coord[0], coord[1], coord[2] ); - const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] ); - listNewNodes.push_back( newNode ); + static_cast( itN->next() ); + TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node ); + if ( nIt == mapNewNodes.end() ) { + nIt = mapNewNodes.insert( make_pair( node, list() )).first; + list& listNewNodes = nIt->second; + + // make new nodes + aX = node->X(); aY = node->Y(); aZ = node->Z(); + + Standard_Real aAngle1x, aAngleT1T0, aTolAng; + gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x; + gp_Ax1 anAx1, anAxT1T0; + gp_Dir aDT1x, aDT0x, aDT1T0; + + aTolAng=1.e-4; + + aV0x = aV0; + aPN0.SetCoord(aX, aY, aZ); + + const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0]; + aP0x = aPP0.Pnt(); + aDT0x= aPP0.Tangent(); + + for ( j = 1; j < aNbTP; ++j ) { + const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j]; + aP1x = aPP1.Pnt(); + aDT1x = aPP1.Tangent(); + aAngle1x = aPP1.Angle(); + + gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0; + // Translation + gp_Vec aV01x( aP0x, aP1x ); + aTrsf.SetTranslation( aV01x ); + + // traslated point + aV1x = aV0x.Transformed( aTrsf ); + aPN1 = aPN0.Transformed( aTrsf ); + + // rotation 1 [ T1,T0 ] + aAngleT1T0=-aDT1x.Angle( aDT0x ); + if (fabs(aAngleT1T0) > aTolAng) { + aDT1T0=aDT1x^aDT0x; + anAxT1T0.SetLocation( aV1x ); + anAxT1T0.SetDirection( aDT1T0 ); + aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 ); + + aPN1 = aPN1.Transformed( aTrsfRotT1T0 ); + } + + // rotation 2 + if ( theHasAngles ) { + anAx1.SetLocation( aV1x ); + anAx1.SetDirection( aDT1x ); + aTrsfRot.SetRotation( anAx1, aAngle1x ); + + aPN1 = aPN1.Transformed( aTrsfRot ); + } + + // make new node + aX = aPN1.X(); + aY = aPN1.Y(); + aZ = aPN1.Z(); + const SMDS_MeshNode* newNode = aMesh->AddNode( aX, aY, aZ ); + listNewNodes.push_back( newNode ); + + aPN0 = aPN1; + aP0x = aP1x; + aV0x = aV1x; + aDT0x = aDT1x; + } } + newNodesItVec.push_back( nIt ); } // make new elements - sweepElement( aMesh, elem, mapNewNodes ); + sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem] ); } + + makeWalls( aMesh, mapNewNodes, newElemsMap, mapElemNewNodes, theElements ); + + return EXTR_OK; } //======================================================================= //function : Transform -//purpose : +//purpose : //======================================================================= void SMESH_MeshEditor::Transform (set & theElems, @@ -1798,11 +2847,11 @@ void SMESH_MeshEditor::Transform (set & theElems, SMDS_ElemIteratorPtr itN = elem->nodesIterator(); while ( itN->more() ) { - // check if a node has been already transormed + // check if a node has been already transformed const SMDS_MeshNode* node = static_cast( itN->next() ); if (nodeMap.find( node ) != nodeMap.end() ) - continue; + continue; double coord[3]; coord[0] = node->X(); @@ -1812,8 +2861,12 @@ void SMESH_MeshEditor::Transform (set & theElems, const SMDS_MeshNode * newNode = node; if ( theCopy ) newNode = aMesh->AddNode( coord[0], coord[1], coord[2] ); - else + else { aMesh->MoveNode( node, coord[0], coord[1], coord[2] ); + // node position on shape becomes invalid + const_cast< SMDS_MeshNode* > ( node )->SetPosition + ( SMDS_SpacePosition::originSpacePosition() ); + } nodeMap.insert( TNodeNodeMap::value_type( node, newNode )); // keep inverse elements @@ -1836,7 +2889,7 @@ void SMESH_MeshEditor::Transform (set & theElems, theElems.insert( *invElemIt ); } - // replicate or reverse elements + // replicate or reverse elements enum { REV_TETRA = 0, // = nbNodes - 4 @@ -1847,12 +2900,12 @@ void SMESH_MeshEditor::Transform (set & theElems, FORWARD = 5 }; int index[][8] = { - { 2, 1, 0, 3, 4, 0, 0, 0 }, // REV_TETRA + { 2, 1, 0, 3, 4, 0, 0, 0 }, // REV_TETRA { 2, 1, 0, 3, 4, 0, 0, 0 }, // REV_PYRAMID - { 2, 1, 0, 5, 4, 3, 0, 0 }, // REV_PENTA - { 2, 1, 0, 3, 0, 0, 0, 0 }, // REV_FACE - { 2, 1, 0, 3, 6, 5, 4, 7 }, // REV_HEXA - { 0, 1, 2, 3, 4, 5, 6, 7 } // FORWARD + { 2, 1, 0, 5, 4, 3, 0, 0 }, // REV_PENTA + { 2, 1, 0, 3, 0, 0, 0, 0 }, // REV_FACE + { 2, 1, 0, 3, 6, 5, 4, 7 }, // REV_HEXA + { 0, 1, 2, 3, 4, 5, 6, 7 } // FORWARD }; for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) @@ -1864,6 +2917,82 @@ void SMESH_MeshEditor::Transform (set & theElems, int nbNodes = elem->NbNodes(); int elemType = elem->GetType(); + if (elem->IsPoly()) { + // Polygon or Polyhedral Volume + switch ( elemType ) { + case SMDSAbs_Face: + { + vector poly_nodes (nbNodes); + int iNode = 0; + SMDS_ElemIteratorPtr itN = elem->nodesIterator(); + while (itN->more()) { + const SMDS_MeshNode* node = + static_cast(itN->next()); + TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node); + if (nodeMapIt == nodeMap.end()) + break; // not all nodes transformed + if (needReverse) { + // reverse mirrored faces and volumes + poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second; + } else { + poly_nodes[iNode] = (*nodeMapIt).second; + } + iNode++; + } + if ( iNode != nbNodes ) + continue; // not all nodes transformed + + if ( theCopy ) { + aMesh->AddPolygonalFace(poly_nodes); + } else { + aMesh->ChangePolygonNodes(elem, poly_nodes); + } + } + break; + case SMDSAbs_Volume: + { + // ATTENTION: Reversing is not yet done!!! + const SMDS_PolyhedralVolumeOfNodes* aPolyedre = + (const SMDS_PolyhedralVolumeOfNodes*) elem; + if (!aPolyedre) { + MESSAGE("Warning: bad volumic element"); + continue; + } + + vector poly_nodes; + vector quantities; + + bool allTransformed = true; + int nbFaces = aPolyedre->NbFaces(); + for (int iface = 1; iface <= nbFaces && allTransformed; iface++) { + int nbFaceNodes = aPolyedre->NbFaceNodes(iface); + for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) { + const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode); + TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node); + if (nodeMapIt == nodeMap.end()) { + allTransformed = false; // not all nodes transformed + } else { + poly_nodes.push_back((*nodeMapIt).second); + } + } + quantities.push_back(nbFaceNodes); + } + if ( !allTransformed ) + continue; // not all nodes transformed + + if ( theCopy ) { + aMesh->AddPolyhedralVolume(poly_nodes, quantities); + } else { + aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities); + } + } + break; + default:; + } + continue; + } + + // Regular elements int* i = index[ FORWARD ]; if ( needReverse && nbNodes > 2) // reverse mirrored faces and volumes if ( elemType == SMDSAbs_Face ) @@ -1887,7 +3016,7 @@ void SMESH_MeshEditor::Transform (set & theElems, if ( iNode != nbNodes ) continue; // not all nodes transformed - if ( theCopy ) + if ( theCopy ) { // add a new element switch ( elemType ) { @@ -1928,17 +3057,26 @@ void SMESH_MeshEditor::Transform (set & theElems, //======================================================================= //function : FindCoincidentNodes //purpose : Return list of group of nodes close to each other within theTolerance +// Search among theNodes or in the whole mesh if theNodes is empty. //======================================================================= -void SMESH_MeshEditor::FindCoincidentNodes (const double theTolerance, - TListOfListOfNodes & theGroupsOfNodes) +void SMESH_MeshEditor::FindCoincidentNodes (set & theNodes, + const double theTolerance, + TListOfListOfNodes & theGroupsOfNodes) { double tol2 = theTolerance * theTolerance; list nodes; - SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator(); - while ( nIt->more() ) - nodes.push_back( nIt->next() ); + if ( theNodes.empty() ) + { // get all nodes in the mesh + SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator(); + while ( nIt->more() ) + nodes.push_back( nIt->next() ); + } + else + { + nodes.insert( nodes.end(), theNodes.begin(), theNodes.end() ); + } list::iterator it2, it1 = nodes.begin(); for ( ; it1 != nodes.end(); it1++ ) @@ -1967,6 +3105,88 @@ void SMESH_MeshEditor::FindCoincidentNodes (const double theTolerance, } } +//======================================================================= +//function : SimplifyFace +//purpose : +//======================================================================= +int SMESH_MeshEditor::SimplifyFace (const vector faceNodes, + vector& poly_nodes, + vector& quantities) const +{ + int nbNodes = faceNodes.size(); + + if (nbNodes < 3) + return 0; + + set nodeSet; + + // get simple seq of nodes + const SMDS_MeshNode* simpleNodes[ nbNodes ]; + int iSimple = 0, nbUnique = 0; + + simpleNodes[iSimple++] = faceNodes[0]; + nbUnique++; + for (int iCur = 1; iCur < nbNodes; iCur++) { + if (faceNodes[iCur] != simpleNodes[iSimple - 1]) { + simpleNodes[iSimple++] = faceNodes[iCur]; + if (nodeSet.insert( faceNodes[iCur] ).second) + nbUnique++; + } + } + int nbSimple = iSimple; + if (simpleNodes[nbSimple - 1] == simpleNodes[0]) { + nbSimple--; + iSimple--; + } + + if (nbUnique < 3) + return 0; + + // separate loops + int nbNew = 0; + bool foundLoop = (nbSimple > nbUnique); + while (foundLoop) { + foundLoop = false; + set loopSet; + for (iSimple = 0; iSimple < nbSimple && !foundLoop; iSimple++) { + const SMDS_MeshNode* n = simpleNodes[iSimple]; + if (!loopSet.insert( n ).second) { + foundLoop = true; + + // separate loop + int iC = 0, curLast = iSimple; + for (; iC < curLast; iC++) { + if (simpleNodes[iC] == n) break; + } + int loopLen = curLast - iC; + if (loopLen > 2) { + // create sub-element + nbNew++; + quantities.push_back(loopLen); + for (; iC < curLast; iC++) { + poly_nodes.push_back(simpleNodes[iC]); + } + } + // shift the rest nodes (place from the first loop position) + for (iC = curLast + 1; iC < nbSimple; iC++) { + simpleNodes[iC - loopLen] = simpleNodes[iC]; + } + nbSimple -= loopLen; + iSimple -= loopLen; + } + } // for (iSimple = 0; iSimple < nbSimple; iSimple++) + } // while (foundLoop) + + if (iSimple > 2) { + nbNew++; + quantities.push_back(iSimple); + for (int i = 0; i < iSimple; i++) + poly_nodes.push_back(simpleNodes[i]); + } + + return nbNew; +} + //======================================================================= //function : MergeNodes //purpose : In each group, the cdr of nodes are substituted by the first one @@ -1995,7 +3215,7 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes) nodeNodeMap.insert( TNodeNodeMap::value_type( nToRemove, nToKeep )); if ( nToRemove != nToKeep ) { rmNodeIds.push_back( nToRemove->GetID() ); - addToSameGroups( nToKeep, nToRemove, aMesh ); + AddToSameGroups( nToKeep, nToRemove, aMesh ); } SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator(); @@ -2003,7 +3223,7 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes) elems.insert( invElemIt->next() ); } } - // Change element nodes or remove an element + // Change element nodes or remove an element set::iterator eIt = elems.begin(); for ( ; eIt != elems.end(); eIt++ ) @@ -2041,6 +3261,88 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes) int nbUniqueNodes = nodeSet.size(); if ( nbNodes != nbUniqueNodes ) // some nodes stick { + // Polygons and Polyhedral volumes + if (elem->IsPoly()) { + + if (elem->GetType() == SMDSAbs_Face) { + // Polygon + vector face_nodes (nbNodes); + int inode = 0; + for (; inode < nbNodes; inode++) { + face_nodes[inode] = curNodes[inode]; + } + + vector polygons_nodes; + vector quantities; + int nbNew = SimplifyFace(face_nodes, polygons_nodes, quantities); + + if (nbNew > 0) { + inode = 0; + for (int iface = 0; iface < nbNew - 1; iface++) { + int nbNodes = quantities[iface]; + vector poly_nodes (nbNodes); + for (int ii = 0; ii < nbNodes; ii++, inode++) { + poly_nodes[ii] = polygons_nodes[inode]; + } + SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes); + if (aShapeId) + aMesh->SetMeshElementOnShape(newElem, aShapeId); + } + aMesh->ChangeElementNodes(elem, &polygons_nodes[inode], quantities[nbNew - 1]); + } else { + rmElemIds.push_back(elem->GetID()); + } + + } else if (elem->GetType() == SMDSAbs_Volume) { + // Polyhedral volume + if (nbUniqueNodes < 4) { + rmElemIds.push_back(elem->GetID()); + } else { + // each face has to be analized in order to check volume validity + const SMDS_PolyhedralVolumeOfNodes* aPolyedre = + static_cast( elem ); + if (aPolyedre) { + int nbFaces = aPolyedre->NbFaces(); + + vector poly_nodes; + vector quantities; + + for (int iface = 1; iface <= nbFaces; iface++) { + int nbFaceNodes = aPolyedre->NbFaceNodes(iface); + vector faceNodes (nbFaceNodes); + + for (int inode = 1; inode <= nbFaceNodes; inode++) { + const SMDS_MeshNode * faceNode = aPolyedre->GetFaceNode(iface, inode); + TNodeNodeMap::iterator nnIt = nodeNodeMap.find(faceNode); + if (nnIt != nodeNodeMap.end()) { // faceNode sticks + faceNode = (*nnIt).second; + } + faceNodes[inode - 1] = faceNode; + } + + SimplifyFace(faceNodes, poly_nodes, quantities); + } + + if (quantities.size() > 3) { + // to be done: remove coincident faces + } + + if (quantities.size() > 3) + aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities); + else + rmElemIds.push_back(elem->GetID()); + + } else { + rmElemIds.push_back(elem->GetID()); + } + } + } else { + } + + continue; + } + + // Regular elements switch ( nbNodes ) { case 2: ///////////////////////////////////// EDGE isOk = false; break; @@ -2087,7 +3389,7 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes) if ( nbRepl == 2 && iRepl[ 1 ] - iRepl [ 0 ] == 3 ) { // a bottom node sticks with a linked top one // 1. - SMDS_MeshElement* newElem = + SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ 3 ], curNodes[ 4 ], curNodes[ 5 ], @@ -2163,7 +3465,7 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes) uniqueNodes[ 3 ] = curNodes[indTop[ 0 ]]; nbUniqueNodes = 4; // tetrahedron 2 - SMDS_MeshElement* newElem = + SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ind[ 0 ]], curNodes[ind[ 3 ]], curNodes[ind[ 2 ]], @@ -2281,11 +3583,42 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes) } // switch ( nbNodes ) } // if ( nbNodes != nbUniqueNodes ) // some nodes stick - - if ( isOk ) - aMesh->ChangeElementNodes( elem, uniqueNodes, nbUniqueNodes ); - else + + if ( isOk ) { + if (elem->IsPoly() && elem->GetType() == SMDSAbs_Volume) { + // Change nodes of polyedre + const SMDS_PolyhedralVolumeOfNodes* aPolyedre = + static_cast( elem ); + if (aPolyedre) { + int nbFaces = aPolyedre->NbFaces(); + + vector poly_nodes; + vector quantities (nbFaces); + + for (int iface = 1; iface <= nbFaces; iface++) { + int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface); + quantities[iface - 1] = nbFaceNodes; + + for (inode = 1; inode <= nbFaceNodes; inode++) { + const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode); + + TNodeNodeMap::iterator nnIt = nodeNodeMap.find( curNode ); + if (nnIt != nodeNodeMap.end()) { // curNode sticks + curNode = (*nnIt).second; + } + poly_nodes.push_back(curNode); + } + } + aMesh->ChangePolyhedronNodes( elem, poly_nodes, quantities ); + } + } else { + // Change regular element or polygon + aMesh->ChangeElementNodes( elem, uniqueNodes, nbUniqueNodes ); + } + } else { + // Remove invalid regular element or invalid polygon rmElemIds.push_back( elem->GetID() ); + } } // loop on elements @@ -2344,45 +3677,67 @@ void SMESH_MeshEditor::MergeEqualElements() } //======================================================================= -//function : findAdjacentFace -//purpose : +//function : FindFaceInSet +//purpose : Return a face having linked nodes n1 and n2 and which is +// - not in avoidSet, +// - in elemSet provided that !elemSet.empty() //======================================================================= -#define CHECKIND(max,val) {if ( (val) >= (max) ) \ -static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1, - const SMDS_MeshNode* n2, - const SMDS_MeshElement* elem) +const SMDS_MeshElement* + SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode* n1, + const SMDS_MeshNode* n2, + const set& elemSet, + const set& avoidSet) + { - SMDS_ElemIteratorPtr invElemIt = n1->facesIterator(); + SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(); while ( invElemIt->more() ) { // loop on inverse elements of n1 - const SMDS_MeshElement* adjElem = invElemIt->next(); - if ( elem != adjElem ) { - // get face nodes and find index of n1 - int i1, nbN = adjElem->NbNodes(), iNode = 0; - const SMDS_MeshNode* faceNodes[ nbN ], *n; - SMDS_ElemIteratorPtr nIt = adjElem->nodesIterator(); - while ( nIt->more() ) { - faceNodes[ iNode ] = static_cast( nIt->next() ); - if ( faceNodes[ iNode++ ] == n1 ) - i1 = iNode - 1; - } - // find a n2 linked to n1 - for ( iNode = 0; iNode < 2; iNode++ ) { - if ( iNode ) // node before n1 - n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ]; - else // node after n1 - n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ]; - if ( n == n2 ) - return adjElem; - } + const SMDS_MeshElement* elem = invElemIt->next(); + if (elem->GetType() != SMDSAbs_Face || + avoidSet.find( elem ) != avoidSet.end() ) + continue; + if ( !elemSet.empty() && elemSet.find( elem ) == elemSet.end()) + continue; + // get face nodes and find index of n1 + int i1, nbN = elem->NbNodes(), iNode = 0; + const SMDS_MeshNode* faceNodes[ nbN ], *n; + SMDS_ElemIteratorPtr nIt = elem->nodesIterator(); + while ( nIt->more() ) { + faceNodes[ iNode ] = static_cast( nIt->next() ); + if ( faceNodes[ iNode++ ] == n1 ) + i1 = iNode - 1; + } + // find a n2 linked to n1 + for ( iNode = 0; iNode < 2; iNode++ ) { + if ( iNode ) // node before n1 + n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ]; + else // node after n1 + n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ]; + if ( n == n2 ) + return elem; } } return 0; } - + +//======================================================================= +//function : findAdjacentFace +//purpose : +//======================================================================= + +static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1, + const SMDS_MeshNode* n2, + const SMDS_MeshElement* elem) +{ + set elemSet, avoidSet; + if ( elem ) + avoidSet.insert ( elem ); + return SMESH_MeshEditor::FindFaceInSet( n1, n2, elemSet, avoidSet ); +} + //======================================================================= //function : findFreeBorder -//purpose : +//purpose : //======================================================================= #define ControlFreeBorder SMESH::Controls::FreeEdges::IsFreeEdge @@ -2528,7 +3883,7 @@ bool SMESH_MeshEditor::CheckFreeBorderNodes(const SMDS_MeshNode* theNode1, //======================================================================= //function : SewFreeBorder -//purpose : +//purpose : //======================================================================= SMESH_MeshEditor::Sew_Error @@ -2538,7 +3893,9 @@ SMESH_MeshEditor::Sew_Error const SMDS_MeshNode* theSideFirstNode, const SMDS_MeshNode* theSideSecondNode, const SMDS_MeshNode* theSideThirdNode, - bool theSideIsFreeBorder) + const bool theSideIsFreeBorder, + const bool toCreatePolygons, + const bool toCreatePolyedrs) { MESSAGE("::SewFreeBorder()"); Sew_Error aResult = SEW_OK; @@ -2560,7 +3917,7 @@ SMESH_MeshEditor::Sew_Error aResult = SEW_BORDER1_NOT_FOUND; } if (theSideIsFreeBorder) - { + { // Free border 2 // -------------- if (!findFreeBorder(theSideFirstNode, theSideSecondNode, theSideThirdNode, @@ -2619,7 +3976,7 @@ SMESH_MeshEditor::Sew_Error toBordSys.SetTransformation( toBordAx ); fromSide2Sys.SetTransformation( fromSideAx, toGlobalAx ); fromSide2Sys.SetScaleFactor( Zs.Magnitude() / Zb.Magnitude() ); - + // move for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) { const SMDS_MeshNode* n = *nBordIt; @@ -2649,7 +4006,7 @@ SMESH_MeshEditor::Sew_Error LinkID_Gen aLinkID_Gen( GetMeshDS() ); set foundSideLinkIDs, checkedLinkIDs; SMDS_VolumeTool volume; - const SMDS_MeshNode* faceNodes[ 4 ]; + //const SMDS_MeshNode* faceNodes[ 4 ]; const SMDS_MeshNode* sideNode; const SMDS_MeshElement* sideElem; @@ -2677,6 +4034,7 @@ SMESH_MeshEditor::Sew_Error const SMDS_MeshElement* elem = invElemIt->next(); // prepare data for a loop on links, of a face or a volume int iPrevNode, iNode = 0, nbNodes = elem->NbNodes(); + const SMDS_MeshNode* faceNodes[ nbNodes ]; bool isVolume = volume.Set( elem ); const SMDS_MeshNode** nodes = isVolume ? volume.GetNodes() : faceNodes; if ( isVolume ) // --volume @@ -2760,7 +4118,7 @@ SMESH_MeshEditor::Sew_Error } while ( sideNode != theSideSecondNode ); - if ( hasVolumes && sideNodes.size () != bordNodes.size() ) { + if ( hasVolumes && sideNodes.size () != bordNodes.size() && !toCreatePolyedrs) { MESSAGE("VOLUME SPLITTING IS FORBIDDEN"); return SEW_VOLUMES_TO_SPLIT; // volume splitting is forbidden } @@ -2784,7 +4142,7 @@ SMESH_MeshEditor::Sew_Error nIt[0]++, nIt[1]++ ) { nodeGroupsToMerge.push_back( list() ); - nodeGroupsToMerge.back().push_back( *nIt[1] ); // to keep + nodeGroupsToMerge.back().push_back( *nIt[1] ); // to keep nodeGroupsToMerge.back().push_back( *nIt[0] ); // tp remove } } @@ -2841,7 +4199,7 @@ SMESH_MeshEditor::Sew_Error double minParam = Min( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]); double maxParam = Max( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]); double minSegLen = Min( nextParam - minParam, maxParam - prevParam ); - + // choose to insert or to merge nodes double du = param[ 1 ][ i[1] ] - param[ 0 ][ i[0] ]; if ( Abs( du ) <= minSegLen * 0.2 ) { @@ -2884,15 +4242,19 @@ SMESH_MeshEditor::Sew_Error list & nodeList = (*insertMapIt).second; const SMDS_MeshNode* n12 = nodeList.front(); nodeList.pop_front(); const SMDS_MeshNode* n22 = nodeList.front(); nodeList.pop_front(); - InsertNodesIntoLink( elem, n12, n22, nodeList ); + InsertNodesIntoLink( elem, n12, n22, nodeList, toCreatePolygons ); // 2. perform insertion into the link of adjacent faces while (true) { const SMDS_MeshElement* adjElem = findAdjacentFace( n12, n22, elem ); if ( adjElem ) - InsertNodesIntoLink( adjElem, n12, n22, nodeList ); + InsertNodesIntoLink( adjElem, n12, n22, nodeList, toCreatePolygons ); else break; } + if (toCreatePolyedrs) { + // perform insertion into the links of adjacent volumes + UpdateVolumes(n12, n22, nodeList); + } // 3. find an element appeared on n1 and n2 after the insertion insertMap.erase( elem ); elem = findAdjacentFace( n1, n2, 0 ); @@ -2932,18 +4294,22 @@ SMESH_MeshEditor::Sew_Error const SMDS_MeshNode* n1 = nodeList.front(); nodeList.pop_front(); const SMDS_MeshNode* n2 = nodeList.front(); nodeList.pop_front(); - InsertNodesIntoLink( elem, n1, n2, nodeList ); + InsertNodesIntoLink( elem, n1, n2, nodeList, toCreatePolygons ); if ( !theSideIsFreeBorder ) { // look for and insert nodes into the faces adjacent to elem while (true) { const SMDS_MeshElement* adjElem = findAdjacentFace( n1, n2, elem ); if ( adjElem ) - InsertNodesIntoLink( adjElem, n1, n2, nodeList ); + InsertNodesIntoLink( adjElem, n1, n2, nodeList, toCreatePolygons ); else break; } } + if (toCreatePolyedrs) { + // perform insertion into the links of adjacent volumes + UpdateVolumes(n1, n2, nodeList); + } } } // end: insert new nodes @@ -2962,14 +4328,15 @@ SMESH_MeshEditor::Sew_Error void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theFace, const SMDS_MeshNode* theBetweenNode1, const SMDS_MeshNode* theBetweenNode2, - list& theNodesToInsert) + list& theNodesToInsert, + const bool toCreatePoly) { if ( theFace->GetType() != SMDSAbs_Face ) return; // find indices of 2 link nodes and of the rest nodes int iNode = 0, il1, il2, i3, i4; il1 = il2 = i3 = i4 = -1; - const SMDS_MeshNode* nodes[ 8 ]; + const SMDS_MeshNode* nodes[ theFace->NbNodes() ]; SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator(); while ( nodeIt->more() ) { const SMDS_MeshNode* n = static_cast( nodeIt->next() ); @@ -2988,11 +4355,12 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theFace, // arrange link nodes to go one after another regarding the face orientation bool reverse = ( Abs( il2 - il1 ) == 1 ? il2 < il1 : il1 < il2 ); + list aNodesToInsert = theNodesToInsert; if ( reverse ) { iNode = il1; il1 = il2; il2 = iNode; - theNodesToInsert.reverse(); + aNodesToInsert.reverse(); } // check that not link nodes of a quadrangles are in good order int nbFaceNodes = theFace->NbNodes(); @@ -3000,15 +4368,61 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theFace, iNode = i3; i3 = i4; i4 = iNode; - } + } + + if (toCreatePoly || theFace->IsPoly()) { + + iNode = 0; + vector poly_nodes (nbFaceNodes + aNodesToInsert.size()); + + // add nodes of face up to first node of link + bool isFLN = false; + nodeIt = theFace->nodesIterator(); + while ( nodeIt->more() && !isFLN ) { + const SMDS_MeshNode* n = static_cast( nodeIt->next() ); + poly_nodes[iNode++] = n; + if (n == nodes[il1]) { + isFLN = true; + } + } + + // add nodes to insert + list::iterator nIt = aNodesToInsert.begin(); + for (; nIt != aNodesToInsert.end(); nIt++) { + poly_nodes[iNode++] = *nIt; + } - // put theNodesToInsert between theBetweenNode1 and theBetweenNode2 - int nbLinkNodes = 2 + theNodesToInsert.size(); + // add nodes of face starting from last node of link + while ( nodeIt->more() ) { + const SMDS_MeshNode* n = static_cast( nodeIt->next() ); + poly_nodes[iNode++] = n; + } + + // edit or replace the face + SMESHDS_Mesh *aMesh = GetMeshDS(); + + if (theFace->IsPoly()) { + aMesh->ChangePolygonNodes(theFace, poly_nodes); + + } else { + int aShapeId = FindShape( theFace ); + + SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes); + if ( aShapeId && newElem ) + aMesh->SetMeshElementOnShape( newElem, aShapeId ); + + aMesh->RemoveElement(theFace); + } + return; + } + + // put aNodesToInsert between theBetweenNode1 and theBetweenNode2 + int nbLinkNodes = 2 + aNodesToInsert.size(); const SMDS_MeshNode* linkNodes[ nbLinkNodes ]; linkNodes[ 0 ] = nodes[ il1 ]; linkNodes[ nbLinkNodes - 1 ] = nodes[ il2 ]; - list::iterator nIt = theNodesToInsert.begin(); - for ( iNode = 1; nIt != theNodesToInsert.end(); nIt++ ) { + list::iterator nIt = aNodesToInsert.begin(); + for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) { linkNodes[ iNode++ ] = *nIt; } // decide how to split a quadrangle: compare possible variants @@ -3053,7 +4467,7 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theFace, // create new elements SMESHDS_Mesh *aMesh = GetMeshDS(); int aShapeId = FindShape( theFace ); - + i1 = 0; i2 = 1; for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ ) { SMDS_MeshElement* newElem = 0; @@ -3079,9 +4493,90 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theFace, aMesh->ChangeElementNodes( theFace, newNodes, iSplit == iBestQuad ? 4 : 3 ); } +//======================================================================= +//function : UpdateVolumes +//purpose : +//======================================================================= +void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode* theBetweenNode1, + const SMDS_MeshNode* theBetweenNode2, + list& theNodesToInsert) +{ + SMDS_ElemIteratorPtr invElemIt = theBetweenNode1->GetInverseElementIterator(); + while (invElemIt->more()) { // loop on inverse elements of theBetweenNode1 + const SMDS_MeshElement* elem = invElemIt->next(); + if (elem->GetType() != SMDSAbs_Volume) + continue; + + // check, if current volume has link theBetweenNode1 - theBetweenNode2 + SMDS_VolumeTool aVolume (elem); + if (!aVolume.IsLinked(theBetweenNode1, theBetweenNode2)) + continue; + + // insert new nodes in all faces of the volume, sharing link theBetweenNode1 - theBetweenNode2 + int iface, nbFaces = aVolume.NbFaces(); + vector poly_nodes; + vector quantities (nbFaces); + + for (iface = 0; iface < nbFaces; iface++) { + int nbFaceNodes = aVolume.NbFaceNodes(iface), nbInserted = 0; + // faceNodes will contain (nbFaceNodes + 1) nodes, last = first + const SMDS_MeshNode** faceNodes = aVolume.GetFaceNodes(iface); + + for (int inode = 0; inode < nbFaceNodes; inode++) { + poly_nodes.push_back(faceNodes[inode]); + + if (nbInserted == 0) { + if (faceNodes[inode] == theBetweenNode1) { + if (faceNodes[inode + 1] == theBetweenNode2) { + nbInserted = theNodesToInsert.size(); + + // add nodes to insert + list::iterator nIt = theNodesToInsert.begin(); + for (; nIt != theNodesToInsert.end(); nIt++) { + poly_nodes.push_back(*nIt); + } + } + } else if (faceNodes[inode] == theBetweenNode2) { + if (faceNodes[inode + 1] == theBetweenNode1) { + nbInserted = theNodesToInsert.size(); + + // add nodes to insert in reversed order + list::iterator nIt = theNodesToInsert.end(); + nIt--; + for (; nIt != theNodesToInsert.begin(); nIt--) { + poly_nodes.push_back(*nIt); + } + poly_nodes.push_back(*nIt); + } + } else { + } + } + } + quantities[iface] = nbFaceNodes + nbInserted; + } + + // Replace or update the volume + SMESHDS_Mesh *aMesh = GetMeshDS(); + + if (elem->IsPoly()) { + aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities); + + } else { + int aShapeId = FindShape( elem ); + + SMDS_MeshElement* newElem = + aMesh->AddPolyhedralVolume(poly_nodes, quantities); + if (aShapeId && newElem) + aMesh->SetMeshElementOnShape(newElem, aShapeId); + + aMesh->RemoveElement(elem); + } + } +} + //======================================================================= //function : SewSideElements -//purpose : +//purpose : //======================================================================= SMESH_MeshEditor::Sew_Error @@ -3191,7 +4686,7 @@ SMESH_MeshEditor::Sew_Error if ( !volSet->empty() ) { //int nodeSetSize = nodeSet->size(); - + // loop on given volumes for ( vIt = volSet->begin(); vIt != volSet->end(); vIt++ ) { SMDS_VolumeTool vol (*vIt); @@ -3210,17 +4705,31 @@ SMESH_MeshEditor::Sew_Error bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second; if ( isNewFace ) { // no such a face is given but it still can exist, check it - if ( nbNodes == 3 ) + if ( nbNodes == 3 ) { aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2] ); - else + } else if ( nbNodes == 4 ) { aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] ); + } else { + vector poly_nodes (nbNodes); + for (int inode = 0; inode < nbNodes; inode++) { + poly_nodes[inode] = fNodes[inode]; + } + aFreeFace = aMesh->FindFace(poly_nodes); + } } if ( !aFreeFace ) { // create a temporary face - if ( nbNodes == 3 ) + if ( nbNodes == 3 ) { aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] ); - else + } else if ( nbNodes == 4 ) { aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] ); + } else { + vector poly_nodes (nbNodes); + for (int inode = 0; inode < nbNodes; inode++) { + poly_nodes[inode] = fNodes[inode]; + } + aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes); + } } if ( aFreeFace ) freeFaceList.push_back( aFreeFace ); @@ -3252,7 +4761,7 @@ SMESH_MeshEditor::Sew_Error maxNbNodes = nbSharedNodes; fIt++; } - else + else freeFaceList.erase( fIt++ ); // here fIt++ occures before erase } if ( freeFaceList.size() > 1 )