From: eap Date: Tue, 25 Aug 2015 14:28:48 +0000 (+0300) Subject: Fix MA construction X-Git-Tag: V7_8_0a2~25^2~2 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=refs%2Fheads%2Fhydro%2Fimps_2015_V760;p=modules%2Fsmesh.git Fix MA construction --- diff --git a/src/SMESHUtils/SMESH_MAT2d.cxx b/src/SMESHUtils/SMESH_MAT2d.cxx index dcde99881..66ac2ff30 100644 --- a/src/SMESHUtils/SMESH_MAT2d.cxx +++ b/src/SMESHUtils/SMESH_MAT2d.cxx @@ -47,6 +47,7 @@ #include #ifdef _DEBUG_ +//#define _MYDEBUG_ #include "SMESH_File.hxx" #include "SMESH_Comment.hxx" #endif @@ -66,8 +67,8 @@ namespace struct InPoint { - int _a, _b; - double _param; + int _a, _b; // coordinates + double _param; // param on EDGE InPoint(int x, int y, double param) : _a(x), _b(y), _param(param) {} InPoint() : _a(0), _b(0), _param(0) {} @@ -76,6 +77,8 @@ namespace size_t index( const vector< InPoint >& inPoints ) const { return this - &inPoints[0]; } bool operator==( const InPoint& other ) const { return _a == other._a && _b == other._b; } + bool operator==( const TVDVertex* v ) const { return ( Abs( _a - v->x() ) < 1. && + Abs( _b - v->y() ) < 1. ); } }; // ------------------------------------------------------------------------------------- @@ -105,10 +108,12 @@ namespace // check if a TVDEdge begins at my end or ends at my start inline bool InSegment::isConnected( const TVDEdge* edge ) { - return ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&& - Abs( edge->vertex0()->y() - _p1->_b ) < 1. ) || - (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&& - Abs( edge->vertex1()->y() - _p0->_b ) < 1. )); + return (( edge->vertex0() && edge->vertex1() ) + && + ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&& + Abs( edge->vertex0()->y() - _p1->_b ) < 1. ) || + (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&& + Abs( edge->vertex1()->y() - _p0->_b ) < 1. ))); } // check if a MA TVDEdge is outside of a domain @@ -147,7 +152,7 @@ namespace // } // ------------------------------------------------------------------------------------- -#ifdef _DEBUG_ +#ifdef _MYDEBUG_ // writes segments into a txt file readable by voronoi_visualizer void inSegmentsToFile( vector< InSegment>& inSegments) { @@ -155,6 +160,7 @@ namespace return; const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAdebug.txt"; SMESH_File file(fileName, false ); + file.remove(); file.openForWriting(); SMESH_Comment text; text << "0\n"; // nb points @@ -180,7 +186,7 @@ namespace if ( !edge->vertex1() ) cout << ") -> ( INF, INF"; else - cout << ") -> (" << edge->vertex1()->x() << ", " << edge->vertex1()->y(); + cout << ") -> ( " << edge->vertex1()->x() << ", " << edge->vertex1()->y(); cout << ")\t cell=" << edge->cell() << " iBnd=" << edge->color() << " twin=" << edge->twin() @@ -253,6 +259,7 @@ namespace boost { namespace { const int theNoBrachID = 0; // std::numeric_limits::max(); + double theScale[2]; // scale used in bndSegsToMesh() // ------------------------------------------------------------------------------------- /*! @@ -298,7 +305,7 @@ namespace if ( _edge ) // pass branch to an opposite BndSeg { size_t oppSegIndex = SMESH_MAT2d::Branch::getBndSegment( _edge->twin() ); - if ( oppSegIndex < bndSegs.size() /*&& bndSegs[ oppSegIndex ]._branchID == theNoBrachID*/ ) + if ( oppSegIndex < bndSegs.size() && bndSegs[ oppSegIndex ]._branchID == theNoBrachID ) bndSegs[ oppSegIndex ]._branchID = -branchID; } } @@ -349,7 +356,64 @@ namespace //================================================================================ /*! - * \brief Computes length of of TVDEdge + * \brief debug: to visually check found MA edges + */ + //================================================================================ + + void bndSegsToMesh( const vector< BndSeg >& bndSegs ) + { +#ifdef _MYDEBUG_ + if ( !getenv("bndSegsToMesh")) return; + map< const TVDVertex *, int > v2Node; + map< const TVDVertex *, int >::iterator v2n; + set< const TVDEdge* > addedEdges; + + const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAedges.py"; + SMESH_File file(fileName, false ); + file.remove(); + file.openForWriting(); + SMESH_Comment text; + text << "import salome, SMESH\n"; + text << "salome.salome_init()\n"; + text << "from salome.smesh import smeshBuilder\n"; + text << "smesh = smeshBuilder.New(salome.myStudy)\n"; + text << "m=smesh.Mesh()\n"; + for ( size_t i = 0; i < bndSegs.size(); ++i ) + { + if ( !bndSegs[i]._edge ) + text << "# " << i << " NULL edge"; + else if ( !bndSegs[i]._edge->vertex0() || + !bndSegs[i]._edge->vertex1() ) + text << "# " << i << " INFINITE edge"; + else if ( addedEdges.insert( bndSegs[i]._edge ).second && + addedEdges.insert( bndSegs[i]._edge->twin() ).second ) + { + v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex0(), v2Node.size() + 1 )).first; + int n0 = v2n->second; + if ( n0 == v2Node.size() ) + text << "n" << n0 << " = m.AddNode( " + << bndSegs[i]._edge->vertex0()->x() / theScale[0] << ", " + << bndSegs[i]._edge->vertex0()->y() / theScale[1] << ", 0 )\n"; + + v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex1(), v2Node.size() + 1 )).first; + int n1 = v2n->second; + if ( n1 == v2Node.size() ) + text << "n" << n1 << " = m.AddNode( " + << bndSegs[i]._edge->vertex1()->x() / theScale[0] << ", " + << bndSegs[i]._edge->vertex1()->y() / theScale[1] << ", 0 )\n"; + + text << "e" << i << " = m.AddEdge([ n" << n0 << ", n" << n1 << " ])\n"; + } + } + text << "\n"; + file.write( text.c_str(), text.size() ); + cout << "Write " << fileName << endl; +#endif + } + + //================================================================================ + /*! + * \brief Computes length of a TVDEdge */ //================================================================================ @@ -459,24 +523,25 @@ namespace gp_Pnt p, pPrev; if ( !c3d.IsNull() ) pPrev = c3d->Value( f ); - for ( int i = 2; i <= discret.NbPoints(); i++ ) // skip the 1st point - { - double u = discret.Parameter(i); - if ( !c3d.IsNull() ) + if ( discret.NbPoints() > 2 ) + for ( int i = 2; i <= discret.NbPoints(); i++ ) // skip the 1st point { - p = c3d->Value( u ); - int nbDiv = int( p.Distance( pPrev ) / minSegLen / theDiscrCoef ); - double dU = ( u - points.back()._u ) / nbDiv; - for ( int iD = 1; iD < nbDiv; ++iD ) + double u = discret.Parameter(i); + if ( !c3d.IsNull() ) { - double uD = points.back()._u + dU; - points.push_back( UVU( c2d->Value( uD ), uD )); + p = c3d->Value( u ); + int nbDiv = int( p.Distance( pPrev ) / minSegLen / theDiscrCoef ); + double dU = ( u - points.back()._u ) / nbDiv; + for ( int iD = 1; iD < nbDiv; ++iD ) + { + double uD = points.back()._u + dU; + points.push_back( UVU( c2d->Value( uD ), uD )); + } + pPrev = p; } - pPrev = p; + points.push_back( UVU( c2d->Value( u ), u )); + uvBox.Add( points.back()._uv ); } - points.push_back( UVU( c2d->Value( u ), u )); - uvBox.Add( points.back()._uv ); - } // if ( !c3d.IsNull() ) // { // vector params; @@ -592,9 +657,52 @@ namespace } } } + // debug + theScale[0] = scale[0]; + theScale[1] = scale[1]; + return true; } + //================================================================================ + /*! + * \brief Update a branch joined to another one + */ + //================================================================================ + + void updateJoinedBranch( vector< const TVDEdge* > & branchEdges, + const size_t newID, + vector< BndSeg > & bndSegs, + const bool reverse) + { + if ( reverse ) + { + for ( size_t i = 0; i < branchEdges.size(); ++i ) + { + size_t seg1 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i] ); + size_t seg2 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i]->twin() ); + bndSegs[ seg1 ]._branchID /= bndSegs[ seg1 ].branchID(); + bndSegs[ seg2 ]._branchID /= bndSegs[ seg2 ].branchID(); + bndSegs[ seg1 ]._branchID *= -newID; + bndSegs[ seg2 ]._branchID *= -newID; + branchEdges[i] = branchEdges[i]->twin(); + } + std::reverse( branchEdges.begin(), branchEdges.end() ); + } + else + { + for ( size_t i = 0; i < branchEdges.size(); ++i ) + { + size_t seg1 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i] ); + size_t seg2 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i]->twin() ); + bndSegs[ seg1 ]._branchID /= bndSegs[ seg1 ].branchID(); + bndSegs[ seg2 ]._branchID /= bndSegs[ seg2 ].branchID(); + bndSegs[ seg1 ]._branchID *= newID; + bndSegs[ seg2 ]._branchID *= newID; + } + } + } + //================================================================================ /*! * \brief Create MA branches and FACE boundary data @@ -608,6 +716,7 @@ namespace //================================================================================ void makeMA( const TVD& vd, + const bool ignoreCorners, vector< InPoint >& inPoints, vector< InSegment > & inSegments, vector< SMESH_MAT2d::Branch >& branch, @@ -688,9 +797,12 @@ namespace continue; inPntChecked[ pInd ] = true; - const TVDEdge* edge = // a TVDEdge passing through an end of inSeg - is2nd ? maEdges.front()->prev() : maEdges.back()->next(); - while ( true ) + const TVDEdge* maE = is2nd ? maEdges.front() : maEdges.back(); + if ( inPnt == ( is2nd ? maE->vertex0() : maE->vertex1() )) + continue; + const TVDEdge* edge = // a secondary TVDEdge connecting inPnt and maE + is2nd ? maE->prev() : maE->next(); + while ( inSeg.isConnected( edge )) { if ( edge->is_primary() ) break; // this should not happen const TVDEdge* edge2 = edge->twin(); // we are in a neighbor cell, add MA edges to inPnt @@ -732,7 +844,7 @@ namespace inPoints[0]._edges.clear(); } - // Divide InSegment's into BndSeg's + // Divide InSegment's into BndSeg's (so that each BndSeg corresponds to one MA edge) vector< BndSeg > bndSegs; bndSegs.reserve( inSegments.size() * 3 ); @@ -750,25 +862,26 @@ namespace inPntChecked[ ip0 ] = false; // segments of InSegment's - size_t nbMaEdges = inSeg._edges.size(); + const size_t nbMaEdges = inSeg._edges.size(); switch ( nbMaEdges ) { case 0: // "around" circle center bndSegs.push_back( BndSeg( &inSeg, 0, inSeg._p1->_param )); break; case 1: bndSegs.push_back( BndSeg( &inSeg, inSeg._edges.back(), inSeg._p1->_param )); break; default: - vector< double > len; - len.push_back(0); - for ( e = inSeg._edges.rbegin(); e != inSeg._edges.rend(); ++e ) - len.push_back( len.back() + length( *e )); - + gp_XY inSegDir( inSeg._p1->_a - inSeg._p0->_a, + inSeg._p1->_b - inSeg._p0->_b ); + const double inSegLen2 = inSegDir.SquareModulus(); e = inSeg._edges.rbegin(); - for ( size_t l = 1; l < len.size(); ++e, ++l ) + for ( size_t iE = 1; iE < nbMaEdges; ++e, ++iE ) { - double dl = len[l] / len.back(); - double u = dl * inSeg._p1->_param + ( 1. - dl ) * inSeg._p0->_param; + gp_XY toMA( (*e)->vertex0()->x() - inSeg._p0->_a, + (*e)->vertex0()->y() - inSeg._p0->_b ); + double r = toMA * inSegDir / inSegLen2; + double u = r * inSeg._p1->_param + ( 1. - r ) * inSeg._p0->_param; bndSegs.push_back( BndSeg( &inSeg, *e, u )); } + bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param )); } // segments around 2nd concave point size_t ip1 = inSeg._p1->index( inPoints ); @@ -783,6 +896,8 @@ namespace for ( size_t i = 0; i < bndSegs.size(); ++i ) bndSegs[i].setIndexToEdge( i ); + bndSegsToMesh( bndSegs ); // debug: visually check found MA edges + // Find TVDEdge's of Branches and associate them with bndSegs @@ -797,7 +912,7 @@ namespace size_t i1st = 0; while ( i1st < bndSegs.size() && !bndSegs[i1st].hasOppositeEdge( noEdgeID )) ++i1st; - bndSegs[i1st].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg + bndSegs[i1st].setBranch( branchID, bndSegs ); // set to the i-th and to the opposite bndSeg branchEdges[ branchID ].push_back( bndSegs[i1st]._edge ); for ( size_t i = i1st+1; i < bndSegs.size(); ++i ) @@ -824,7 +939,7 @@ namespace endType.insert( make_pair( bndSegs[i]._edge->vertex1(), SMESH_MAT2d::BE_BRANCH_POINT )); } - bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg + bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and to the opposite bndSeg if ( bndSegs[i].hasOppositeEdge( noEdgeID )) branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge ); } @@ -858,17 +973,92 @@ namespace br2.clear(); } + // remove branches ending at BE_ON_VERTEX + + vector isBranchRemoved( branchEdges.size(), false ); + + if ( ignoreCorners && branchEdges.size() > 2 && !branchEdges[2].empty() ) + { + // find branches to remove + map< const TVDVertex*, SMESH_MAT2d::BranchEndType >::iterator v2et; + for ( size_t iB = 1; iB < branchEdges.size(); ++iB ) + { + if ( branchEdges[iB].empty() ) + continue; + const TVDVertex* v0 = branchEdges[iB][0]->vertex1(); + const TVDVertex* v1 = branchEdges[iB].back()->vertex0(); + v2et = endType.find( v0 ); + if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX ) + isBranchRemoved[ iB ] = true; + v2et = endType.find( v1 ); + if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX ) + isBranchRemoved[ iB ] = true; + } + // try to join not removed branches into one + for ( size_t iB = 1; iB < branchEdges.size(); ++iB ) + { + if ( branchEdges[iB].empty() || isBranchRemoved[iB] ) + continue; + const TVDVertex* v0 = branchEdges[iB][0]->vertex1(); + const TVDVertex* v1 = branchEdges[iB].back()->vertex0(); + v2et = endType.find( v0 ); + if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT ) + v0 = 0; + v2et = endType.find( v1 ); + if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT ) + v1 = 0; + if ( !v0 && !v1 ) + continue; + + size_t iBrToJoin = 0; + for ( size_t iB2 = 1; iB2 < branchEdges.size(); ++iB2 ) + { + if ( branchEdges[iB2].empty() || isBranchRemoved[iB2] || iB == iB2 ) + continue; + const TVDVertex* v02 = branchEdges[iB2][0]->vertex1(); + const TVDVertex* v12 = branchEdges[iB2].back()->vertex0(); + if ( v0 == v02 || v0 == v12 || v1 == v02 || v1 == v12 ) + { + if ( iBrToJoin > 0 ) + { + iBrToJoin = 0; + break; // more than 2 not removed branches meat at a TVDVertex + } + iBrToJoin = iB2; + } + } + if ( iBrToJoin > 0 ) + { + vector& branch = branchEdges[ iBrToJoin ]; + const TVDVertex* v02 = branch[0]->vertex1(); + const TVDVertex* v12 = branch.back()->vertex0(); + updateJoinedBranch( branch, iB, bndSegs, /*reverse=*/(v0 == v02 || v1 == v12 )); + if ( v0 == v02 || v0 == v12 ) + branchEdges[iB].insert( branchEdges[iB].begin(), branch.begin(), branch.end() ); + else + branchEdges[iB].insert( branchEdges[iB].end(), branch.begin(), branch.end() ); + branch.clear(); + } + } // loop on branchEdges + } // if ( ignoreCorners ) + // associate branchIDs and the input branch vector (arg) - vector< const SMESH_MAT2d::Branch* > branchByID( branchEdges.size() ); + vector< SMESH_MAT2d::Branch* > branchByID( branchEdges.size(), 0 ); int nbBranches = 0; for ( size_t i = 0; i < branchEdges.size(); ++i ) { nbBranches += ( !branchEdges[i].empty() ); } branch.resize( nbBranches ); - for ( size_t iBr = 0, brID = 0; brID < branchEdges.size(); ++brID ) + size_t iBr = 0; + for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // 1st - not removed + { + if ( !branchEdges[ brID ].empty() && !isBranchRemoved[ brID ]) + branchByID[ brID ] = & branch[ iBr++ ]; + } + for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // then - removed { - if ( !branchEdges[ brID ].empty() ) + if ( !branchEdges[ brID ].empty() && isBranchRemoved[ brID ]) branchByID[ brID ] = & branch[ iBr++ ]; } @@ -912,28 +1102,24 @@ namespace { edgeInd += dInd; - if ( edgeInd < 0 || - edgeInd >= (int) branchEdges[ brID ].size() || - branchEdges[ brID ][ edgeInd ] != bndSegs[ i ]._edge ) + if (( edgeInd < 0 || + edgeInd >= (int) branchEdges[ brID ].size() ) || + ( branchEdges[ brID ][ edgeInd ] != bndSegs[ i ]._edge && + branchEdges[ brID ][ edgeInd ]->twin() != bndSegs[ i ]._edge )) { - if ( bndSegs[ i ]._branchID < 0 && - branchEdges[ brID ].back() == bndSegs[ i ]._edge ) + if ( bndSegs[ i ]._branchID < 0 ) { - edgeInd = branchEdges[ brID ].size() - 1; dInd = -1; + for ( edgeInd = branchEdges[ brID ].size() - 1; edgeInd > 0; --edgeInd ) + if ( branchEdges[ brID ][ edgeInd ]->twin() == bndSegs[ i ]._edge ) + break; } - else if ( bndSegs[ i ]._branchID > 0 && - branchEdges[ brID ].front() == bndSegs[ i ]._edge ) + else // bndSegs[ i ]._branchID > 0 { - edgeInd = 0; dInd = +1; - } - else - { for ( edgeInd = 0; edgeInd < branchEdges[ brID ].size(); ++edgeInd ) if ( branchEdges[ brID ][ edgeInd ] == bndSegs[ i ]._edge ) break; - dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1; } } } @@ -941,15 +1127,10 @@ namespace { // no MA edge, bndSeg corresponds to an end point of a branch if ( bndPoints._maEdges.empty() ) - { - // should not get here according to algo design edgeInd = 0; - } else - { edgeInd = branchEdges[ brID ].size(); - dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1; - } + dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1; } bndPoints._maEdges.push_back( make_pair( br, ( 1 + edgeInd ) * dInd )); @@ -958,19 +1139,40 @@ namespace iSeg = iSegEnd; - } // loop on all bndSegs + } // loop on all bndSegs to construct Boundary + // Initialize branches + // find a not removed branch + size_t iBrNorRemoved = 0; + for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) + if ( !branchEdges[brID].empty() && !isBranchRemoved[brID] ) + { + iBrNorRemoved = brID; + break; + } // fill the branches with MA edges - for ( size_t iBr = 0, brID = 0; brID < branchEdges.size(); ++brID ) + for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) if ( !branchEdges[brID].empty() ) { - branch[ iBr ].init( branchEdges[brID], & boundary, endType ); - iBr++; + branchByID[ brID ]->init( branchEdges[brID], & boundary, endType ); + } + // mark removed branches + for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) + if ( isBranchRemoved[brID] && iBrNorRemoved > 0 ) + { + SMESH_MAT2d::Branch* branch = branchByID[ brID ]; + SMESH_MAT2d::Branch* mainBranch = branchByID[ iBrNorRemoved ]; + bool is1stBrPnt = ( branch->getEnd(0)->_type == SMESH_MAT2d::BE_BRANCH_POINT ); + const TVDVertex* branchVextex = + is1stBrPnt ? branch->getEnd(0)->_vertex : branch->getEnd(1)->_vertex; + SMESH_MAT2d::BranchPoint bp = mainBranch->getPoint( branchVextex ); + branch->setRemoved( bp ); } // set branches to branch ends for ( size_t i = 0; i < branch.size(); ++i ) - branch[i].setBranchesToEnds( branch ); + if ( !branch[i].isRemoved() ) + branch[i].setBranchesToEnds( branch ); // fill branchPnt arg map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* > v2end; @@ -1013,13 +1215,30 @@ SMESH_MAT2d::MedialAxis::MedialAxis(const TopoDS_Face& face, if ( !makeInputData( face, edges, minSegLen, inPoints, inSegments, _scale )) return; - //inSegmentsToFile( inSegments ); + inSegmentsToFile( inSegments ); // build voronoi diagram construct_voronoi( inSegments.begin(), inSegments.end(), &_vd ); // make MA data - makeMA( _vd, inPoints, inSegments, _branch, _branchPnt, _boundary ); + makeMA( _vd, ignoreCorners, inPoints, inSegments, _branch, _branchPnt, _boundary ); + + // count valid branches + _nbBranches = _branch.size(); + for ( size_t i = 0; i < _branch.size(); ++i ) + if ( _branch[i].isRemoved() ) + --_nbBranches; +} + +//================================================================================ +/*! + * \brief Returns the i-th branch + */ +//================================================================================ + +const SMESH_MAT2d::Branch* SMESH_MAT2d::MedialAxis::getBranch(size_t i) const +{ + return i < _nbBranches ? &_branch[i] : 0; } //================================================================================ @@ -1028,16 +1247,16 @@ SMESH_MAT2d::MedialAxis::MedialAxis(const TopoDS_Face& face, */ //================================================================================ -void SMESH_MAT2d::MedialAxis::getPoints( const Branch& branch, +void SMESH_MAT2d::MedialAxis::getPoints( const Branch* branch, std::vector< gp_XY >& points) const { - branch.getPoints( points, _scale ); + branch->getPoints( points, _scale ); } //================================================================================ /*! * \brief Returns a BranchPoint corresponding to a given point on a geom EDGE - * \param [in] iGeomEdge - index of geom EDGE within a vector passed at construction + * \param [in] iEdge - index of geom EDGE within a vector passed at MA construction * \param [in] u - parameter of the point on EDGE curve * \param [out] p - the found BranchPoint * \return bool - is OK @@ -1052,7 +1271,7 @@ bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge, return false; const BndPoints& points = _pointsPerEdge[ iEdge ]; - const bool edgeReverse = ( points._params[0] > points._params.back() ); + const bool edgeReverse = ( points._params[0] > points._params.back() ); if ( u < ( edgeReverse ? points._params.back() : points._params[0] )) u = edgeReverse ? points._params.back() : points._params[0]; @@ -1071,14 +1290,30 @@ bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge, while ( points._params[i ] > u ) --i; while ( points._params[i+1] < u ) ++i; } + double edgeParam = ( u - points._params[i] ) / ( points._params[i+1] - points._params[i] ); + if ( !points._maEdges[ i ].second ) // no branch at the EDGE end, look for a closest branch + { + if ( i < points._maEdges.size() / 2 ) // near 1st point + { + while ( i < points._maEdges.size()-1 && !points._maEdges[ i ].second ) + ++i; + edgeParam = edgeReverse; + } + else // near last point + { + while ( i > 0 && !points._maEdges[ i ].second ) + --i; + edgeParam = !edgeReverse; + } + } const std::pair< const Branch*, int >& maE = points._maEdges[ i ]; bool maReverse = ( maE.second < 0 ); - p._branch = maE.first; - p._iEdge = maE.second - 1; // countered from 1 to store sign - p._edgeParam = maReverse ? ( 1. - edgeParam ) : edgeParam; + p._branch = maE.first; + p._iEdge = ( maReverse ? -maE.second : maE.second ) - 1; // countered from 1 to store sign + p._edgeParam = ( maE.first && maReverse ) ? ( 1. - edgeParam ) : edgeParam; return true; } @@ -1093,16 +1328,36 @@ bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge, */ //================================================================================ -bool SMESH_MAT2d::Boundary::IsConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const +bool SMESH_MAT2d::Boundary::isConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const { if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() ) return false; const BndPoints& points = _pointsPerEdge[ iEdge ]; - if ( points._params.size() >= iSeg+1 ) + if ( points._params.size() <= iSeg+1 ) + return false; + + return Abs( points._params[ iSeg ] - points._params[ iSeg+1 ]) < 1e-20; +} + +//================================================================================ +/*! + * \brief Moves (changes _param) a given BoundaryPoint to a closest EDGE end + */ +//================================================================================ + +bool SMESH_MAT2d::Boundary::moveToClosestEdgeEnd( BoundaryPoint& bp ) const +{ + if ( bp._edgeIndex >= _pointsPerEdge.size() ) return false; - return Abs( points._params[ iEdge ] - points._params[ iEdge+1 ]) < 1e-20; + const BndPoints& points = _pointsPerEdge[ bp._edgeIndex ]; + if ( bp._param - points._params[0] < points._params.back() - bp._param ) + bp._param = points._params[0]; + else + bp._param = points._params.back(); + + return true; } //================================================================================ @@ -1179,7 +1434,7 @@ void SMESH_MAT2d::Branch::init( vector& maEdges, //================================================================================ /*! - * \brief fill BranchEnd::_branches of its ends + * \brief fills BranchEnd::_branches of its ends */ //================================================================================ @@ -1197,6 +1452,47 @@ void SMESH_MAT2d::Branch::setBranchesToEnds( const vector< Branch >& branches ) } } +//================================================================================ +/*! + * \brief returns a BranchPoint corresponding to a TVDVertex + */ +//================================================================================ + +SMESH_MAT2d::BranchPoint SMESH_MAT2d::Branch::getPoint( const TVDVertex* vertex ) const +{ + BranchPoint p; + p._branch = this; + p._iEdge = 0; + + if ( vertex == _maEdges[0]->vertex1() ) + { + p._edgeParam = 0; + } + else + { + for ( ; p._iEdge < _maEdges.size(); ++p._iEdge ) + if ( vertex == _maEdges[ p._iEdge ]->vertex0() ) + { + p._edgeParam = _params[ p._iEdge ]; + break; + } + } + return p; +} + +//================================================================================ +/*! + * \brief Sets a proxy point for a removed branch + * \param [in] proxyPoint - a point of another branch to which all points of this + * branch are mapped + */ +//================================================================================ + +void SMESH_MAT2d::Branch::setRemoved( const BranchPoint& proxyPoint ) +{ + _proxyPoint = proxyPoint; +} + //================================================================================ /*! * \brief Returns points on two EDGEs, equidistant from a given point of this Branch @@ -1213,7 +1509,7 @@ bool SMESH_MAT2d::Branch::getBoundaryPoints(double param, { if ( param < _params[0] || param > _params.back() ) return false; - + // look for an index of a MA edge by param double ip = param * _params.size(); size_t i = size_t( Min( int( _maEdges.size()-1), int( ip ))); @@ -1242,8 +1538,13 @@ bool SMESH_MAT2d::Branch::getBoundaryPoints(std::size_t iMAEdge, BoundaryPoint& bp1, BoundaryPoint& bp2 ) const { + if ( isRemoved() ) + return _proxyPoint._branch->getBoundaryPoints( _proxyPoint, bp1, bp2 ); + if ( iMAEdge > _maEdges.size() ) return false; + if ( iMAEdge == _maEdges.size() ) + iMAEdge = _maEdges.size() - 1; size_t iGeom1 = getGeomEdge( _maEdges[ iMAEdge ] ); size_t iGeom2 = getGeomEdge( _maEdges[ iMAEdge ]->twin() ); @@ -1275,8 +1576,16 @@ bool SMESH_MAT2d::Branch::getBoundaryPoints(const BranchPoint& p, bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const { + if ( this != p._branch && p._branch ) + return p._branch->getParameter( p, u ); + + if ( isRemoved() ) + return _proxyPoint._branch->getParameter( _proxyPoint, u ); + if ( p._iEdge > _params.size()-1 ) return false; + if ( p._iEdge == _params.size()-1 ) + return u = 1.; u = ( _params[ p._iEdge ] * ( 1 - p._edgeParam ) + _params[ p._iEdge+1 ] * p._edgeParam ); @@ -1366,8 +1675,8 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >& size_t iSeg1 = getBndSegment( maEdges[ i-1 ] ); size_t iSeg2 = getBndSegment( maEdges[ i ] ); - bool isConcaPrev = _boundary->IsConcaveSegment( edgeIDs1.back(), iSeg1 ); - bool isConcaNext = _boundary->IsConcaveSegment( ie1, iSeg2 ); + bool isConcaPrev = _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 ); + bool isConcaNext = _boundary->isConcaveSegment( ie1, iSeg2 ); if ( !isConcaNext && !isConcaPrev ) return false; @@ -1381,7 +1690,7 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >& while ( iNext < maEdges.size() ) { iSeg2 = getBndSegment( maEdges[ iNext ] ); - if ( _boundary->IsConcaveSegment( ie1, iSeg2 )) + if ( _boundary->isConcaveSegment( ie1, iSeg2 )) ++iNext; else break; @@ -1403,7 +1712,7 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >& } if ( vertexFound ) { - i = --iNext; + iPrev = i = --iNext; // not to add a BP in the moddle isConcaveV = true; } } @@ -1413,7 +1722,7 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >& while ( iPrev-1 >= 0 ) { iSeg1 = getBndSegment( maEdges[ iPrev-1 ] ); - if ( _boundary->IsConcaveSegment( edgeIDs1.back(), iSeg1 )) + if ( _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 )) --iPrev; else break; @@ -1423,7 +1732,7 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >& if ( iPrev < i-1 || iNext > i ) { // no VERTEX on the opposite EDGE, put the Branch Point in the middle - double par1 = _params[ iPrev ], par2 = _params[ iNext ]; + double par1 = _params[ iPrev+1 ], par2 = _params[ iNext ]; double midPar = 0.5 * ( par1 + par2 ); divisionPnt._iEdge = iPrev; while ( _params[ divisionPnt._iEdge + 1 ] < midPar ) diff --git a/src/SMESHUtils/SMESH_MAT2d.hxx b/src/SMESHUtils/SMESH_MAT2d.hxx index 2b60605bc..2ba0066c5 100644 --- a/src/SMESHUtils/SMESH_MAT2d.hxx +++ b/src/SMESHUtils/SMESH_MAT2d.hxx @@ -88,6 +88,7 @@ namespace SMESH_MAT2d /*! * \brief Branch is a set of MA edges enclosed between branch points and/or MA ends. * It's main feature is to return two BoundaryPoint's per a point on it. + * Points on a Branch are defined by [0,1] parameter */ class SMESHUtils_EXPORT Branch { @@ -114,15 +115,20 @@ namespace SMESH_MAT2d std::vector< std::size_t >& edgeIDs2, std::vector< BranchPoint >& divPoints) const; - // construction + bool isRemoved() const { return _proxyPoint._branch; } + + public: // internal: construction + void init( std::vector& maEdges, const Boundary* boundary, std::map< const TVDVertex*, BranchEndType > endType); void setBranchesToEnds( const std::vector< Branch >& branches); + BranchPoint getPoint( const TVDVertex* vertex ) const; + void setRemoved( const BranchPoint& proxyPoint ); - static void setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge ); - static std::size_t getGeomEdge( const TVDEdge* maEdge ); - static void setBndSegment( std::size_t segIndex, const TVDEdge* maEdge ); + static void setGeomEdge ( std::size_t geomIndex, const TVDEdge* maEdge ); + static std::size_t getGeomEdge ( const TVDEdge* maEdge ); + static void setBndSegment( std::size_t segIndex, const TVDEdge* maEdge ); static std::size_t getBndSegment( const TVDEdge* maEdge ); private: @@ -142,6 +148,7 @@ namespace SMESH_MAT2d const Boundary* _boundary; // face boundary BranchEnd _endPoint1; BranchEnd _endPoint2; + BranchPoint _proxyPoint; }; //------------------------------------------------------------------------------------- @@ -173,7 +180,9 @@ namespace SMESH_MAT2d bool getBranchPoint( const std::size_t iEdge, double u, BranchPoint& p ) const; - bool IsConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const; + bool isConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const; + + bool moveToClosestEdgeEnd( BoundaryPoint& bp ) const; private: std::vector< BndPoints > _pointsPerEdge; @@ -201,11 +210,12 @@ namespace SMESH_MAT2d const std::vector< TopoDS_Edge >& edges, const double minSegLen, const bool ignoreCorners = false ); - const Boundary& getBoundary() const { return _boundary; } - const std::vector< Branch >& getBranches() const { return _branch; } + std::size_t nbBranches() const { return _nbBranches; } + const Branch* getBranch(size_t i) const; const std::vector< const BranchEnd* >& getBranchPoints() const { return _branchPnt; } + const Boundary& getBoundary() const { return _boundary; } - void getPoints( const Branch& branch, std::vector< gp_XY >& points) const; + void getPoints( const Branch* branch, std::vector< gp_XY >& points) const; Adaptor3d_Curve* make3DCurve(const Branch& branch) const; private: @@ -214,6 +224,7 @@ namespace SMESH_MAT2d TopoDS_Face _face; TVD _vd; std::vector< Branch > _branch; + std::size_t _nbBranches; // removed branches ignored std::vector< const BranchEnd* > _branchPnt; Boundary _boundary; double _scale[2]; diff --git a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx index a6d372d5f..665a84efa 100644 --- a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx +++ b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx @@ -542,11 +542,11 @@ namespace TopoDS_Edge makeEdgeFromMA( SMESH_MesherHelper& theHelper, const SMESH_MAT2d::MedialAxis& theMA ) { - if ( theMA.getBranches().size() != 1 ) + if ( theMA.nbBranches() != 1 ) return TopoDS_Edge(); vector< gp_XY > uv; - theMA.getPoints( theMA.getBranches()[0], uv ); + theMA.getPoints( theMA.getBranch(0), uv ); if ( uv.size() < 2 ) return TopoDS_Edge(); @@ -775,7 +775,7 @@ namespace double uMA; SMESH_MAT2d::BoundaryPoint bp[2]; - const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0]; + const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); for ( size_t i = 0; i < theDivPoints.size(); ++i ) { @@ -830,7 +830,7 @@ namespace const vector& theSinuEdges, const size_t theSinuSide0Size) { - if ( theMA.getBranches().size() != 1 ) + if ( theMA.nbBranches() != 1 ) return false; // normalize theMAParams @@ -858,7 +858,7 @@ namespace // hasComputed = true; } - const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0]; + const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); SMESH_MAT2d::BoundaryPoint bp[2]; vector< std::size_t > edgeIDs1, edgeIDs2; @@ -1248,4 +1248,3 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::Evaluate(SMESH_Mesh & theMesh, { return StdMeshers_Quadrangle_2D::Evaluate(theMesh,theShape,theResMap); } -