From 54d669640def1c7dd6461432564f2c78439fe9ca Mon Sep 17 00:00:00 2001 From: eap Date: Tue, 17 Oct 2017 16:18:05 +0300 Subject: [PATCH] GPUSPHGUI: add FillHole() operation which is needed for NETGEN 2D remesher --- idl/SMESH_MeshEditor.idl | 26 +- src/SMESHUtils/CMakeLists.txt | 5 +- src/SMESHUtils/SMESH_FillHole.cxx | 500 +++++++++++++++++++++++++++ src/SMESHUtils/SMESH_FreeBorders.cxx | 140 +++++++- src/SMESHUtils/SMESH_MeshAlgos.cxx | 265 ++++++++++---- src/SMESHUtils/SMESH_MeshAlgos.hxx | 48 ++- src/SMESH_I/SMESH_2smeshpy.cxx | 2 +- src/SMESH_I/SMESH_MeshEditor_i.cxx | 147 +++++++- src/SMESH_I/SMESH_MeshEditor_i.hxx | 78 +++-- 9 files changed, 1099 insertions(+), 112 deletions(-) create mode 100644 src/SMESHUtils/SMESH_FillHole.cxx diff --git a/idl/SMESH_MeshEditor.idl b/idl/SMESH_MeshEditor.idl index ea19b5ebe..fd5362cda 100644 --- a/idl/SMESH_MeshEditor.idl +++ b/idl/SMESH_MeshEditor.idl @@ -751,7 +751,31 @@ module SMESH * Return point state in a closed 2D mesh in terms of TopAbs_State enumeration. * TopAbs_UNKNOWN state means that either mesh is wrong or the analysis fails. */ - short GetPointState(in double x, in double y, in double z) + short GetPointState(in double x, in double y, in double z) + raises (SALOME::SALOME_Exception); + + /*! + * Check if a 2D mesh is manifold + */ + boolean IsManifold() + raises (SALOME::SALOME_Exception); + + /*! + * Check if orientation of 2D elements is coherent + */ + boolean IsCoherentOrientation2D() + raises (SALOME::SALOME_Exception); + + /*! + * Returns all or only closed FreeBorder's. + */ + ListOfFreeBorders FindFreeBorders(in boolean closedOnly) + raises (SALOME::SALOME_Exception); + + /*! + * Fill with 2D elements a hole defined by a FreeBorder. + */ + void FillHole(in FreeBorder hole) raises (SALOME::SALOME_Exception); /*! diff --git a/src/SMESHUtils/CMakeLists.txt b/src/SMESHUtils/CMakeLists.txt index c84715eed..6750f063c 100644 --- a/src/SMESHUtils/CMakeLists.txt +++ b/src/SMESHUtils/CMakeLists.txt @@ -47,7 +47,7 @@ SET(_link_LIBRARIES ${CAS_TKMesh} ${Boost_LIBRARIES} SMDS -) + ) # --- headers --- @@ -68,7 +68,7 @@ SET(SMESHUtils_HEADERS SMESH_MAT2d.hxx SMESH_ControlPnt.hxx SMESH_Delaunay.hxx -) + ) # --- sources --- @@ -86,6 +86,7 @@ SET(SMESHUtils_SOURCES SMESH_ControlPnt.cxx SMESH_DeMerge.cxx SMESH_Delaunay.cxx + SMESH_FillHole.cxx ) # --- rules --- diff --git a/src/SMESHUtils/SMESH_FillHole.cxx b/src/SMESHUtils/SMESH_FillHole.cxx new file mode 100644 index 000000000..4ded1127f --- /dev/null +++ b/src/SMESHUtils/SMESH_FillHole.cxx @@ -0,0 +1,500 @@ +// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// File : SMESH_FillHole.cxx +// Created : Tue Sep 26 15:11:17 2017 +// Author : Edward AGAPOV (eap) +// + +#include "SMESH_MeshAlgos.hxx" + +#include "SMESH_Comment.hxx" + +#include "SMESH_TypeDefs.hxx" +#include "SMDS_Mesh.hxx" + +#include + +#include +#include + +#include + +namespace +{ + bool isSmallAngle( double cos2 ) + { + // cosine of min angle at which adjacent faces are considered overlapping + const double theMinCos2 = 0.996 * 0.996; // ~5 degrees + return ( cos2 > theMinCos2 ); + } + + struct BEdge; + typedef std::multimap< double, BEdge* > TAngleMap; + typedef std::map< const SMDS_MeshElement*, int > TFaceIndMap; + + //-------------------------------------------------------------------------------- + /*! + * \brief Edge of a free border + */ + struct BEdge + { + const SMDS_MeshNode* myNode1; + const SMDS_MeshNode* myNode2; + const SMDS_MeshElement* myFace; // face adjacent to the border + + gp_XYZ myFaceNorm; + gp_XYZ myDir; // myNode1 -> myNode2 + double myDirCoef; // 1. or -1, to make myDir oriented as myNodes in myFace + double myLength; // between nodes + double myAngleWithPrev; // between myDir and -myPrev->myDir + TAngleMap::iterator myAngleMapPos; + double myOverlapAngle; // angle delta due to overlapping + const SMDS_MeshNode* myNode1Shift; // nodes created to avoid overlapping of faces + const SMDS_MeshNode* myNode2Shift; + + BEdge* myPrev; // neighbors in the border + BEdge* myNext; + + BEdge(): myNode1Shift(0), myNode2Shift(0) {} + void Init( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, + const SMDS_MeshElement* f=0, + const SMDS_MeshNode* nf1=0, const SMDS_MeshNode* nf2=0 ); + void ComputeAngle( bool reverseAngle = false ); + void ShiftOverlapped( const SMDS_MeshNode* oppNode, + const TFaceIndMap& capFaceWithBordInd, + SMDS_Mesh& mesh, + std::vector& newFaces); + void MakeShiftfFaces( SMDS_Mesh& mesh, + std::vector& newFaces, + const bool isReverse ); + gp_XYZ GetInFaceDir() const { return myFaceNorm ^ myDir * myDirCoef; } + void InsertSelf(TAngleMap& edgesByAngle, bool isReverseFaces, bool reBind, bool useOverlap ) + { + if ( reBind ) edgesByAngle.erase( myAngleMapPos ); + double key = (( isReverseFaces ? 2 * M_PI - myAngleWithPrev : myAngleWithPrev ) + + myOverlapAngle * useOverlap ); + myAngleMapPos = edgesByAngle.insert( std::make_pair( key, this )); + } + + // traits used by boost::intrusive::circular_list_algorithms + typedef BEdge node; + typedef BEdge * node_ptr; + typedef const BEdge * const_node_ptr; + static node_ptr get_next(const_node_ptr n) { return n->myNext; } + static void set_next(node_ptr n, node_ptr next) { n->myNext = next; } + static node_ptr get_previous(const_node_ptr n) { return n->myPrev; } + static void set_previous(node_ptr n, node_ptr prev){ n->myPrev = prev; } + }; + + //================================================================================ + /*! + * \brief Initialize a border edge data + */ + //================================================================================ + + void BEdge::Init( const SMDS_MeshNode* n1, + const SMDS_MeshNode* n2, + const SMDS_MeshElement* newFace, // new cap face + const SMDS_MeshNode* nf1, + const SMDS_MeshNode* nf2 ) + { + myNode1 = n1; + myNode2 = n2; + myDir = SMESH_NodeXYZ( n2 ) - SMESH_NodeXYZ( n1 ); + myLength = myDir.Modulus(); + if ( myLength > std::numeric_limits::min() ) + myDir /= myLength; + + myFace = newFace; + if ( !myFace ) + { + TIDSortedElemSet elemSet, avoidSet; + int ind1, ind2; + myFace = SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet, &ind1, &ind2 ); + if ( !myFace ) + throw SALOME_Exception( SMESH_Comment("No face sharing nodes #") + << myNode1->GetID() << " and #" << myNode2->GetID()); + avoidSet.insert( myFace ); + if ( SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet )) + throw SALOME_Exception( SMESH_Comment("No free border between nodes #") + << myNode1->GetID() << " and #" << myNode2->GetID()); + + myDirCoef = SMESH_MeshAlgos::IsRightOrder( myFace, myNode1, myNode2 ) ? 1. : -1.; + } + + if (! SMESH_MeshAlgos::FaceNormal( myFace, myFaceNorm, /*normalized=*/false )) + { + SMDS_ElemIteratorPtr fIt = myNode1->GetInverseElementIterator( SMDSAbs_Face ); + while ( fIt->more() ) + if ( SMESH_MeshAlgos::FaceNormal( fIt->next(), myFaceNorm, /*normalized=*/false )) + break; + } + + if ( newFace ) + { + myFace = 0; + myDirCoef = SMESH_MeshAlgos::IsRightOrder( newFace, nf1, nf2 ) ? 1. : -1.; + if ( myPrev->myNode2 == n1 ) + myNode1Shift = myPrev->myNode2Shift; + if ( myNext->myNode1 == n2 ) + myNode2Shift = myNext->myNode1Shift; + } + else if ( myDirCoef * myPrev->myDirCoef < 0 ) // different orientation of faces + { + myFaceNorm *= -1; + myDirCoef *= -1; + } + + } + + //================================================================================ + /*! + * \brief Compute myAngleWithPrev + */ + //================================================================================ + + void BEdge::ComputeAngle( bool theReverseAngle ) + { + myAngleWithPrev = ACos( myDir.Dot( myPrev->myDir.Reversed() )); + + bool isObtuse; + gp_XYZ inNewFaceDir = myDir - myPrev->myDir; + double dot1 = myDir.Dot( myPrev->myFaceNorm ); + double dot2 = myPrev->myDir.Dot( myFaceNorm ); + bool isOverlap1 = ( inNewFaceDir * myPrev->GetInFaceDir() > 0 ); + bool isOverlap2 = ( inNewFaceDir * GetInFaceDir() > 0 ); + if ( !myPrev->myFace ) + isObtuse = isOverlap1; + else if ( !myFace ) + isObtuse = isOverlap2; + else + { + isObtuse = ( dot1 > 0 || dot2 < 0 ); // suppose face normals point outside the border + if ( theReverseAngle ) + isObtuse = !isObtuse; + } + if ( isObtuse ) + { + myAngleWithPrev = 2 * M_PI - myAngleWithPrev; + } + + // if ( ! isObtuse ) + // isObtuse = + // isSmallAngle( 1 - myDir.CrossSquareMagnitude( myPrev->myDir )); // edges co-directed + + myOverlapAngle = 0; + //if ( !isObtuse ) + { + // check if myFace and a triangle built on this and prev edges overlap + if ( isOverlap1 ) + { + double cos2 = dot1 * dot1 / myPrev->myFaceNorm.SquareModulus(); + myOverlapAngle += 0.5 * M_PI * ( 1 - cos2 ); + } + if ( isOverlap2 ) + { + double cos2 = dot2 * dot2 / myFaceNorm.SquareModulus(); + myOverlapAngle += 0.5 * M_PI * ( 1 - cos2 ); + } + } + } + + //================================================================================ + /*! + * \brief Check if myFace is overlapped by a triangle formed by myNode's and a + * given node. If so, create shifted nodes to avoid overlapping + */ + //================================================================================ + + void BEdge::ShiftOverlapped( const SMDS_MeshNode* theOppNode, + const TFaceIndMap& theCapFaceWithBordInd, + SMDS_Mesh& theMesh, + std::vector& theNewFaces ) + { + if ( myNode1Shift && myNode2Shift ) + return; + + gp_XYZ inNewFaceDir = SMESH_NodeXYZ( theOppNode ) - SMESH_NodeXYZ( myNode1 ); + double dot = inNewFaceDir.Dot( myFaceNorm ); + double cos2 = dot * dot / myFaceNorm.SquareModulus() / inNewFaceDir.SquareModulus(); + bool isOverlap = ( isSmallAngle( 1 - cos2 ) && GetInFaceDir() * inNewFaceDir > 0 ); + + if ( isOverlap ) + { + gp_XYZ shift = myFaceNorm / myLength / 4; + if ( myFace ) + shift.Reverse(); + if ( !myNode1Shift ) + { + gp_XYZ p = SMESH_NodeXYZ( myNode1 ) + shift; + myNode1Shift = theMesh.AddNode( p.X(), p.Y(), p.Z() ); + myPrev->myNode2Shift = myNode1Shift; + } + if ( !myNode2Shift ) + { + gp_XYZ p = SMESH_NodeXYZ( myNode2 ) + shift; + myNode2Shift = theMesh.AddNode( p.X(), p.Y(), p.Z() ); + myNext->myNode1Shift = myNode2Shift; + } + + // MakeShiftfFaces() for already created cap faces + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + const SMDS_MeshNode* ns = is2nd ? myNode2Shift : myNode1Shift; + const SMDS_MeshNode* n = is2nd ? myNode2 : myNode1; + if ( !ns ) continue; + + SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face ); + while ( fIt->more() ) + { + const SMDS_MeshElement* f = fIt->next(); + if ( !f->isMarked() ) continue; + + TFaceIndMap::const_iterator f2i = theCapFaceWithBordInd.find( f ); + if ( f2i == theCapFaceWithBordInd.end() ) + continue; + const SMDS_MeshNode* nf1 = f->GetNode( f2i->second ); + const SMDS_MeshNode* nf2 = f->GetNode(( f2i->second+1 ) % f->NbNodes() ); + if ( nf1 == n || nf2 == n ) + { + BEdge tmpE; + tmpE.myPrev = tmpE.myNext = this; + tmpE.Init( nf1, nf2, f, nf1, nf2 ); + if ( !tmpE.myNode1Shift && !tmpE.myNode2Shift ) + tmpE.Init( nf2, nf1, f, nf2, nf1 ); + tmpE.myFace = f; + tmpE.MakeShiftfFaces( theMesh, theNewFaces, tmpE.myDirCoef < 0 ); + } + std::vector< const SMDS_MeshNode* > nodes( f->begin_nodes(), f->end_nodes() ); + nodes[ f->GetNodeIndex( n ) ] = ns; + theMesh.ChangeElementNodes( f, &nodes[0], nodes.size() ); + } + } + } + } + + //================================================================================ + /*! + * \brief Create a triangle + */ + //================================================================================ + + const SMDS_MeshElement* MakeTria( SMDS_Mesh& mesh, + const SMDS_MeshNode* n1, + const SMDS_MeshNode* n2, + const SMDS_MeshNode* n3, + const bool isReverse ) + { + if ( isReverse ) + return mesh.AddFace( n1, n3, n2 ); + return mesh.AddFace( n1, n2, n3 ); + } + + //================================================================================ + /*! + * \brief Create a quadrangle + */ + //================================================================================ + + // const SMDS_MeshElement* MakeQuad( SMDS_Mesh& mesh, + // const SMDS_MeshNode* n1, + // const SMDS_MeshNode* n2, + // const SMDS_MeshNode* n3, + // const SMDS_MeshNode* n4, + // const bool isReverse ) + // { + // if ( isReverse ) + // return mesh.AddFace( n4, n3, n2, n1 ); + // return mesh.AddFace( n1, n2, n3, n4 ); + // } + + //================================================================================ + /*! + * \brief Create faces on myNode* and myNode*Shift + */ + //================================================================================ + + void BEdge::MakeShiftfFaces(SMDS_Mesh& mesh, + std::vector& newFaces, + const bool isReverse ) + { + if ( !myFace ) + return; + if ( myNode1Shift && myNode2Shift ) + { + newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode2Shift, isReverse )); + newFaces.push_back( MakeTria( mesh, myNode1, myNode2Shift, myNode1Shift, isReverse )); + } + else if ( myNode1Shift ) + { + newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode1Shift, isReverse )); + } + else if ( myNode2Shift ) + { + newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode2Shift, isReverse )); + } + } + +} // namespace + +//================================================================================ +/*! + * \brief Fill with 2D elements a hole defined by a TFreeBorder + */ +//================================================================================ + +void SMESH_MeshAlgos::FillHole(const SMESH_MeshAlgos::TFreeBorder & theFreeBorder, + SMDS_Mesh& theMesh, + std::vector& theNewFaces) +{ + if ( theFreeBorder.size() < 4 || // at least 3 nodes + theFreeBorder[0] != theFreeBorder.back() ) // the hole must be closed + return; + + // prepare data of the border + + ObjectPool< BEdge > edgeAllocator; + boost::intrusive::circular_list_algorithms< BEdge > circularList; + BEdge* edge; + BEdge* edge0 = edgeAllocator.getNew(); + BEdge* edgePrev = edge0; + circularList.init_header( edge0 ); + edge0->Init( theFreeBorder[0], theFreeBorder[1], 0 ); + Bnd_B3d box; + box.Add( SMESH_NodeXYZ( edge0->myNode1 )); + for ( size_t i = 2; i < theFreeBorder.size(); ++i ) + { + edge = edgeAllocator.getNew(); + circularList.link_after( edgePrev, edge ); + edge->Init( theFreeBorder[i-1], theFreeBorder[i] ); + edge->ComputeAngle(); + edgePrev = edge; + box.Add( SMESH_NodeXYZ( edge->myNode1 )); + } + edge0->ComputeAngle(); + + // check if face normals point outside the border + + gp_XYZ hSize = 0.5 * ( box.CornerMax() - box.CornerMin() ); + const double hDelta = 1e-6 * hSize.Modulus(); + hSize -= gp_XYZ( hDelta, hDelta, hDelta ); + if ( hSize.X() < 0 ) hSize.SetX(hDelta); + if ( hSize.Y() < 0 ) hSize.SetY(hDelta); + if ( hSize.Z() < 0 ) hSize.SetZ(hDelta); + box.SetHSize( hSize ); // decrease the box by hDelta + + size_t nbEdges = theFreeBorder.size() - 1; + edge = edge0; + int nbRev = 0, nbFrw = 0; + double angTol = M_PI - ( nbEdges - 2 ) * M_PI / nbEdges, sumDirCoeff = 0; + for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext ) + { + if ( box.IsOut( SMESH_NodeXYZ( edge->myNode1 )) && + edge->myOverlapAngle < 0.1 * M_PI ) + { + nbRev += edge->myAngleWithPrev > M_PI + angTol; + nbFrw += edge->myAngleWithPrev < M_PI - angTol; + } + sumDirCoeff += edge->myDirCoef; + + // unmark all adjacent faces, new faces will be marked + SMDS_ElemIteratorPtr fIt = edge->myNode1->GetInverseElementIterator( SMDSAbs_Face ); + while ( fIt->more() ) + fIt->next()->setIsMarked( false ); + } + bool isReverseAngle = ( nbRev > nbFrw ); // true == face normals point inside the border + //std::cout << "nbRev="<< nbRev << ", nbFrw="<< nbFrw<myNext ) + edge->InsertSelf( edgesByAngle, isReverseAngle, /*reBind=*/false, useOverlap ); + + // create triangles to fill the hole + + //compare order of nodes in the edges with their order in faces + bool isReverse = sumDirCoeff > 0.5 * nbEdges; + + // faces filling the hole (cap faces) and indices of border edges in them + TFaceIndMap capFaceWithBordInd; + + theNewFaces.reserve( nbEdges - 2 ); + std::vector< const SMDS_MeshNode* > nodes(3); + while ( edgesByAngle.size() > 2 ) + { + TAngleMap::iterator a2e = edgesByAngle.begin(); + if ( useOverlap && a2e->first > M_PI - angTol ) // all new triangles need shift + { + // re-sort the edges + useOverlap = false; + edge = a2e->second; + nbEdges = edgesByAngle.size(); + edgesByAngle.clear(); + for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext ) + edge->InsertSelf( edgesByAngle, isReverseAngle, /*reBind=*/false, useOverlap ); + a2e = edgesByAngle.begin(); + } + edge = a2e->second; + edgePrev = edge->myPrev; + + // create shift nodes and faces + edgePrev->ShiftOverlapped( edge->myNode2, capFaceWithBordInd, theMesh, theNewFaces ); + edge->ShiftOverlapped( edgePrev->myNode1, capFaceWithBordInd, theMesh, theNewFaces ); + edge ->MakeShiftfFaces( theMesh, theNewFaces, isReverse ); + edgePrev->MakeShiftfFaces( theMesh, theNewFaces, isReverse ); + + // make a cap face + //nodes.resize( 3 ); + nodes[0] = edgePrev->myNode1Shift ? edgePrev->myNode1Shift : edgePrev->myNode1; + nodes[1] = edgePrev->myNode2Shift ? edgePrev->myNode2Shift : edgePrev->myNode2; + nodes[2] = edge->myNode2Shift ? edge->myNode2Shift : edge->myNode2; + theNewFaces.push_back( MakeTria( theMesh, nodes[0], nodes[1], nodes[2], isReverse )); + // std::cout << nodes[1]->GetID() << " " << nodes[0]->GetID() << " " << nodes[2]->GetID() + // << " " << edge->myAngleWithPrev << std::endl; + + // remember a border edge within the new cap face + theNewFaces.back()->setIsMarked( true ); + if ( edgePrev->myFace ) + capFaceWithBordInd.insert( std::make_pair( theNewFaces.back(), isReverse ? 2 : 0 )); + if ( edge->myFace ) + capFaceWithBordInd.insert( std::make_pair( theNewFaces.back(), 1 )); + + // remove edgePrev from the list and update + edgesByAngle.erase( edgePrev->myAngleMapPos ); + circularList.unlink( edgePrev ); // remove edgePrev from the border + + edge->Init( edgePrev->myNode1, edge->myNode2, theNewFaces.back(), nodes[0], nodes[2] ); + edge->ComputeAngle( isReverseAngle ); + edge->InsertSelf( edgesByAngle, /*isReverse=*/false, /*reBind=*/true, useOverlap ); + edge->myNext->ComputeAngle( isReverseAngle ); + edge->myNext->InsertSelf( edgesByAngle, /*isReverse=*/false, /*reBind=*/true, useOverlap ); + // std::cout << "A " << edge->myNode1->GetID() << " " << edge->myAngleWithPrev + // << " " << edge->myNext->myNode1->GetID() << " " << edge->myNext->myAngleWithPrev + // << std::endl; + } + edge = edgesByAngle.begin()->second; + edge-> MakeShiftfFaces( theMesh, theNewFaces, isReverse ); + edge->myNext->MakeShiftfFaces( theMesh, theNewFaces, isReverse ); +} diff --git a/src/SMESHUtils/SMESH_FreeBorders.cxx b/src/SMESHUtils/SMESH_FreeBorders.cxx index ae6efa6db..d09860465 100644 --- a/src/SMESHUtils/SMESH_FreeBorders.cxx +++ b/src/SMESHUtils/SMESH_FreeBorders.cxx @@ -131,8 +131,11 @@ namespace if ( myID < 0 ) { myID = id; - if ( myNext ) - myNext->SetID( id + 1 ); + + for ( BEdge* be = myNext; be && be->myID < 0; be = be->myNext ) + { + be->myID = ++id; + } } } //================================================================================ @@ -821,3 +824,136 @@ void SMESH_MeshAlgos::FindCoincidentFreeBorders(SMDS_Mesh& mesh, } // SMESH_MeshAlgos::FindCoincidentFreeBorders() +//================================================================================ +/* + * Returns all TFreeBorder's. Optionally check if the mesh is manifold + * and if faces are correctly oriented. + */ +//================================================================================ + +void SMESH_MeshAlgos::FindFreeBorders(SMDS_Mesh& theMesh, + TFreeBorderVec & theFoundFreeBordes, + const bool theClosedOnly, + bool* theIsManifold, + bool* theIsGoodOri) +{ + bool isManifold = true; + + // find free links + typedef NCollection_DataMap TLink2FaceMap; + TLink2FaceMap linkMap; + int nbSharedLinks = 0; + SMDS_FaceIteratorPtr faceIt = theMesh.facesIterator(); + while ( faceIt->more() ) + { + const SMDS_MeshElement* face = faceIt->next(); + if ( !face ) continue; + + const SMDS_MeshNode* n0 = face->GetNode( face->NbNodes() - 1 ); + SMDS_NodeIteratorPtr nodeIt = face->interlacedNodesIterator(); + while ( nodeIt->more() ) + { + const SMDS_MeshNode* n1 = nodeIt->next(); + SMESH_TLink link( n0, n1 ); + if ( const SMDS_MeshElement** faceInMap = linkMap.ChangeSeek( link )) + { + if ( *faceInMap ) + { + if ( theIsGoodOri && *theIsGoodOri && !IsRightOrder( *faceInMap, n1, n0 )) + *theIsGoodOri = false; + } + else + { + isManifold = false; + } + nbSharedLinks += bool( *faceInMap ); + *faceInMap = 0; + } + else + { + linkMap.Bind( link, face ); + } + n0 = n1; + } + } + if ( theIsManifold ) + *theIsManifold = isManifold; + + if ( linkMap.Extent() == nbSharedLinks ) + return; + + // form free borders + std::set < BNode > bNodes; + std::vector< BEdge > bEdges( linkMap.Extent() - nbSharedLinks ); + + TLink2FaceMap::Iterator linkIt( linkMap ); + for ( int iEdge = 0; linkIt.More(); linkIt.Next() ) + { + if ( !linkIt.Value() ) continue; + const SMESH_TLink & link = linkIt.Key(); + std::set< BNode >::iterator n1 = bNodes.insert( BNode( link.node1() )).first; + std::set< BNode >::iterator n2 = bNodes.insert( BNode( link.node2() )).first; + bEdges[ iEdge ].Set( &*n1, &*n2, linkIt.Value(), iEdge+1 ); + n1->AddLinked( & bEdges[ iEdge ] ); + n2->AddLinked( & bEdges[ iEdge ] ); + ++iEdge; + } + linkMap.Clear(); + + // assign IDs to borders + std::vector< BEdge* > borders; // 1st of connected (via myPrev and myNext) edges + std::set< BNode >::iterator bn = bNodes.begin(); + for ( ; bn != bNodes.end(); ++bn ) + { + for ( size_t i = 0; i < bn->myLinkedEdges.size(); ++i ) + { + if ( bn->myLinkedEdges[i]->myBorderID < 0 ) + { + BEdge* be = bn->myLinkedEdges[i]; + int borderID = borders.size(); + borders.push_back( be ); + for ( ; be && be->myBorderID < 0; be = be->myNext ) + { + be->myBorderID = borderID; + be->Orient(); + } + bool isClosed = ( be == bn->myLinkedEdges[i] ); + if ( !isClosed && theClosedOnly ) + { + borders.pop_back(); + continue; + } + be = bn->myLinkedEdges[i]->myPrev; + for ( ; be && be->myBorderID < 0; be = be->myPrev ) + { + be->myBorderID = borderID; + be->Orient(); + } + if ( !isClosed ) + while ( borders.back()->myPrev ) + borders.back() = borders.back()->myPrev; + } + } + } + theFoundFreeBordes.resize( borders.size() ); + for ( size_t i = 0; i < borders.size(); ++i ) + { + TFreeBorder & bordNodes = theFoundFreeBordes[ i ]; + BEdge* be = borders[i]; + + size_t cnt = 1; + for ( be = be->myNext; be && be != borders[i]; be = be->myNext ) + ++cnt; + bordNodes.resize( cnt + 1 ); + + BEdge* beLast; + for ( be = borders[i], cnt = 0; + be && cnt < bordNodes.size()-1; + be = be->myNext, ++cnt ) + { + bordNodes[ cnt ] = be->myBNode1->Node(); + beLast = be; + } + bordNodes.back() = beLast->myBNode2->Node(); + } +} diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index 8464e245f..dbf355e11 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.cxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.cxx @@ -35,6 +35,8 @@ #include "SMDS_VolumeTool.hxx" #include "SMESH_OctreeNode.hxx" +#include + #include #include #include @@ -228,12 +230,12 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(), double tolerance = NodeRadius ); - void prepare(); // !!!call it before calling the following methods!!! void getElementsNearPoint( const gp_Pnt& point, vector& foundElems ); void getElementsNearLine ( const gp_Ax1& line, vector& foundElems); void getElementsInSphere ( const gp_XYZ& center, const double radius, vector& foundElems); + ElementBndBoxTree* getLeafAtPoint( const gp_XYZ& point ); protected: ElementBndBoxTree() {} @@ -337,19 +339,6 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() } } - //================================================================================ - /*! - * \brief Un-mark all elements - */ - //================================================================================ - - void ElementBndBoxTree::prepare() - { - // TElementBoxPool& elBoPool = getElementBoxPool(); - // for ( size_t i = 0; i < elBoPool.nbElements(); ++i ) - // const_cast< ElementBox* >( elBoPool[ i ])->_isMarked = false; - } - //================================================================================ /*! * \brief Return elements which can include the point @@ -471,6 +460,30 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() } } + //================================================================================ + /*! + * \brief Return a leaf including a point + */ + //================================================================================ + + ElementBndBoxTree* ElementBndBoxTree::getLeafAtPoint( const gp_XYZ& point ) + { + if ( getBox()->IsOut( point )) + return 0; + + if ( isLeaf() ) + { + return this; + } + else + { + for (int i = 0; i < 8; i++) + if ( ElementBndBoxTree* l = ((ElementBndBoxTree*) myChildren[i])->getLeafAtPoint( point )) + return l; + } + return 0; + } + //================================================================================ /*! * \brief Construct the element box @@ -539,13 +552,16 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher virtual const SMDS_MeshElement* FindClosestTo( const gp_Pnt& point, SMDSAbs_ElementType type ); - void GetElementsNearLine( const gp_Ax1& line, - SMDSAbs_ElementType type, - vector< const SMDS_MeshElement* >& foundElems); - void GetElementsInSphere( const gp_XYZ& center, - const double radius, - SMDSAbs_ElementType type, - vector< const SMDS_MeshElement* >& foundElems); + virtual void GetElementsNearLine( const gp_Ax1& line, + SMDSAbs_ElementType type, + vector< const SMDS_MeshElement* >& foundElems); + virtual void GetElementsInSphere( const gp_XYZ& center, + const double radius, + SMDSAbs_ElementType type, + vector< const SMDS_MeshElement* >& foundElems); + virtual gp_XYZ Project(const gp_Pnt& point, + SMDSAbs_ElementType type, + const SMDS_MeshElement** closestElem); double getTolerance(); bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face, const double tolerance, double & param); @@ -834,10 +850,6 @@ FindElementsByPoint(const gp_Pnt& point, { _ebbTree[_elementType] = new ElementBndBoxTree( *_mesh, type, _meshPartIt, tolerance ); } - else - { - _ebbTree[ type ]->prepare(); - } vector< const SMDS_MeshElement* > suspectElems; _ebbTree[ type ]->getElementsNearPoint( point, suspectElems ); vector< const SMDS_MeshElement* >::iterator elem = suspectElems.begin(); @@ -863,13 +875,13 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt& point, const SMDS_MeshElement* closestElem = 0; _elementType = type; - if ( type == SMDSAbs_Face || type == SMDSAbs_Volume ) + if ( type == SMDSAbs_Face || + type == SMDSAbs_Volume || + type == SMDSAbs_Edge ) { ElementBndBoxTree*& ebbTree = _ebbTree[ type ]; if ( !ebbTree ) ebbTree = new ElementBndBoxTree( *_mesh, type, _meshPartIt ); - else - ebbTree->prepare(); vector suspectElems; ebbTree->getElementsNearPoint( point, suspectElems ); @@ -885,7 +897,6 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt& point, radius = ebbTree->maxSize() / pow( 2., getTreeHeight()) / 2; while ( suspectElems.empty() ) { - ebbTree->prepare(); ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems ); radius *= 1.1; } @@ -956,8 +967,6 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point) ElementBndBoxTree*& ebbTree = _ebbTree[ SMDSAbs_Face ]; if ( !ebbTree ) ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt ); - else - ebbTree->prepare(); // Algo: analyse transition of a line starting at the point through mesh boundary; // try three lines parallel to axis of the coordinate system and perform rough @@ -974,7 +983,6 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point) gp_Lin line ( lineAxis ); vector suspectFaces; // faces possibly intersecting the line - if ( axis > 0 ) ebbTree->prepare(); ebbTree->getElementsNearLine( lineAxis, suspectFaces ); // Intersect faces with the line @@ -1187,8 +1195,6 @@ void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1& ElementBndBoxTree*& ebbTree = _ebbTree[ type ]; if ( !ebbTree ) ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt ); - else - ebbTree->prepare(); ebbTree->getElementsNearLine( line, foundElems ); } @@ -1208,12 +1214,59 @@ void SMESH_ElementSearcherImpl::GetElementsInSphere( const gp_XYZ& ElementBndBoxTree*& ebbTree = _ebbTree[ type ]; if ( !ebbTree ) ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt ); - else - ebbTree->prepare(); ebbTree->getElementsInSphere( center, radius, foundElems ); } +//======================================================================= +/* + * \brief Return a projection of a given point to a mesh. + * Optionally return the closest element + */ +//======================================================================= + +gp_XYZ SMESH_ElementSearcherImpl::Project(const gp_Pnt& point, + SMDSAbs_ElementType type, + const SMDS_MeshElement** closestElem) +{ + _elementType = type; + if ( _mesh->GetMeshInfo().NbElements( _elementType ) == 0 ) + throw SALOME_Exception( LOCALIZED( "No elements of given type in the mesh" )); + + ElementBndBoxTree*& ebbTree = _ebbTree[ _elementType ]; + if ( !ebbTree ) + ebbTree = new ElementBndBoxTree( *_mesh, _elementType ); + + gp_XYZ p = point.XYZ(); + ElementBndBoxTree* ebbLeaf = ebbTree->getLeafAtPoint( p ); + const Bnd_B3d* box = ebbLeaf->getBox(); + double radius = ( box->CornerMax() - box->CornerMin() ).Modulus(); + + vector< const SMDS_MeshElement* > elems; + ebbTree->getElementsInSphere( p, radius, elems ); + while ( elems.empty() ) + { + radius *= 1.5; + ebbTree->getElementsInSphere( p, radius, elems ); + } + gp_XYZ proj, bestProj; + const SMDS_MeshElement* elem = 0; + double minDist = 2 * radius; + for ( size_t i = 0; i < elems.size(); ++i ) + { + double d = SMESH_MeshAlgos::GetDistance( elems[i], p, &proj ); + if ( d < minDist ) + { + bestProj = proj; + elem = elems[i]; + minDist = d; + } + } + if ( closestElem ) *closestElem = elem; + + return bestProj; +} + //======================================================================= /*! * \brief Return true if the point is IN or ON of the element @@ -1461,17 +1514,19 @@ namespace //======================================================================= double SMESH_MeshAlgos::GetDistance( const SMDS_MeshElement* elem, - const gp_Pnt& point ) + const gp_Pnt& point, + gp_XYZ* closestPnt ) { switch ( elem->GetType() ) { case SMDSAbs_Volume: - return GetDistance( dynamic_cast( elem ), point); + return GetDistance( dynamic_cast( elem ), point, closestPnt ); case SMDSAbs_Face: - return GetDistance( dynamic_cast( elem ), point); + return GetDistance( dynamic_cast( elem ), point, closestPnt ); case SMDSAbs_Edge: - return GetDistance( dynamic_cast( elem ), point); + return GetDistance( dynamic_cast( elem ), point, closestPnt ); case SMDSAbs_Node: + if ( closestPnt ) *closestPnt = SMESH_TNodeXYZ( elem ); return point.Distance( SMESH_TNodeXYZ( elem )); default:; } @@ -1487,9 +1542,10 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshElement* elem, //======================================================================= double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face, - const gp_Pnt& point ) + const gp_Pnt& point, + gp_XYZ* closestPnt ) { - double badDistance = -1; + const double badDistance = -1; if ( !face ) return badDistance; // coordinates of nodes (medium nodes, if any, ignored) @@ -1533,7 +1589,7 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face, trsf.Transforms( tmpPnt ); gp_XY point2D( tmpPnt.X(), tmpPnt.Z() ); - // loop on segments of the face to analyze point position ralative to the face + // loop on edges of the face to analyze point position ralative to the face set< PointPos > pntPosSet; for ( size_t i = 1; i < xy.size(); ++i ) { @@ -1543,31 +1599,40 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face, // compute distance PointPos pos = *pntPosSet.begin(); - // cout << "Face " << face->GetID() << " DIST: "; switch ( pos._name ) { - case POS_LEFT: { - // point is most close to a segment - gp_Vec p0p1( point, xyz[ pos._index ] ); - gp_Vec p1p2( xyz[ pos._index ], xyz[ pos._index+1 ]); // segment vector - p1p2.Normalize(); - double projDist = p0p1 * p1p2; // distance projected to the segment - gp_Vec projVec = p1p2 * projDist; - gp_Vec distVec = p0p1 - projVec; - // cout << distVec.Magnitude() << ", SEG " << face->GetNode(pos._index)->GetID() - // << " - " << face->GetNodeWrap(pos._index+1)->GetID() << endl; - return distVec.Magnitude(); + case POS_LEFT: + { + // point is most close to an edge + gp_Vec edge( xyz[ pos._index ], xyz[ pos._index+1 ]); + gp_Vec n1p ( xyz[ pos._index ], point ); + double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge + // projection of the point on the edge + gp_XYZ proj = ( 1. - u ) * xyz[ pos._index ] + u * xyz[ pos._index+1 ]; + if ( closestPnt ) *closestPnt = proj; + return point.Distance( proj ); } - case POS_RIGHT: { + case POS_RIGHT: + { // point is inside the face - double distToFacePlane = tmpPnt.Y(); - // cout << distToFacePlane << ", INSIDE " << endl; - return Abs( distToFacePlane ); + double distToFacePlane = Abs( tmpPnt.Y() ); + if ( closestPnt ) + { + if ( distToFacePlane < std::numeric_limits::min() ) { + *closestPnt = point.XYZ(); + } + else { + tmpPnt.SetY( 0 ); + trsf.Inverted().Transforms( tmpPnt ); + *closestPnt = tmpPnt; + } + } + return distToFacePlane; } - case POS_VERTEX: { + case POS_VERTEX: + { // point is most close to a node gp_Vec distVec( point, xyz[ pos._index ]); - // cout << distVec.Magnitude() << " VERTEX " << face->GetNode(pos._index)->GetID() << endl; return distVec.Magnitude(); } default:; @@ -1581,7 +1646,9 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face, */ //======================================================================= -double SMESH_MeshAlgos::GetDistance( const SMDS_MeshEdge* seg, const gp_Pnt& point ) +double SMESH_MeshAlgos::GetDistance( const SMDS_MeshEdge* seg, + const gp_Pnt& point, + gp_XYZ* closestPnt ) { double dist = Precision::Infinite(); if ( !seg ) return dist; @@ -1597,16 +1664,25 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshEdge* seg, const gp_Pnt& poi { gp_Vec edge( xyz[i-1], xyz[i] ); gp_Vec n1p ( xyz[i-1], point ); - double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge + double d, u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge if ( u <= 0. ) { - dist = Min( dist, n1p.SquareMagnitude() ); + if (( d = n1p.SquareMagnitude() ) < dist ) { + dist = d; + if ( closestPnt ) *closestPnt = xyz[i-1]; + } } else if ( u >= 1. ) { - dist = Min( dist, point.SquareDistance( xyz[i] )); + if (( d = point.SquareDistance( xyz[i] )) < dist ) { + dist = d; + if ( closestPnt ) *closestPnt = xyz[i]; + } } else { - gp_XYZ proj = ( 1. - u ) * xyz[i-1] + u * xyz[i]; // projection of the point on the edge - dist = Min( dist, point.SquareDistance( proj )); + gp_XYZ proj = xyz[i-1] + u * edge.XYZ(); // projection of the point on the edge + if (( d = point.SquareDistance( proj )) < dist ) { + dist = d; + if ( closestPnt ) *closestPnt = proj; + } } } return Sqrt( dist ); @@ -1620,7 +1696,9 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshEdge* seg, const gp_Pnt& poi */ //======================================================================= -double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt& point ) +double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume, + const gp_Pnt& point, + gp_XYZ* closestPnt ) { SMDS_VolumeTool vTool( volume ); vTool.SetExternalNormal(); @@ -1628,6 +1706,8 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt double n[3], bc[3]; double minDist = 1e100, dist; + gp_XYZ closeP = point.XYZ(); + bool isOut = false; for ( int iF = 0; iF < vTool.NbFaces(); ++iF ) { // skip a facet with normal not "looking at" the point @@ -1644,23 +1724,34 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt case 3: { SMDS_FaceOfNodes tmpFace( nodes[0], nodes[ 1*iQ ], nodes[ 2*iQ ] ); - dist = GetDistance( &tmpFace, point ); + dist = GetDistance( &tmpFace, point, closestPnt ); break; } case 4: { SMDS_FaceOfNodes tmpFace( nodes[0], nodes[ 1*iQ ], nodes[ 2*iQ ], nodes[ 3*iQ ]); - dist = GetDistance( &tmpFace, point ); + dist = GetDistance( &tmpFace, point, closestPnt ); break; } default: vector nvec( nodes, nodes + vTool.NbFaceNodes( iF )); SMDS_PolygonalFaceOfNodes tmpFace( nvec ); - dist = GetDistance( &tmpFace, point ); + dist = GetDistance( &tmpFace, point, closestPnt ); + } + if ( dist < minDist ) + { + minDist = dist; + isOut = true; + if ( closestPnt ) closeP = *closestPnt; } - minDist = Min( minDist, dist ); } - return minDist; + if ( isOut ) + { + if ( closestPnt ) *closestPnt = closeP; + return minDist; + } + + return 0; // point is inside the volume } //================================================================================ @@ -1804,6 +1895,34 @@ vector< const SMDS_MeshNode*> SMESH_MeshAlgos::GetCommonNodes(const SMDS_MeshEle common.push_back( e1->GetNode( i )); return common; } +//================================================================================ +/*! + * \brief Return true if node1 encounters first in the face and node2, after + */ +//================================================================================ + +bool SMESH_MeshAlgos::IsRightOrder( const SMDS_MeshElement* face, + const SMDS_MeshNode* node0, + const SMDS_MeshNode* node1 ) +{ + int i0 = face->GetNodeIndex( node0 ); + int i1 = face->GetNodeIndex( node1 ); + if ( face->IsQuadratic() ) + { + if ( face->IsMediumNode( node0 )) + { + i0 -= ( face->NbNodes()/2 - 1 ); + i1 *= 2; + } + else + { + i1 -= ( face->NbNodes()/2 - 1 ); + i0 *= 2; + } + } + int diff = i1 - i0; + return ( diff == 1 ) || ( diff == -face->NbNodes()+1 ); +} //======================================================================= /*! diff --git a/src/SMESHUtils/SMESH_MeshAlgos.hxx b/src/SMESHUtils/SMESH_MeshAlgos.hxx index 1114bfe02..10af7a636 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.hxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.hxx @@ -100,6 +100,15 @@ struct SMESHUtils_EXPORT SMESH_ElementSearcher * \brief Find out if the given point is out of closed 2D mesh. */ virtual TopAbs_State GetPointState(const gp_Pnt& point) = 0; + + /*! + * \brief Return a projection of a given point to a 2D mesh. + * Optionally return the closest face + */ + virtual gp_XYZ Project(const gp_Pnt& point, + SMDSAbs_ElementType type, + const SMDS_MeshElement** closestFace= 0) = 0; + virtual ~SMESH_ElementSearcher(); }; @@ -112,16 +121,16 @@ namespace SMESH_MeshAlgos bool IsOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol ); SMESHUtils_EXPORT - double GetDistance( const SMDS_MeshElement* elem, const gp_Pnt& point ); + double GetDistance( const SMDS_MeshElement* elem, const gp_Pnt& point, gp_XYZ* closestPnt = 0 ); SMESHUtils_EXPORT - double GetDistance( const SMDS_MeshEdge* edge, const gp_Pnt& point ); + double GetDistance( const SMDS_MeshEdge* edge, const gp_Pnt& point, gp_XYZ* closestPnt = 0 ); SMESHUtils_EXPORT - double GetDistance( const SMDS_MeshFace* face, const gp_Pnt& point ); + double GetDistance( const SMDS_MeshFace* face, const gp_Pnt& point, gp_XYZ* closestPnt = 0 ); SMESHUtils_EXPORT - double GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt& point ); + double GetDistance( const SMDS_MeshVolume* volume, const gp_Pnt& point, gp_XYZ* closestPnt = 0 ); SMESHUtils_EXPORT void GetBarycentricCoords( const gp_XY& point, @@ -153,6 +162,14 @@ namespace SMESH_MeshAlgos SMESHUtils_EXPORT std::vector< const SMDS_MeshNode*> GetCommonNodes(const SMDS_MeshElement* e1, const SMDS_MeshElement* e2); + /*! + * \brief Return true if node1 encounters first in the face and node2, after. + * The nodes are supposed to be neighbor nodes in the face. + */ + SMESHUtils_EXPORT + bool IsRightOrder( const SMDS_MeshElement* face, + const SMDS_MeshNode* node0, + const SMDS_MeshNode* node1 ); /*! * \brief Return SMESH_NodeSearcher. The caller is responsible for deleteing it @@ -204,7 +221,28 @@ namespace SMESH_MeshAlgos void FindCoincidentFreeBorders(SMDS_Mesh& mesh, double tolerance, CoincidentFreeBorders & foundFreeBordes); - + /*! + * Returns all or only closed TFreeBorder's. + * Optionally check if the mesh is manifold and if faces are correctly oriented. + * + * (Implemented in ./SMESH_FreeBorders.cxx) + */ + SMESHUtils_EXPORT + void FindFreeBorders(SMDS_Mesh& mesh, + TFreeBorderVec & foundFreeBordes, + const bool closedOnly, + bool* isManifold = 0, + bool* isGoodOri = 0); + /*! + * Fill a hole defined by a TFreeBorder with 2D elements. + * + * (Implemented in ./SMESH_FillHole.cxx) + */ + SMESHUtils_EXPORT + void FillHole(const TFreeBorder & freeBorder, + SMDS_Mesh& mesh, + std::vector& newFaces); + /*! * \brief Find nodes whose merge makes the element invalid diff --git a/src/SMESH_I/SMESH_2smeshpy.cxx b/src/SMESH_I/SMESH_2smeshpy.cxx index 5a9e2e519..38fc1414f 100644 --- a/src/SMESH_I/SMESH_2smeshpy.cxx +++ b/src/SMESH_I/SMESH_2smeshpy.cxx @@ -2430,7 +2430,7 @@ void _pyMeshEditor::Process( const Handle(_pyCommand)& theCommand) "ExtrusionAlongPathX","ExtrusionAlongPathObject1D","ExtrusionAlongPathObject2D", "ExtrusionSweepObjects","RotationSweepObjects","ExtrusionAlongPathObjects", "Mirror","MirrorObject","Translate","TranslateObject","Rotate","RotateObject", - "FindCoincidentNodes","MergeNodes","FindEqualElements", + "FindCoincidentNodes","MergeNodes","FindEqualElements","FillHole", "MergeElements","MergeEqualElements","SewFreeBorders","SewConformFreeBorders", "FindCoincidentFreeBorders", "SewCoincidentFreeBorders", "SewBorderToSide","SewSideElements","ChangeElemNodes","GetLastCreatedNodes", diff --git a/src/SMESH_I/SMESH_MeshEditor_i.cxx b/src/SMESH_I/SMESH_MeshEditor_i.cxx index 6dfd0d61e..bca571b9a 100644 --- a/src/SMESH_I/SMESH_MeshEditor_i.cxx +++ b/src/SMESH_I/SMESH_MeshEditor_i.cxx @@ -4006,7 +4006,7 @@ SMESH_MeshEditor_i::ScaleMakeMesh(SMESH::SMESH_IDSource_ptr theObject, // and then "GetGroups" using SMESH_Mesh::GetGroups() TPythonDump pydump; // to prevent dump at mesh creation - mesh = makeMesh( theMeshName ); + mesh = makeMesh( theMeshName ); mesh_i = SMESH::DownCast( mesh ); if ( mesh_i ) @@ -4593,6 +4593,151 @@ CORBA::Short SMESH_MeshEditor_i::GetPointState(CORBA::Double x, return 0; } +//======================================================================= +//function : IsManifold +//purpose : Check if a 2D mesh is manifold +//======================================================================= + +CORBA::Boolean SMESH_MeshEditor_i::IsManifold() + throw (SALOME::SALOME_Exception) +{ + bool isManifold = true; + + SMESH_TRY; + SMESH_MeshAlgos::TFreeBorderVec foundFreeBordes; + SMESH_MeshAlgos::FindFreeBorders( *getMeshDS(), + foundFreeBordes, + /*closedOnly=*/true, + &isManifold ); + SMESH_CATCH( SMESH::throwCorbaException ); + + return isManifold; +} + +//======================================================================= +//function : IsCoherentOrientation2D +//purpose : Check if orientation of 2D elements is coherent +//======================================================================= + +CORBA::Boolean SMESH_MeshEditor_i::IsCoherentOrientation2D() + throw (SALOME::SALOME_Exception) +{ + bool isGoodOri = true; + + SMESH_TRY; + SMESH_MeshAlgos::TFreeBorderVec foundFreeBordes; + SMESH_MeshAlgos::FindFreeBorders( *getMeshDS(), + foundFreeBordes, + /*closedOnly=*/true, + /*isManifold=*/0, + &isGoodOri); + SMESH_CATCH( SMESH::throwCorbaException ); + + return isGoodOri; +} + +//======================================================================= +//function : FindFreeBorders +//purpose : Returns all or only closed FreeBorder's. +//======================================================================= + +SMESH::ListOfFreeBorders* SMESH_MeshEditor_i::FindFreeBorders(CORBA::Boolean closedOnly) + throw (SALOME::SALOME_Exception) +{ + SMESH::ListOfFreeBorders_var resBorders = new SMESH::ListOfFreeBorders; + SMESH_TRY; + + SMESH_MeshAlgos::TFreeBorderVec foundFreeBordes; + SMESH_MeshAlgos::FindFreeBorders( *getMeshDS(), foundFreeBordes, closedOnly ); + + resBorders->length( foundFreeBordes.size() ); + for ( size_t i = 0; i < foundFreeBordes.size(); ++i ) + { + const SMESH_MeshAlgos::TFreeBorder& bordNodes = foundFreeBordes[i]; + SMESH::FreeBorder& bordOut = resBorders[i]; + bordOut.nodeIDs.length( bordNodes.size() ); + for ( size_t iN = 0; iN < bordNodes.size(); ++iN ) + bordOut.nodeIDs[ iN ] = bordNodes[ iN ]->GetID(); + } + + SMESH_CATCH( SMESH::throwCorbaException ); + + return resBorders._retn(); +} + +//======================================================================= +//function : FillHole +//purpose : Fill with 2D elements a hole defined by a FreeBorder. +//======================================================================= + +void SMESH_MeshEditor_i::FillHole(const SMESH::FreeBorder& theHole) + throw (SALOME::SALOME_Exception) +{ + initData(); + + if ( theHole.nodeIDs.length() < 4 ) + THROW_SALOME_CORBA_EXCEPTION("A hole should be bound by at least 3 nodes", SALOME::BAD_PARAM); + if ( theHole.nodeIDs[0] != theHole.nodeIDs[ theHole.nodeIDs.length()-1 ] ) + THROW_SALOME_CORBA_EXCEPTION("Not closed hole boundary. " + "First and last nodes must be same", SALOME::BAD_PARAM); + + SMESH_MeshAlgos::TFreeBorder bordNodes; + bordNodes.resize( theHole.nodeIDs.length() ); + for ( size_t iN = 0; iN < theHole.nodeIDs.length(); ++iN ) + { + bordNodes[ iN ] = getMeshDS()->FindNode( theHole.nodeIDs[ iN ]); + if ( !bordNodes[ iN ] ) + THROW_SALOME_CORBA_EXCEPTION(SMESH_Comment("Node #") << theHole.nodeIDs[ iN ] + << " does not exist", SALOME::BAD_PARAM); + } + + SMESH_TRY; + + MeshEditor_I::TPreviewMesh* previewMesh = 0; + SMDS_Mesh* meshDS = getMeshDS(); + if ( myIsPreviewMode ) + { + // copy faces sharing nodes of theHole + TIDSortedElemSet holeFaces; + previewMesh = getPreviewMesh( SMDSAbs_Face ); + for ( size_t i = 0; i < bordNodes.size(); ++i ) + { + SMDS_ElemIteratorPtr fIt = bordNodes[i]->GetInverseElementIterator( SMDSAbs_Face ); + while ( fIt->more() ) + { + const SMDS_MeshElement* face = fIt->next(); + if ( holeFaces.insert( face ).second ) + previewMesh->Copy( face ); + } + bordNodes[i] = previewMesh->GetMeshDS()->FindNode( bordNodes[i]->GetID() ); + ASSERT( bordNodes[i] ); + } + meshDS = previewMesh->GetMeshDS(); + } + + std::vector newFaces; + SMESH_MeshAlgos::FillHole( bordNodes, *meshDS, newFaces ); + + if ( myIsPreviewMode ) + { + previewMesh->Clear(); + for ( size_t i = 0; i < newFaces.size(); ++i ) + previewMesh->Copy( newFaces[i] ); + } + else + { + getEditor().ClearLastCreated(); + SMESH_SequenceOfElemPtr& aSeq = + const_cast( getEditor().GetLastCreatedElems() ); + for ( size_t i = 0; i < newFaces.size(); ++i ) + aSeq.Append( newFaces[i] ); + + TPythonDump() << this << ".FillHole( SMESH.FreeBorder(" << theHole.nodeIDs << " ))"; + } + + SMESH_CATCH( SMESH::throwCorbaException ); +} + //======================================================================= //function : convError //purpose : diff --git a/src/SMESH_I/SMESH_MeshEditor_i.hxx b/src/SMESH_I/SMESH_MeshEditor_i.hxx index 74eaa36ad..5eaa5ee79 100644 --- a/src/SMESH_I/SMESH_MeshEditor_i.hxx +++ b/src/SMESH_I/SMESH_MeshEditor_i.hxx @@ -527,7 +527,7 @@ public: SMESH::ElementType type) throw (SALOME::SALOME_Exception); /*! - * Searching among the given elements, return elements of given type + * Searching among the given elements, return elements of given type * where the given point is IN or ON. * 'ALL' type means elements of any type excluding nodes */ @@ -545,6 +545,30 @@ public: CORBA::Short GetPointState(CORBA::Double x, CORBA::Double y, CORBA::Double z) throw (SALOME::SALOME_Exception); + /*! + * Check if a 2D mesh is manifold + */ + CORBA::Boolean IsManifold() + throw (SALOME::SALOME_Exception); + + /*! + * Check if orientation of 2D elements is coherent + */ + CORBA::Boolean IsCoherentOrientation2D() + throw (SALOME::SALOME_Exception); + + /*! + * Returns all or only closed FreeBorder's. + */ + SMESH::ListOfFreeBorders* FindFreeBorders(CORBA::Boolean closedOnly) + throw (SALOME::SALOME_Exception); + + /*! + * Fill with 2D elements a hole defined by a FreeBorder. + */ + void FillHole(const SMESH::FreeBorder& hole) + throw (SALOME::SALOME_Exception); + SMESH::CoincidentFreeBorders* FindCoincidentFreeBorders(CORBA::Double tolerance); CORBA::Short SewCoincidentFreeBorders(const SMESH::CoincidentFreeBorders& freeBorders, CORBA::Boolean createPolygons, @@ -552,35 +576,35 @@ public: throw (SALOME::SALOME_Exception); SMESH::SMESH_MeshEditor::Sew_Error - SewFreeBorders(CORBA::Long FirstNodeID1, - CORBA::Long SecondNodeID1, - CORBA::Long LastNodeID1, - CORBA::Long FirstNodeID2, - CORBA::Long SecondNodeID2, - CORBA::Long LastNodeID2, - CORBA::Boolean CreatePolygons, - CORBA::Boolean CreatePolyedrs) throw (SALOME::SALOME_Exception); + SewFreeBorders(CORBA::Long FirstNodeID1, + CORBA::Long SecondNodeID1, + CORBA::Long LastNodeID1, + CORBA::Long FirstNodeID2, + CORBA::Long SecondNodeID2, + CORBA::Long LastNodeID2, + CORBA::Boolean CreatePolygons, + CORBA::Boolean CreatePolyedrs) throw (SALOME::SALOME_Exception); SMESH::SMESH_MeshEditor::Sew_Error - SewConformFreeBorders(CORBA::Long FirstNodeID1, - CORBA::Long SecondNodeID1, - CORBA::Long LastNodeID1, - CORBA::Long FirstNodeID2, - CORBA::Long SecondNodeID2) throw (SALOME::SALOME_Exception); + SewConformFreeBorders(CORBA::Long FirstNodeID1, + CORBA::Long SecondNodeID1, + CORBA::Long LastNodeID1, + CORBA::Long FirstNodeID2, + CORBA::Long SecondNodeID2) throw (SALOME::SALOME_Exception); SMESH::SMESH_MeshEditor::Sew_Error - SewBorderToSide(CORBA::Long FirstNodeIDOnFreeBorder, - CORBA::Long SecondNodeIDOnFreeBorder, - CORBA::Long LastNodeIDOnFreeBorder, - CORBA::Long FirstNodeIDOnSide, - CORBA::Long LastNodeIDOnSide, - CORBA::Boolean CreatePolygons, - CORBA::Boolean CreatePolyedrs) throw (SALOME::SALOME_Exception); + SewBorderToSide(CORBA::Long FirstNodeIDOnFreeBorder, + CORBA::Long SecondNodeIDOnFreeBorder, + CORBA::Long LastNodeIDOnFreeBorder, + CORBA::Long FirstNodeIDOnSide, + CORBA::Long LastNodeIDOnSide, + CORBA::Boolean CreatePolygons, + CORBA::Boolean CreatePolyedrs) throw (SALOME::SALOME_Exception); SMESH::SMESH_MeshEditor::Sew_Error - SewSideElements(const SMESH::long_array& IDsOfSide1Elements, - const SMESH::long_array& IDsOfSide2Elements, - CORBA::Long NodeID1OfSide1ToMerge, - CORBA::Long NodeID1OfSide2ToMerge, - CORBA::Long NodeID2OfSide1ToMerge, - CORBA::Long NodeID2OfSide2ToMerge) throw (SALOME::SALOME_Exception); + SewSideElements(const SMESH::long_array& IDsOfSide1Elements, + const SMESH::long_array& IDsOfSide2Elements, + CORBA::Long NodeID1OfSide1ToMerge, + CORBA::Long NodeID1OfSide2ToMerge, + CORBA::Long NodeID2OfSide1ToMerge, + CORBA::Long NodeID2OfSide2ToMerge) throw (SALOME::SALOME_Exception); /*! * Set new nodes for given element.