From: eap Date: Thu, 20 Apr 2017 11:16:23 +0000 (+0300) Subject: 54122: Bad quality prismatic mesh X-Git-Tag: Before_multi_study_removal_06072017~8 X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=commitdiff_plain;h=f6b5d2f920970263bd4aa50e6ef7116d5c1b1625 54122: Bad quality prismatic mesh Extract class Delaunay from class Morph + Print warning when creating an algorithm with unexpected arguments (expected ones are algo type and geometry) -- smeshBuilder.py --- diff --git a/src/SMESHUtils/CMakeLists.txt b/src/SMESHUtils/CMakeLists.txt index 3641d913c..c84715eed 100644 --- a/src/SMESHUtils/CMakeLists.txt +++ b/src/SMESHUtils/CMakeLists.txt @@ -67,6 +67,7 @@ SET(SMESHUtils_HEADERS SMESH_MeshAlgos.hxx SMESH_MAT2d.hxx SMESH_ControlPnt.hxx + SMESH_Delaunay.hxx ) # --- sources --- @@ -84,6 +85,7 @@ SET(SMESHUtils_SOURCES SMESH_FreeBorders.cxx SMESH_ControlPnt.cxx SMESH_DeMerge.cxx + SMESH_Delaunay.cxx ) # --- rules --- diff --git a/src/SMESHUtils/SMESH_Delaunay.cxx b/src/SMESHUtils/SMESH_Delaunay.cxx new file mode 100644 index 000000000..d5bc2f663 --- /dev/null +++ b/src/SMESHUtils/SMESH_Delaunay.cxx @@ -0,0 +1,296 @@ +// 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_Delaunay.cxx +// Created : Wed Apr 19 15:41:15 2017 +// Author : Edward AGAPOV (eap) + +#include "SMESH_Delaunay.hxx" +#include "SMESH_MeshAlgos.hxx" + +#include +#include + +//================================================================================ +/*! + * \brief Construct a Delaunay triangulation of given boundary nodes + * \param [in] boundaryNodes - vector of nodes of a wire + * \param [in] face - the face + * \param [in] faceID - the face ID + * \param [in] nbNodesToVisit - nb of non-marked nodes on the face + */ +//================================================================================ + +SMESH_Delaunay::SMESH_Delaunay(const std::vector< const UVPtStructVec* > & boundaryNodes, + const TopoDS_Face& face, + const int faceID) + : _face( face ), _faceID( faceID ), _scale( 1., 1. ) +{ + // compute _scale + BRepAdaptor_Surface surf( face ); + if ( surf.GetType() != GeomAbs_Plane ) + { + const int nbDiv = 100; + const double uRange = surf.LastUParameter() - surf.FirstUParameter(); + const double vRange = surf.LastVParameter() - surf.FirstVParameter(); + const double dU = uRange / nbDiv; + const double dV = vRange / nbDiv; + double u = surf.FirstUParameter(), v = surf.FirstVParameter(); + gp_Pnt p0U = surf.Value( u, v ), p0V = p0U; + double lenU = 0, lenV = 0; + for ( ; u < surf.LastUParameter(); u += dU, v += dV ) + { + gp_Pnt p1U = surf.Value( u, surf.FirstVParameter() ); + lenU += p1U.Distance( p0U ); + p0U = p1U; + gp_Pnt p1V = surf.Value( surf.FirstUParameter(), v ); + lenV += p1V.Distance( p0V ); + p0V = p1V; + } + _scale.SetCoord( lenU / uRange, lenV / vRange ); + } + + // count boundary points + int iP = 1, nbP = 0; + for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW ) // loop on wires + { + nbP += boundaryNodes[iW]->size(); + if ( boundaryNodes[iW]->front().node == boundaryNodes[iW]->back().node ) + --nbP; // 1st and last points coincide + } + _bndNodes.resize( nbP ); + + // fill boundary points + BRepMesh::Array1OfVertexOfDelaun bndVert( 1, 1 + nbP ); + BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier ); + for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW ) + { + const UVPtStructVec& bndPnt = *boundaryNodes[iW]; + int i = 0, nb = bndPnt.size(); + if ( bndPnt[0].node == bndPnt.back().node ) + --nb; + for ( ; i < nb; ++i, ++iP ) + { + _bndNodes[ iP-1 ] = bndPnt[i].node; + bndPnt[i].node->setIsMarked( true ); + + v.ChangeCoord() = bndPnt[i].UV().Multiplied( _scale ); + bndVert( iP ) = v; + } + } + + // triangulate the srcFace in 2D + BRepMesh_Delaun Delaunay( bndVert ); + _triaDS = Delaunay.Result(); +} + +//================================================================================ +/*! + * \brief Prepare to the exploration of nodes + */ +//================================================================================ + +void SMESH_Delaunay::InitTraversal(const int nbNodesToVisit) +{ + _nbNodesToVisit = (size_t) nbNodesToVisit; + _nbVisitedNodes = _iBndNode = 0; + _noTriQueue.clear(); +} + +//================================================================================ +/*! + * \brief Return a node with its Barycentric Coordinates within the triangle + * defined by its node indices (zero based) + * \param [out] bc - Barycentric Coordinates of the returned node + * \param [out] triaNodes - indices of triangle nodes + * \return const SMDS_MeshNode* - the next node or NULL + */ +//================================================================================ + +const SMDS_MeshNode* SMESH_Delaunay::NextNode( double bc[3], int triaNodes[3] ) +{ + while ( _nbVisitedNodes < _nbNodesToVisit ) + { + while ( !_noTriQueue.empty() ) + { + const SMDS_MeshNode* node = _noTriQueue.front().first; + const BRepMesh_Triangle* tria = _noTriQueue.front().second; + _noTriQueue.pop_front(); + if ( node->isMarked() ) + continue; + ++_nbVisitedNodes; + node->setIsMarked( true ); + + // find a Delaunay triangle containing the src node + gp_XY uv = getNodeUV( _face, node ); + tria = FindTriangle( uv, tria, bc, triaNodes ); + if ( tria ) + { + addCloseNodes( node, tria, _faceID, _noTriQueue ); + return node; + } + } + for ( ; _iBndNode < _bndNodes.size() && _noTriQueue.empty(); ++_iBndNode ) + { + if ( const BRepMesh_Triangle* tria = GetTriangleNear( _iBndNode )) + addCloseNodes( _bndNodes[ _iBndNode ], tria, _faceID, _noTriQueue ); + } + if ( _noTriQueue.empty() ) + break; + } + + // if ( _nbVisitedNodes < _nbNodesToVisit ) + // _nbVisitedNodes = std::numeric_limits::max(); + return NULL; +} + +//================================================================================ +/*! + * \brief Find a Delaunay triangle containing a given 2D point and return + * barycentric coordinates within the found triangle + */ +//================================================================================ + +const BRepMesh_Triangle* SMESH_Delaunay::FindTriangle( const gp_XY& UV, + const BRepMesh_Triangle* tria, + double bc[3], + int triaNodes[3] ) +{ + int nodeIDs[3]; + gp_XY nodeUVs[3]; + int linkIDs[3]; + Standard_Boolean ori[3]; + + gp_XY uv = UV.Multiplied( _scale ); + + while ( tria ) + { + // check if the uv is in tria + + _triaDS->ElementNodes( *tria, nodeIDs ); + nodeUVs[0] = _triaDS->GetNode( nodeIDs[0] ).Coord(); + nodeUVs[1] = _triaDS->GetNode( nodeIDs[1] ).Coord(); + nodeUVs[2] = _triaDS->GetNode( nodeIDs[2] ).Coord(); + + SMESH_MeshAlgos::GetBarycentricCoords( uv, + nodeUVs[0], nodeUVs[1], nodeUVs[2], + bc[0], bc[1] ); + if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 ) + { + bc[2] = 1 - bc[0] - bc[1]; + triaNodes[0] = nodeIDs[0] - 1; + triaNodes[1] = nodeIDs[1] - 1; + triaNodes[2] = nodeIDs[2] - 1; + return tria; + } + + // look for a neighbor triangle, which is adjacent to a link intersected + // by a segment( triangle center -> uv ) + + gp_XY gc = ( nodeUVs[0] + nodeUVs[1] + nodeUVs[2] ) / 3.; + gp_XY seg = uv - gc; + + tria->Edges( linkIDs, ori ); + int triaID = _triaDS->IndexOf( *tria ); + tria = 0; + + for ( int i = 0; i < 3; ++i ) + { + const BRepMesh_PairOfIndex & triIDs = _triaDS->ElementsConnectedTo( linkIDs[i] ); + if ( triIDs.Extent() < 2 ) + continue; // no neighbor triangle + + // check if a link intersects gc2uv + const BRepMesh_Edge & link = _triaDS->GetLink( linkIDs[i] ); + const BRepMesh_Vertex & n1 = _triaDS->GetNode( link.FirstNode() ); + const BRepMesh_Vertex & n2 = _triaDS->GetNode( link.LastNode() ); + gp_XY uv1 = n1.Coord(); + gp_XY lin = n2.Coord() - uv1; // link direction + + double crossSegLin = seg ^ lin; + if ( Abs( crossSegLin ) < std::numeric_limits::min() ) + continue; // parallel + + double uSeg = ( uv1 - gc ) ^ lin / crossSegLin; + if ( 0. <= uSeg && uSeg <= 1. ) + { + tria = & _triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID ))); + break; + } + } + } + return tria; +} + +//================================================================================ +/*! + * \brief Return a triangle sharing a given boundary node + * \param [in] iBndNode - index of the boundary node + * \return const BRepMesh_Triangle* - a found triangle + */ +//================================================================================ + +const BRepMesh_Triangle* SMESH_Delaunay::GetTriangleNear( int iBndNode ) +{ + const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndNode + 1 ); + const BRepMesh_PairOfIndex & triaIds = _triaDS->ElementsConnectedTo( linkIds.First() ); + const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(1) ); + return &tria; +} + +//================================================================================ +/*! + * \brief Return UV of the i-th source boundary node (zero based) + */ +//================================================================================ + +gp_XY SMESH_Delaunay::GetBndUV(const int iNode) const +{ + return _triaDS->GetNode( iNode+1 ).Coord(); +} + +//================================================================================ +/*! + * \brief Add non-marked nodes surrounding a given one to a queue + */ +//================================================================================ + +void SMESH_Delaunay::addCloseNodes( const SMDS_MeshNode* node, + const BRepMesh_Triangle* tria, + const int faceID, + TNodeTriaList & _noTriQueue ) +{ + // find in-FACE nodes + SMDS_ElemIteratorPtr elems = node->GetInverseElementIterator(SMDSAbs_Face); + while ( elems->more() ) + { + const SMDS_MeshElement* elem = elems->next(); + if ( elem->getshapeId() == faceID ) + { + for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i ) + { + const SMDS_MeshNode* n = elem->GetNode( i ); + if ( !n->isMarked() /*&& n->getshapeId() == faceID*/ ) + _noTriQueue.push_back( std::make_pair( n, tria )); + } + } + } +} diff --git a/src/SMESHUtils/SMESH_Delaunay.hxx b/src/SMESHUtils/SMESH_Delaunay.hxx new file mode 100644 index 000000000..8d04f506a --- /dev/null +++ b/src/SMESHUtils/SMESH_Delaunay.hxx @@ -0,0 +1,112 @@ +// 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_Delaunay.hxx +// Created : Wed Apr 19 15:42:54 2017 +// Author : Edward AGAPOV (eap) + + +#ifndef __SMESH_Delaunay_HXX__ +#define __SMESH_Delaunay_HXX__ + +#include "SMESH_TypeDefs.hxx" + +#include + +/*! + * \brief Create a Delaunay triangulation of nodes on a face boundary + * and provide exploration of nodes shared by elements lying on + * the face. For a returned node, also return a Delaunay triangle + * the node lies in and its Barycentric Coordinates within the triangle. + * Only non-marked nodes are visited. Boundary nodes given at the construction + * are not returned. + * + * For usage, this class needs to be subclassed to implement getNodeUV(); + */ +class SMESHUtils_EXPORT SMESH_Delaunay +{ + public: + + // construct a Delaunay triangulation of given boundary nodes + SMESH_Delaunay(const std::vector< const UVPtStructVec* > & boundaryNodes, + const TopoDS_Face& face, + const int faceID); + + virtual ~SMESH_Delaunay() {} + + // prepare to the exploration of nodes + void InitTraversal(const int nbNodesToVisit = -1); + + // return a node with its Barycentric Coordinates within the triangle + // defined by its node indices (zero based) + const SMDS_MeshNode* NextNode( double bc[3], int triaNodes[3] ); + + // return nb of nodes returned by NextNode() + int NbVisitedNodes() const { return _nbVisitedNodes; } + + + // find a triangle containing an UV, starting from a given triangle; + // return barycentric coordinates of the UV and the found triangle (indices are zero based). + const BRepMesh_Triangle* FindTriangle( const gp_XY& uv, + const BRepMesh_Triangle* bmTria, + double bc[3], + int triaNodes[3]); + + // return any Delaunay triangle neighboring a given boundary node (zero based) + const BRepMesh_Triangle* GetTriangleNear( int iBndNode ); + + // return source boundary nodes + const std::vector< const SMDS_MeshNode* >& GetBndNodes() const { return _bndNodes; } + + // return UV of the i-th source boundary node (zero based) + gp_XY GetBndUV(const int iNode) const; + + // return scale factor to convert real UV to/from UV used for Delauney meshing: + // delauney_UV = real_UV * scale + const gp_XY& GetScale() const { return _scale; } + + + protected: + + // container of a node and a triangle serving as a start while searching a + // triangle including the node UV + typedef std::list< std::pair< const SMDS_MeshNode*, const BRepMesh_Triangle* > > TNodeTriaList; + + // return UV of a node on the face + virtual gp_XY getNodeUV( const TopoDS_Face& face, const SMDS_MeshNode* node ) const = 0; + + // add non-marked nodes surrounding a given one to a queue + static void addCloseNodes( const SMDS_MeshNode* node, + const BRepMesh_Triangle* bmTria, + const int faceID, + TNodeTriaList & noTriQueue ); + + const TopoDS_Face& _face; + int _faceID; + std::vector< const SMDS_MeshNode* > _bndNodes; + gp_XY _scale; + Handle(BRepMesh_DataStructureOfDelaun) _triaDS; + size_t _nbNodesToVisit, _nbVisitedNodes, _iBndNode; + TNodeTriaList _noTriQueue; + +}; + +#endif diff --git a/src/SMESH_SWIG/smeshBuilder.py b/src/SMESH_SWIG/smeshBuilder.py index 900aecf18..88d2b0b1a 100644 --- a/src/SMESH_SWIG/smeshBuilder.py +++ b/src/SMESH_SWIG/smeshBuilder.py @@ -5182,10 +5182,11 @@ omniORB.registerObjref(SMESH._objref_SMESH_Pattern._NP_RepositoryId, Pattern) ## Private class used to bind methods creating algorithms to the class Mesh # class algoCreator: - def __init__(self): + def __init__(self, method): self.mesh = None self.defaultAlgoType = "" self.algoTypeToClass = {} + self.method = method # Store a python class of algorithm def add(self, algoClass): @@ -5199,25 +5200,48 @@ class algoCreator: # Create a copy of self and assign mesh to the copy def copy(self, mesh): - other = algoCreator() + other = algoCreator( self.method ) other.defaultAlgoType = self.defaultAlgoType - other.algoTypeToClass = self.algoTypeToClass + other.algoTypeToClass = self.algoTypeToClass other.mesh = mesh return other # Create an instance of algorithm def __call__(self,algo="",geom=0,*args): - algoType = self.defaultAlgoType - for arg in args + (algo,geom): - if isinstance( arg, geomBuilder.GEOM._objref_GEOM_Object ): - geom = arg - if isinstance( arg, str ) and arg: + algoType = "" + shape = 0 + if isinstance( algo, str ): + algoType = algo + elif ( isinstance( algo, geomBuilder.GEOM._objref_GEOM_Object ) and \ + not isinstance( geom, geomBuilder.GEOM._objref_GEOM_Object )): + shape = algo + elif algo: + args += (algo,) + + if isinstance( geom, geomBuilder.GEOM._objref_GEOM_Object ): + shape = geom + elif not algoType and isinstance( geom, str ): + algoType = geom + elif geom: + args += (geom,) + for arg in args: + if isinstance( arg, geomBuilder.GEOM._objref_GEOM_Object ) and not shape: + shape = arg + elif isinstance( arg, str ) and not algoType: algoType = arg + else: + import traceback, sys + msg = "Warning. Unexpected argument in mesh.%s() ---> %s" % ( self.method, arg ) + sys.stderr.write( msg + '\n' ) + tb = traceback.extract_stack(None,2) + traceback.print_list( [tb[0]] ) + if not algoType: + algoType = self.defaultAlgoType if not algoType and self.algoTypeToClass: algoType = self.algoTypeToClass.keys()[0] if self.algoTypeToClass.has_key( algoType ): #print "Create algo",algoType - return self.algoTypeToClass[ algoType ]( self.mesh, geom ) + return self.algoTypeToClass[ algoType ]( self.mesh, shape ) raise RuntimeError, "No class found for algo type %s" % algoType return None @@ -5299,7 +5323,7 @@ for pluginName in os.environ[ "SMESH_MeshersList" ].split( ":" ): if type( algo ).__name__ == 'classobj' and hasattr( algo, "meshMethod" ): #print " meshMethod:" , str(algo.meshMethod) if not hasattr( Mesh, algo.meshMethod ): - setattr( Mesh, algo.meshMethod, algoCreator() ) + setattr( Mesh, algo.meshMethod, algoCreator( algo.meshMethod )) pass getattr( Mesh, algo.meshMethod ).add( algo ) pass diff --git a/src/StdMeshers/StdMeshers_Prism_3D.cxx b/src/StdMeshers/StdMeshers_Prism_3D.cxx index 5ae0a4b0b..0ee4041ce 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.cxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.cxx @@ -52,7 +52,6 @@ #include #include #include -#include #include #include #include @@ -1177,6 +1176,9 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) { // use transformation (issue 0020680, IPAL0052499) or a "straight line" approach StdMeshers_Sweeper sweeper; + sweeper.myHelper = myHelper; + sweeper.myBotFace = thePrism.myBottom; + sweeper.myTopFace = thePrism.myTop; // load boundary nodes into sweeper bool dummy; @@ -1213,15 +1215,15 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) { double tol = getSweepTolerance( thePrism ); bool allowHighBndError = !isSimpleBottom( thePrism ); - myUseBlock = !sweeper.ComputeNodes( *myHelper, tol, allowHighBndError ); + myUseBlock = !sweeper.ComputeNodesByTrsf( tol, allowHighBndError ); } else if ( sweeper.CheckSameZ() ) { - myUseBlock = !sweeper.ComputeNodesOnStraightSameZ( *myHelper ); + myUseBlock = !sweeper.ComputeNodesOnStraightSameZ(); } else { - myUseBlock = !sweeper.ComputeNodesOnStraight( *myHelper, thePrism.myBottom, thePrism.myTop ); + myUseBlock = !sweeper.ComputeNodesOnStraight(); } myHelper->SetElementsOnShape( false ); } @@ -4956,13 +4958,13 @@ void StdMeshers_Sweeper::applyBoundaryError(const vector< gp_XYZ >& bndPoints, //================================================================================ /*! - * \brief Create internal nodes of the prism + * \brief Create internal nodes of the prism by computing an affine transformation + * from layer to layer */ //================================================================================ -bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, - const double tol, - const bool allowHighBndError) +bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, + const bool allowHighBndError) { const size_t zSize = myBndColumns[0]->size(); const size_t zSrc = 0, zTgt = zSize-1; @@ -5229,7 +5231,7 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, for ( size_t z = zSrc + 1; z < zTgt; ++z ) // vertical loop on layers { const gp_XYZ & xyz = intPntsOfLayer[ z ][ iP ]; - if ( !( nodeCol[ z ] = helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ))) + if ( !( nodeCol[ z ] = myHelper->AddNode( xyz.X(), xyz.Y(), xyz.Z() ))) return false; } } @@ -5301,7 +5303,7 @@ bool StdMeshers_Sweeper::CheckSameZ() */ //================================================================================ -bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper ) +bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ() { TZColumn& z = myZColumns[0]; @@ -5313,7 +5315,7 @@ bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper for ( size_t iZ = 0; iZ < z.size(); ++iZ ) { gp_XYZ p = n0 * ( 1 - z[iZ] ) + n1 * z[iZ]; - nodes[ iZ+1 ] = helper.AddNode( p.X(), p.Y(), p.Z() ); + nodes[ iZ+1 ] = myHelper->AddNode( p.X(), p.Y(), p.Z() ); } } @@ -5327,144 +5329,51 @@ bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper */ //================================================================================ -bool StdMeshers_Sweeper::ComputeNodesOnStraight( SMESH_MesherHelper& helper, - const TopoDS_Face& botFace, - const TopoDS_Face& topFace ) +bool StdMeshers_Sweeper::ComputeNodesOnStraight() { - // get data to create a Morph - UVPtStructVec botUV( myBndColumns.size() + 1 ); - UVPtStructVec topUV( myBndColumns.size() + 1 ); - for ( size_t i = 0; i < myBndColumns.size(); ++i ) - { - TNodeColumn& nodes = *myBndColumns[i]; - botUV[i].node = nodes[0]; - botUV[i].SetUV( helper.GetNodeUV( botFace, nodes[0] )); - topUV[i].node = nodes.back(); - topUV[i].SetUV( helper.GetNodeUV( topFace, nodes.back() )); - botUV[i].node->setIsMarked( true ); - } - botUV.back() = botUV[0]; - topUV.back() = topUV[0]; - - TopoDS_Edge dummyE; - TSideVector botWires( 1, StdMeshers_FaceSide::New( botUV, botFace, dummyE, helper.GetMesh() )); - TSideVector topWires( 1, StdMeshers_FaceSide::New( topUV, topFace, dummyE, helper.GetMesh() )); - - // use Morph to make delauney mesh on the FACEs. Locating of a node within a - // delauney triangle will be used to get a weighted Z. - NSProjUtils::Morph botDelauney( botWires ); - NSProjUtils::Morph topDelauney( topWires ); - - if ( helper.GetIsQuadratic() ) - { - // mark all medium nodes of faces on botFace to avoid their treating - SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( botFace ); - SMDS_ElemIteratorPtr eIt = smDS->GetElements(); - while ( eIt->more() ) - { - const SMDS_MeshElement* e = eIt->next(); - for ( int i = e->NbCornerNodes(), nb = e->NbNodes(); i < nb; ++i ) - e->GetNode( i )->setIsMarked( true ); - } - } + prepareTopBotDelaunay(); - // map to get a node column by a bottom node - TColStd_DataMapOfIntegerInteger iNode2iCol( myIntColumns.size() ); - - // un-mark nodes to treat (internal bottom nodes); later we will mark treated nodes - for ( size_t i = 0; i < myIntColumns.size(); ++i ) - { - const SMDS_MeshNode* botNode = myIntColumns[i]->front(); - botNode->setIsMarked( false ); - iNode2iCol.Bind( botNode->GetID(), i ); - } - - const int botFaceID = helper.GetMesh()->GetSubMesh( botFace )->GetId(); const SMDS_MeshNode *botNode, *topNode; - const BRepMesh_Triangle *botTria, *topTria; + const BRepMesh_Triangle *topTria; double botBC[3], topBC[3]; // barycentric coordinates int botTriaNodes[3], topTriaNodes[3]; bool checkUV = true; - // a queue of bottom nodes with starting delauney triangles - NSProjUtils::Morph::TNodeTriaList botNoTriQueue; + int nbInternalNodes = myIntColumns.size(); + myBotDelaunay->InitTraversal( nbInternalNodes ); - size_t iBndN = 1; // index of a bottom boundary node - int nbNodesToProcess = myIntColumns.size(); - while ( nbNodesToProcess > 0 ) + while (( botNode = myBotDelaunay->NextNode( botBC, botTriaNodes ))) { - while ( !botNoTriQueue.empty() ) // treat all nodes in the queue - { - botNode = botNoTriQueue.front().first; - botTria = botNoTriQueue.front().second; - botNoTriQueue.pop_front(); - if ( botNode->isMarked() ) - continue; - --nbNodesToProcess; - botNode->setIsMarked( true ); - - TNodeColumn* column = myIntColumns[ iNode2iCol( botNode->GetID() )]; - - // find a delauney triangle containing the botNode - gp_XY botUV = helper.GetNodeUV( botFace, botNode, NULL, &checkUV ); - botUV *= botDelauney.GetScale(); - botTria = botDelauney.FindTriangle( botUV, botTria, botBC, botTriaNodes ); - if ( !botTria ) - return false; - - // find a delauney triangle containing the topNode - topNode = column->back(); - gp_XY topUV = helper.GetNodeUV( topFace, topNode, NULL, &checkUV ); - topUV *= topDelauney.GetScale(); - // get a starting triangle basing on that top and bot boundary nodes have same index - topTria = topDelauney.GetTriangleNear( botTriaNodes[0] ); - topTria = topDelauney.FindTriangle( topUV, topTria, topBC, topTriaNodes ); - if ( !topTria ) - return false; - - // create nodes along a line - SMESH_NodeXYZ botP( botNode ), topP( topNode); - for ( size_t iZ = 0; iZ < myZColumns[0].size(); ++iZ ) - { - // use barycentric coordinates as weight of Z of boundary columns - double botZ = 0, topZ = 0; - for ( int i = 0; i < 3; ++i ) - { - botZ += botBC[i] * myZColumns[ botTriaNodes[i]-1 ][ iZ ]; - topZ += topBC[i] * myZColumns[ topTriaNodes[i]-1 ][ iZ ]; - } - double rZ = double( iZ + 1 ) / ( myZColumns[0].size() + 1 ); - double z = botZ * ( 1 - rZ ) + topZ * rZ; - gp_XYZ p = botP * ( 1 - z ) + topP * z; - (*column)[ iZ+1 ] = helper.AddNode( p.X(), p.Y(), p.Z() ); - } - - // add neighbor nodes to the queue - botDelauney.AddCloseNodes( botNode, botTria, botFaceID, botNoTriQueue ); - } + TNodeColumn* column = myIntColumns[ myNodeID2ColID( botNode->GetID() )]; + + // find a Delaunay triangle containing the topNode + topNode = column->back(); + gp_XY topUV = myHelper->GetNodeUV( myTopFace, topNode, NULL, &checkUV ); + // get a starting triangle basing on that top and bot boundary nodes have same index + topTria = myTopDelaunay->GetTriangleNear( botTriaNodes[0] ); + topTria = myTopDelaunay->FindTriangle( topUV, topTria, topBC, topTriaNodes ); + if ( !topTria ) + return false; - if ( nbNodesToProcess > 0 ) // fill the queue + // create nodes along a line + SMESH_NodeXYZ botP( botNode ), topP( topNode); + for ( size_t iZ = 0; iZ < myZColumns[0].size(); ++iZ ) { - // assure that all bot nodes are visited - for ( ; iBndN-1 < myBndColumns.size() && botNoTriQueue.empty(); ++iBndN ) - { - botTria = botDelauney.GetTriangleNear( iBndN ); - const SMDS_MeshNode* bndNode = botDelauney.GetBndNodes()[ iBndN ]; - botDelauney.AddCloseNodes( bndNode, botTria, botFaceID, botNoTriQueue ); - } - if ( botNoTriQueue.empty() ) + // use barycentric coordinates as weight of Z of boundary columns + double botZ = 0, topZ = 0; + for ( int i = 0; i < 3; ++i ) { - for ( size_t i = 0; i < myIntColumns.size(); ++i ) - { - botNode = myIntColumns[i]->front(); - if ( !botNode->isMarked() ) - botNoTriQueue.push_back( make_pair( botNode, botTria )); - } + botZ += botBC[i] * myZColumns[ botTriaNodes[i] ][ iZ ]; + topZ += topBC[i] * myZColumns[ topTriaNodes[i] ][ iZ ]; } + double rZ = double( iZ + 1 ) / ( myZColumns[0].size() + 1 ); + double z = botZ * ( 1 - rZ ) + topZ * rZ; + gp_XYZ p = botP * ( 1 - z ) + topP * z; + (*column)[ iZ+1 ] = myHelper->AddNode( p.X(), p.Y(), p.Z() ); } } - return true; + return myBotDelaunay->NbVisitedNodes() == nbInternalNodes; } //================================================================================ @@ -5490,3 +5399,58 @@ void StdMeshers_Sweeper::fillZColumn( TZColumn& zColumn, zColumn[i] = ( line * vec ) / len2; // param [0,1] on the line } } + +//================================================================================ +/*! + * \brief Initialize *Delaunay members + */ +//================================================================================ + +void StdMeshers_Sweeper::prepareTopBotDelaunay() +{ + UVPtStructVec botUV( myBndColumns.size() ); + UVPtStructVec topUV( myBndColumns.size() ); + for ( size_t i = 0; i < myBndColumns.size(); ++i ) + { + TNodeColumn& nodes = *myBndColumns[i]; + botUV[i].node = nodes[0]; + botUV[i].SetUV( myHelper->GetNodeUV( myBotFace, nodes[0] )); + topUV[i].node = nodes.back(); + topUV[i].SetUV( myHelper->GetNodeUV( myTopFace, nodes.back() )); + botUV[i].node->setIsMarked( true ); + } + TopoDS_Edge dummyE; + SMESH_Mesh* mesh = myHelper->GetMesh(); + TSideVector botWires( 1, StdMeshers_FaceSide::New( botUV, myBotFace, dummyE, mesh )); + TSideVector topWires( 1, StdMeshers_FaceSide::New( topUV, myTopFace, dummyE, mesh )); + + // Delaunay mesh on the FACEs. + bool checkUV = false; + myBotDelaunay.reset( new NSProjUtils::Delaunay( botWires, checkUV )); + myTopDelaunay.reset( new NSProjUtils::Delaunay( topWires, checkUV )); + + if ( myHelper->GetIsQuadratic() ) + { + // mark all medium nodes of faces on botFace to avoid their treating + SMESHDS_SubMesh* smDS = myHelper->GetMeshDS()->MeshElements( myBotFace ); + SMDS_ElemIteratorPtr eIt = smDS->GetElements(); + while ( eIt->more() ) + { + const SMDS_MeshElement* e = eIt->next(); + for ( int i = e->NbCornerNodes(), nb = e->NbNodes(); i < nb; ++i ) + e->GetNode( i )->setIsMarked( true ); + } + } + + // map to get a node column by a bottom node + myNodeID2ColID.Clear(/*doReleaseMemory=*/false); + myNodeID2ColID.ReSize( myIntColumns.size() ); + + // un-mark nodes to treat (internal bottom nodes) to be returned by myBotDelaunay + for ( size_t i = 0; i < myIntColumns.size(); ++i ) + { + const SMDS_MeshNode* botNode = myIntColumns[i]->front(); + botNode->setIsMarked( false ); + myNodeID2ColID.Bind( botNode->GetID(), i ); + } +} diff --git a/src/StdMeshers/StdMeshers_Prism_3D.hxx b/src/StdMeshers/StdMeshers_Prism_3D.hxx index ccdf5457f..51a9c4e71 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.hxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.hxx @@ -31,16 +31,18 @@ #include "SMESHDS_Mesh.hxx" #include "SMESH_Block.hxx" +#include "SMESH_Comment.hxx" #include "SMESH_Mesh.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_TypeDefs.hxx" -#include "SMESH_Comment.hxx" #include "SMESH_subMesh.hxx" +#include "StdMeshers_ProjectionUtils.hxx" #include #include #include #include +#include #include #include #include @@ -53,10 +55,6 @@ namespace Prism_3D struct TNode; struct TPrismTopo; } -namespace StdMeshers_ProjectionUtils -{ - class TrsfFinder3D; -} class SMESHDS_SubMesh; class TopoDS_Edge; @@ -417,23 +415,29 @@ private: */ struct StdMeshers_Sweeper { - std::vector< TNodeColumn* > myBndColumns; // boundary nodes - std::vector< TNodeColumn* > myIntColumns; // internal nodes + SMESH_MesherHelper* myHelper; + TopoDS_Face myBotFace; + TopoDS_Face myTopFace; + + std::vector< TNodeColumn* > myBndColumns; // boundary nodes + std::vector< TNodeColumn* > myIntColumns; // internal nodes typedef std::vector< double > TZColumn; - std::vector< TZColumn > myZColumns; // Z distribution of boundary nodes + std::vector< TZColumn > myZColumns; // Z distribution of boundary nodes + + StdMeshers_ProjectionUtils::DelaunayPtr myTopDelaunay; + StdMeshers_ProjectionUtils::DelaunayPtr myBotDelaunay; + TColStd_DataMapOfIntegerInteger myNodeID2ColID; - bool ComputeNodes( SMESH_MesherHelper& helper, - const double tol, - const bool allowHighBndError ); + + bool ComputeNodesByTrsf( const double tol, + const bool allowHighBndError ); bool CheckSameZ(); - bool ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper ); + bool ComputeNodesOnStraightSameZ(); - bool ComputeNodesOnStraight( SMESH_MesherHelper& helper, - const TopoDS_Face& bottom, - const TopoDS_Face& top); + bool ComputeNodesOnStraight(); private: @@ -459,6 +463,8 @@ private: static void fillZColumn( TZColumn& zColumn, TNodeColumn& nodes ); + + void prepareTopBotDelaunay(); }; // =============================================== diff --git a/src/StdMeshers/StdMeshers_ProjectionUtils.cxx b/src/StdMeshers/StdMeshers_ProjectionUtils.cxx index cda4346fa..11f2b5e97 100644 --- a/src/StdMeshers/StdMeshers_ProjectionUtils.cxx +++ b/src/StdMeshers/StdMeshers_ProjectionUtils.cxx @@ -478,6 +478,26 @@ namespace { return !allBndEdges.empty(); } + /*! + * \brief Convertor used in Delaunay constructor + */ + struct SideVector2UVPtStructVec + { + std::vector< const UVPtStructVec* > _uvVecs; + + SideVector2UVPtStructVec( const TSideVector& wires ) + { + _uvVecs.resize( wires.size() ); + for ( size_t i = 0; i < wires.size(); ++i ) + _uvVecs[ i ] = & wires[i]->GetUVPtStruct(); + } + + operator const std::vector< const UVPtStructVec* > & () const + { + return _uvVecs; + } + }; + } // namespace //======================================================================= @@ -2795,126 +2815,6 @@ namespace StdMeshers_ProjectionUtils return true; } - //================================================================================ - /*! - * \brief Add in-FACE nodes surrounding a given node to a queue - */ - //================================================================================ - - void Morph::AddCloseNodes( const SMDS_MeshNode* srcNode, - const BRepMesh_Triangle* bmTria, - const int srcFaceID, - TNodeTriaList & noTriQueue ) - { - // find in-FACE nodes - SMDS_ElemIteratorPtr elems = srcNode->GetInverseElementIterator(SMDSAbs_Face); - while ( elems->more() ) - { - const SMDS_MeshElement* elem = elems->next(); - if ( elem->getshapeId() == srcFaceID ) - { - for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i ) - { - const SMDS_MeshNode* n = elem->GetNode( i ); - if ( !n->isMarked() /*&& n->getshapeId() == srcFaceID*/ ) - noTriQueue.push_back( make_pair( n, bmTria )); - } - } - } - } - - //================================================================================ - /*! - * \brief Find a delauney triangle containing a given 2D point and return - * barycentric coordinates within the found triangle - */ - //================================================================================ - - const BRepMesh_Triangle* Morph::FindTriangle( const gp_XY& uv, - const BRepMesh_Triangle* bmTria, - double bc[3], - int triaNodes[3] ) - { - int nodeIDs[3]; - gp_XY nodeUVs[3]; - int linkIDs[3]; - Standard_Boolean ori[3]; - - while ( bmTria ) - { - // check bmTria - - _triaDS->ElementNodes( *bmTria, nodeIDs ); - nodeUVs[0] = _triaDS->GetNode( nodeIDs[0] ).Coord(); - nodeUVs[1] = _triaDS->GetNode( nodeIDs[1] ).Coord(); - nodeUVs[2] = _triaDS->GetNode( nodeIDs[2] ).Coord(); - - SMESH_MeshAlgos::GetBarycentricCoords( uv, - nodeUVs[0], nodeUVs[1], nodeUVs[2], - bc[0], bc[1] ); - if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 ) - { - bc[2] = 1 - bc[0] - bc[1]; - triaNodes[0] = nodeIDs[0]; - triaNodes[1] = nodeIDs[1]; - triaNodes[2] = nodeIDs[2]; - return bmTria; - } - - // look for a neighbor triangle, which is adjacent to a link intersected - // by a segment( triangle center -> uv ) - - gp_XY gc = ( nodeUVs[0] + nodeUVs[1] + nodeUVs[2] ) / 3.; - gp_XY seg = uv - gc; - - bmTria->Edges( linkIDs, ori ); - int triaID = _triaDS->IndexOf( *bmTria ); - bmTria = 0; - - for ( int i = 0; i < 3; ++i ) - { - const BRepMesh_PairOfIndex & triIDs = _triaDS->ElementsConnectedTo( linkIDs[i] ); - if ( triIDs.Extent() < 2 ) - continue; // no neighbor triangle - - // check if a link intersects gc2uv - const BRepMesh_Edge & link = _triaDS->GetLink( linkIDs[i] ); - const BRepMesh_Vertex & n1 = _triaDS->GetNode( link.FirstNode() ); - const BRepMesh_Vertex & n2 = _triaDS->GetNode( link.LastNode() ); - gp_XY uv1 = n1.Coord(); - gp_XY lin = n2.Coord() - uv1; // link direction - - double crossSegLin = seg ^ lin; - if ( Abs( crossSegLin ) < std::numeric_limits::min() ) - continue; // parallel - - double uSeg = ( uv1 - gc ) ^ lin / crossSegLin; - if ( 0. <= uSeg && uSeg <= 1. ) - { - bmTria = & _triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID ))); - break; - } - } - } - return bmTria; - } - - //================================================================================ - /*! - * \brief Return a triangle sharing a given boundary node - * \param [in] iBndNode - index of the boundary node - * \return const BRepMesh_Triangle* - a found triangle - */ - //================================================================================ - - const BRepMesh_Triangle* Morph::GetTriangleNear( int iBndNode ) - { - const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndNode ); - const BRepMesh_PairOfIndex & triaIds = _triaDS->ElementsConnectedTo( linkIds.First() ); - const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(1) ); - return &tria; - } - //================================================================================ /*! * \brief triangulate the srcFace in 2D @@ -2922,58 +2822,10 @@ namespace StdMeshers_ProjectionUtils */ //================================================================================ - Morph::Morph(const TSideVector& srcWires) + Morph::Morph(const TSideVector& srcWires): + _delaunay( srcWires, /*checkUV=*/true ) { _srcSubMesh = srcWires[0]->GetMesh()->GetSubMesh( srcWires[0]->Face() ); - - // compute _scale - { - BRepAdaptor_Surface surf( srcWires[0]->Face() ); - const int nbDiv = 100; - const double uRange = surf.LastUParameter() - surf.FirstUParameter(); - const double vRange = surf.LastVParameter() - surf.FirstVParameter(); - const double dU = uRange / nbDiv; - const double dV = vRange / nbDiv; - double u = surf.FirstUParameter(), v = surf.FirstVParameter(); - gp_Pnt p0U = surf.Value( u, v ), p0V = p0U; - double lenU = 0, lenV = 0; - for ( ; u < surf.LastUParameter(); u += dU, v += dV ) - { - gp_Pnt p1U = surf.Value( u, surf.FirstVParameter() ); - lenU += p1U.Distance( p0U ); - p0U = p1U; - gp_Pnt p1V = surf.Value( surf.FirstUParameter(), v ); - lenV += p1V.Distance( p0V ); - p0V = p1V; - } - _scale.SetCoord( lenU / uRange, lenV / vRange ); - } - - // count boundary points - int iP = 1, nbP = 0; - for ( size_t iW = 0; iW < srcWires.size(); ++iW ) - nbP += srcWires[iW]->NbPoints() - 1; // 1st and last points coincide - - _bndSrcNodes.resize( nbP + 1 ); _bndSrcNodes[0] = 0; - - // fill boundary points - BRepMesh::Array1OfVertexOfDelaun srcVert( 1, 1 + nbP ); - BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier ); - for ( size_t iW = 0; iW < srcWires.size(); ++iW ) - { - const UVPtStructVec& srcPnt = srcWires[iW]->GetUVPtStruct(); - for ( int i = 0, nb = srcPnt.size() - 1; i < nb; ++i, ++iP ) - { - _bndSrcNodes[ iP ] = srcPnt[i].node; - srcPnt[i].node->setIsMarked( true ); - - v.ChangeCoord() = srcPnt[i].UV().Multiplied( _scale ); - srcVert( iP ) = v; - } - } - // triangulate the srcFace in 2D - BRepMesh_Delaun delauney( srcVert ); - _triaDS = delauney.Result(); } //================================================================================ @@ -2994,32 +2846,27 @@ namespace StdMeshers_ProjectionUtils const TNodeNodeMap& src2tgtNodes, const bool moveAll) { - // get tgt boundary points corresponding to _bndSrcNodes + // get tgt boundary points corresponding to src boundary nodes size_t nbP = 0; for ( size_t iW = 0; iW < tgtWires.size(); ++iW ) nbP += tgtWires[iW]->NbPoints() - 1; // 1st and last points coincide - if ( nbP != _bndSrcNodes.size() - 1 ) + if ( nbP != _delaunay.GetBndNodes().size() ) return false; - BRepMesh::Array1OfVertexOfDelaun tgtVert( 1, 1 + nbP ); - BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier ); - for ( size_t iW = 0, iP = 1; iW < tgtWires.size(); ++iW ) + std::vector< gp_XY > tgtUV( nbP ); + for ( size_t iW = 0, iP = 0; iW < tgtWires.size(); ++iW ) { const UVPtStructVec& tgtPnt = tgtWires[iW]->GetUVPtStruct(); for ( int i = 0, nb = tgtPnt.size() - 1; i < nb; ++i, ++iP ) { - v.ChangeCoord() = tgtPnt[i].UV().Multiplied( _scale ); - tgtVert( iP ) = v; + tgtUV[ iP ] = tgtPnt[i].UV(); } } - const TopoDS_Face& srcFace = TopoDS::Face( _srcSubMesh->GetSubShape() ); - const int srcFaceID = _srcSubMesh->GetId(); SMESHDS_Mesh* tgtMesh = tgtHelper.GetMeshDS(); const SMDS_MeshNode *srcNode, *tgtNode; - const BRepMesh_Triangle *bmTria; - // un-mark internal src nodes; later we will mark moved nodes + // un-mark internal src nodes in order iterate them using _delaunay int nbSrcNodes = 0; SMDS_NodeIteratorPtr nIt = _srcSubMesh->GetSubMeshDS()->GetNodes(); if ( !nIt || !nIt->more() ) return true; @@ -3038,82 +2885,75 @@ namespace StdMeshers_ProjectionUtils // Move tgt nodes double bc[3]; // barycentric coordinates - int nodeIDs[3]; - bool checkUV = true; + int nodeIDs[3]; // nodes of a delaunay triangle const SMDS_FacePosition* pos; - // a queue of nodes with starting triangles - TNodeTriaList noTriQueue; - size_t iBndSrcN = 1; + _delaunay.InitTraversal( nbSrcNodes ); - while ( nbSrcNodes > 0 ) + while (( srcNode = _delaunay.NextNode( bc, nodeIDs ))) { - while ( !noTriQueue.empty() ) - { - srcNode = noTriQueue.front().first; - bmTria = noTriQueue.front().second; - noTriQueue.pop_front(); - if ( srcNode->isMarked() ) - continue; - --nbSrcNodes; - srcNode->setIsMarked( true ); - - // find a delauney triangle containing the src node - gp_XY uv = tgtHelper.GetNodeUV( srcFace, srcNode, NULL, &checkUV ); - uv *= _scale; - bmTria = FindTriangle( uv, bmTria, bc, nodeIDs ); - if ( !bmTria ) - continue; - - // compute new coordinates for a corresponding tgt node - gp_XY uvNew( 0., 0. ), nodeUV; - for ( int i = 0; i < 3; ++i ) - uvNew += bc[i] * tgtVert( nodeIDs[i]).Coord(); - uvNew.SetCoord( uvNew.X() / _scale.X(), - uvNew.Y() / _scale.Y() ); - gp_Pnt xyz = tgtSurface->Value( uvNew ); - - // find and move tgt node - TNodeNodeMap::const_iterator n2n = src2tgtNodes.find( srcNode ); - if ( n2n == src2tgtNodes.end() ) continue; - tgtNode = n2n->second; - tgtMesh->MoveNode( tgtNode, xyz.X(), xyz.Y(), xyz.Z() ); - - if (( pos = dynamic_cast< const SMDS_FacePosition* >( tgtNode->GetPosition() ))) - const_cast( pos )->SetParameters( uvNew.X(), uvNew.Y() ); - - AddCloseNodes( srcNode, bmTria, srcFaceID, noTriQueue ); - } + // compute new coordinates for a corresponding tgt node + gp_XY uvNew( 0., 0. ), nodeUV; + for ( int i = 0; i < 3; ++i ) + uvNew += bc[i] * tgtUV[ nodeIDs[i]]; + gp_Pnt xyz = tgtSurface->Value( uvNew ); - if ( nbSrcNodes > 0 ) - { - // assure that all src nodes are visited - for ( ; iBndSrcN < _bndSrcNodes.size() && noTriQueue.empty(); ++iBndSrcN ) - { - const BRepMesh_Triangle* tria = GetTriangleNear( iBndSrcN ); - AddCloseNodes( _bndSrcNodes[ iBndSrcN ], tria, srcFaceID, noTriQueue ); - } - if ( noTriQueue.empty() ) - { - SMDS_NodeIteratorPtr nIt = _srcSubMesh->GetSubMeshDS()->GetNodes(); - while ( nIt->more() ) - { - srcNode = nIt->next(); - if ( !srcNode->isMarked() ) - noTriQueue.push_back( make_pair( srcNode, bmTria )); - } - } - } + // find and move tgt node + TNodeNodeMap::const_iterator n2n = src2tgtNodes.find( srcNode ); + if ( n2n == src2tgtNodes.end() ) continue; + tgtNode = n2n->second; + tgtMesh->MoveNode( tgtNode, xyz.X(), xyz.Y(), xyz.Z() ); + + if (( pos = dynamic_cast< const SMDS_FacePosition* >( tgtNode->GetPosition() ))) + const_cast( pos )->SetParameters( uvNew.X(), uvNew.Y() ); + + --nbSrcNodes; } - return true; + return nbSrcNodes == 0; } // Morph::Perform - gp_XY Morph::GetBndUV(const int iNode) const + //======================================================================= + //function : Delaunay + //purpose : construct from face sides + //======================================================================= + + Delaunay::Delaunay( const TSideVector& wires, bool checkUV ): + SMESH_Delaunay( SideVector2UVPtStructVec( wires ), + TopoDS::Face( wires[0]->FaceHelper()->GetSubShape() ), + wires[0]->FaceHelper()->GetSubShapeID() ) { - return _triaDS->GetNode( iNode ).Coord(); + _wire = wires[0]; // keep a wire to assure _helper to keep alive + _helper = _wire->FaceHelper(); + _checkUVPtr = checkUV ? & _checkUV : 0; } + //======================================================================= + //function : Delaunay + //purpose : construct from UVPtStructVec's + //======================================================================= + + Delaunay::Delaunay( const std::vector< const UVPtStructVec* > & boundaryNodes, + SMESH_MesherHelper& faceHelper, + bool checkUV): + SMESH_Delaunay( boundaryNodes, + TopoDS::Face( faceHelper.GetSubShape() ), + faceHelper.GetSubShapeID() ) + { + _helper = & faceHelper; + _checkUVPtr = checkUV ? & _checkUV : 0; + } + + //======================================================================= + //function : getNodeUV + //purpose : + //======================================================================= + + gp_XY Delaunay::getNodeUV( const TopoDS_Face& face, const SMDS_MeshNode* node ) const + { + return _helper->GetNodeUV( face, node, 0, _checkUVPtr ); + } + } // namespace StdMeshers_ProjectionUtils diff --git a/src/StdMeshers/StdMeshers_ProjectionUtils.hxx b/src/StdMeshers/StdMeshers_ProjectionUtils.hxx index 6c94e45d0..aa3828443 100644 --- a/src/StdMeshers/StdMeshers_ProjectionUtils.hxx +++ b/src/StdMeshers/StdMeshers_ProjectionUtils.hxx @@ -30,10 +30,10 @@ #include "SMESH_StdMeshers.hxx" -#include "StdMeshers_FaceSide.hxx" #include "SMDS_MeshElement.hxx" +#include "SMESH_Delaunay.hxx" +#include "StdMeshers_FaceSide.hxx" -#include #include #include #include @@ -54,6 +54,7 @@ class SMESH_Mesh; class SMESH_subMesh; class TopoDS_Shape; +//----------------------------------------------------------------------------------------- /*! * \brief Struct used instead of a sole TopTools_DataMapOfShapeShape to avoid * problems with bidirectional bindings @@ -94,6 +95,7 @@ namespace StdMeshers_ProjectionUtils TIDCompare> TNodeNodeMap; + //----------------------------------------------------------------------------------------- /*! * \brief Finds transformation between two sets of 2D points using * a least square approximation @@ -114,6 +116,7 @@ namespace StdMeshers_ProjectionUtils bool IsIdentity() const { return ( _trsf.Form() == gp_Identity ); } }; + //----------------------------------------------------------------------------------------- /*! * \brief Finds transformation between two sets of 3D points using * a least square approximation @@ -139,16 +142,44 @@ namespace StdMeshers_ProjectionUtils bool Invert(); }; + //----------------------------------------------------------------------------------------- + /*! + * \brief Create a Delaunay triangulation of nodes on a face boundary + * and provide exploration of nodes shared by elements lying on + * the face. For a returned node, also return a Delaunay triangle + * the node lies in and its Barycentric Coordinates within the triangle. + * Only non-marked nodes are visited. + * + * The main methods are defined in ../SMESHUtils/SMESH_Delaunay.hxx + */ + class Delaunay : public SMESH_Delaunay + { + public: + + Delaunay( const TSideVector& wires, bool checkUV = false ); + + Delaunay( const std::vector< const UVPtStructVec* > & boundaryNodes, + SMESH_MesherHelper& faceHelper, + bool checkUV = false); + + protected: + virtual gp_XY getNodeUV( const TopoDS_Face& face, const SMDS_MeshNode* node ) const; + + private: + SMESH_MesherHelper* _helper; + StdMeshers_FaceSidePtr _wire; + bool *_checkUVPtr, _checkUV; + }; + typedef boost::shared_ptr< Delaunay > DelaunayPtr; + + //----------------------------------------------------------------------------------------- /*! * \brief Morph mesh on the target FACE to lie within FACE boundary w/o distortion */ class Morph { - std::vector< const SMDS_MeshNode* > _bndSrcNodes; - Handle(BRepMesh_DataStructureOfDelaun) _triaDS; - SMESH_subMesh* _srcSubMesh; - bool _moveAll; - gp_XY _scale; + Delaunay _delaunay; + SMESH_subMesh* _srcSubMesh; public: Morph(const TSideVector& srcWires); @@ -158,37 +189,9 @@ namespace StdMeshers_ProjectionUtils Handle(ShapeAnalysis_Surface) tgtSurface, const TNodeNodeMap& src2tgtNodes, const bool moveAll); - - // find a triangle containing an UV starting from a given triangle; - // return barycentric coordinates of the UV in the found triangle - const BRepMesh_Triangle* FindTriangle( const gp_XY& uv, - const BRepMesh_Triangle* bmTria, - double bc[3], - int triaNodes[3]); - - // return any delauney triangle neighboring a given boundary node - const BRepMesh_Triangle* GetTriangleNear( int iBndNode ); - - // return source boundary nodes. 0-th node is NULL so that indices of - // boundary nodes correspond to indices used by Delauney mesh - const std::vector< const SMDS_MeshNode* >& GetBndNodes() const { return _bndSrcNodes; } - - // return UV of the i-th source boundary node - gp_XY GetBndUV(const int iNode) const; - - // return scale factor to convert real UV to/from UV used for Delauney meshing: - // delauney_UV = real_UV * scale - const gp_XY& GetScale() const { return _scale; } - - typedef std::list< std::pair< const SMDS_MeshNode*, const BRepMesh_Triangle* > > TNodeTriaList; - - // add non-marked nodes surrounding a given one to a list - static void AddCloseNodes( const SMDS_MeshNode* node, - const BRepMesh_Triangle* bmTria, - const int faceID, - TNodeTriaList & noTriQueue ); }; + //----------------------------------------------------------------------------------------- /*! * \brief Looks for association of all sub-shapes of two shapes * \param theShape1 - shape 1