X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_QuadFromMedialAxis_1D2D.cxx;h=0b9005bcaa07fadefbd9c679a0774c0900a7bf3c;hp=4ac71c0b818a8b09e31b0d0d8d7a618978c24738;hb=2f529dcd2629679dadcca3047583bfcf28ca7b1a;hpb=23d90107acec5e54ded86d9f309fe5cb42720b78 diff --git a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx index 4ac71c0b8..0b9005bca 100644 --- a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx +++ b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE +// 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 @@ -25,6 +25,7 @@ #include "StdMeshers_QuadFromMedialAxis_1D2D.hxx" +#include "SMESHDS_Mesh.hxx" #include "SMESH_Block.hxx" #include "SMESH_Gen.hxx" #include "SMESH_MAT2d.hxx" @@ -33,10 +34,14 @@ #include "SMESH_MesherHelper.hxx" #include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" +#include "SMESH_subMeshEventListener.hxx" #include "StdMeshers_FaceSide.hxx" +#include "StdMeshers_LayerDistribution.hxx" +#include "StdMeshers_NumberOfLayers.hxx" #include "StdMeshers_Regular_1D.hxx" #include "StdMeshers_ViscousLayers2D.hxx" +#include #include #include #include @@ -45,6 +50,7 @@ #include #include #include +#include #include #include #include @@ -56,6 +62,8 @@ #include #include +using namespace std; + //================================================================================ /*! * \brief 1D algo @@ -63,16 +71,78 @@ class StdMeshers_QuadFromMedialAxis_1D2D::Algo1D : public StdMeshers_Regular_1D { public: - Algo1D(int studyId, SMESH_Gen* gen): - StdMeshers_Regular_1D( gen->GetANewId(), studyId, gen ) + Algo1D(SMESH_Gen* gen): + StdMeshers_Regular_1D( gen->GetANewId(), gen ) { } void SetSegmentLength( double len ) { + SMESH_Algo::_usedHypList.clear(); _value[ BEG_LENGTH_IND ] = len; - _value[ PRECISION_IND ] = 1e-7; + _value[ PRECISION_IND ] = 1e-7; _hypType = LOCAL_LENGTH; } + void SetRadialDistribution( const SMESHDS_Hypothesis* hyp ) + { + SMESH_Algo::_usedHypList.clear(); + if ( !hyp ) + return; + + if ( const StdMeshers_NumberOfLayers* nl = + dynamic_cast< const StdMeshers_NumberOfLayers* >( hyp )) + { + _ivalue[ NB_SEGMENTS_IND ] = nl->GetNumberOfLayers(); + _ivalue[ DISTR_TYPE_IND ] = 0; + _hypType = NB_SEGMENTS; + } + if ( const StdMeshers_LayerDistribution* ld = + dynamic_cast< const StdMeshers_LayerDistribution* >( hyp )) + { + if ( SMESH_Hypothesis* h = ld->GetLayerDistribution() ) + { + SMESH_Algo::_usedHypList.clear(); + SMESH_Algo::_usedHypList.push_back( h ); + } + } + } + void ComputeDistribution(SMESH_MesherHelper& theHelper, + const gp_Pnt& thePnt1, + const gp_Pnt& thePnt2, + list< double >& theParams) + { + SMESH_Mesh& mesh = *theHelper.GetMesh(); + TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( thePnt1, thePnt2 ); + + SMESH_Hypothesis::Hypothesis_Status aStatus; + CheckHypothesis( mesh, edge, aStatus ); + + theParams.clear(); + BRepAdaptor_Curve C3D(edge); + double f = C3D.FirstParameter(), l = C3D.LastParameter(), len = thePnt1.Distance( thePnt2 ); + if ( !StdMeshers_Regular_1D::computeInternalParameters( mesh, C3D, len, f, l, theParams, false)) + { + for ( size_t i = 1; i < 15; ++i ) + theParams.push_back( i/15. ); // ???? + } + else + { + for (list::iterator itU = theParams.begin(); itU != theParams.end(); ++itU ) + *itU /= len; + } + } + virtual const list & + GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool) + { + return SMESH_Algo::_usedHypList; + } + virtual bool CheckHypothesis(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + SMESH_Hypothesis::Hypothesis_Status& aStatus) + { + if ( !SMESH_Algo::_usedHypList.empty() ) + return StdMeshers_Regular_1D::CheckHypothesis( aMesh, aShape, aStatus ); + return true; + } }; //================================================================================ @@ -82,20 +152,21 @@ public: //================================================================================ StdMeshers_QuadFromMedialAxis_1D2D::StdMeshers_QuadFromMedialAxis_1D2D(int hypId, - int studyId, SMESH_Gen* gen) - : StdMeshers_Quadrangle_2D(hypId, studyId, gen), + : StdMeshers_Quadrangle_2D(hypId, gen), _regular1D( 0 ) { _name = "QuadFromMedialAxis_1D2D"; _shapeType = (1 << TopAbs_FACE); _onlyUnaryInput = true; // FACE by FACE so far _requireDiscreteBoundary = false; // make 1D by myself - _supportSubmeshes = true; // make 1D by myself + _supportSubmeshes = true; // make 1D by myself _neededLowerHyps[ 1 ] = true; // suppress warning on hiding a global 1D algo _neededLowerHyps[ 2 ] = true; // suppress warning on hiding a global 2D algo _compatibleHypothesis.clear(); - //_compatibleHypothesis.push_back("ViscousLayers2D"); + _compatibleHypothesis.push_back("ViscousLayers2D"); + _compatibleHypothesis.push_back("LayerDistribution2D"); + _compatibleHypothesis.push_back("NumberOfLayers2D"); } //================================================================================ @@ -120,6 +191,12 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::CheckHypothesis(SMESH_Mesh& aMe const TopoDS_Shape& aShape, Hypothesis_Status& aStatus) { + aStatus = HYP_OK; + + // get one main optional hypothesis + const list & hyps = GetUsedHypothesis(aMesh, aShape); + _hyp2D = hyps.empty() ? 0 : hyps.front(); + return true; // does not require hypothesis } @@ -133,13 +210,14 @@ namespace */ struct SinuousFace { - FaceQuadStruct::Ptr _quad; - vector< TopoDS_Edge > _edges; - vector< TopoDS_Edge > _sinuSide[2], _shortSide[2]; - vector< TopoDS_Edge > _sinuEdges; - int _nbWires; - list< int > _nbEdgesInWire; - TMergeMap _nodesToMerge; + FaceQuadStruct::Ptr _quad; + vector< TopoDS_Edge > _edges; + vector< TopoDS_Edge > _sinuSide[2], _shortSide[2]; + vector< TopoDS_Edge > _sinuEdges; + vector< Handle(Geom_Curve) > _sinuCurves; + int _nbWires; + list< int > _nbEdgesInWire; + TMergeMap _nodesToMerge; SinuousFace( const TopoDS_Face& f ): _quad( new FaceQuadStruct ) { @@ -151,6 +229,7 @@ namespace _quad->face = f; } const TopoDS_Face& Face() const { return _quad->face; } + bool IsRing() const { return _shortSide[0].empty() && !_sinuSide[0].empty(); } }; //================================================================================ @@ -165,6 +244,43 @@ namespace } }; + //================================================================================ + /*! + * \brief Event listener which removes mesh from EDGEs when 2D hyps change + */ + struct EdgeCleaner : public SMESH_subMeshEventListener + { + int _prevAlgoEvent; + EdgeCleaner(): + SMESH_subMeshEventListener( /*isDeletable=*/true, + "StdMeshers_QuadFromMedialAxis_1D2D::EdgeCleaner") + { + _prevAlgoEvent = -1; + } + virtual void ProcessEvent(const int event, + const int eventType, + SMESH_subMesh* faceSubMesh, + SMESH_subMeshEventListenerData* data, + const SMESH_Hypothesis* hyp) + { + if ( eventType == SMESH_subMesh::ALGO_EVENT ) + { + _prevAlgoEvent = event; + return; + } + // SMESH_subMesh::COMPUTE_EVENT + if ( _prevAlgoEvent == SMESH_subMesh::REMOVE_HYP || + _prevAlgoEvent == SMESH_subMesh::REMOVE_ALGO || + _prevAlgoEvent == SMESH_subMesh::MODIF_HYP ) + { + SMESH_subMeshIteratorPtr smIt = faceSubMesh->getDependsOnIterator(/*includeSelf=*/false); + while ( smIt->more() ) + smIt->next()->ComputeStateEngine( SMESH_subMesh::CLEAN ); + } + _prevAlgoEvent = -1; + } + }; + //================================================================================ /*! * \brief Return a member of a std::pair @@ -235,8 +351,8 @@ namespace algos[i] = sm->GetAlgo(); } - const int nbSegDflt = mesh->GetGen()->GetDefaultNbSegments(); - double minSegLen = Precision::Infinite(); + int nbSegDflt = mesh->GetGen() ? mesh->GetGen()->GetDefaultNbSegments() : 15; + double minSegLen = Precision::Infinite(); for ( size_t i = 0; i < theEdges.size(); ++i ) { @@ -265,7 +381,8 @@ namespace tmpMesh.ShapeToMesh( TopoDS_Shape()); tmpMesh.ShapeToMesh( theEdges[i] ); try { - mesh->GetGen()->Compute( tmpMesh, theEdges[i], true, true ); // make nodes on VERTEXes + if ( !mesh->GetGen() ) continue; // tmp mesh + mesh->GetGen()->Compute( tmpMesh, theEdges[i], SMESH_Gen::SHAPE_ONLY_UPWARD ); // make nodes on VERTEXes if ( !algo->Compute( tmpMesh, theEdges[i] )) continue; } @@ -400,6 +517,7 @@ namespace size_t nbOutEdges = theSinuFace._nbEdgesInWire.front(); theSinuEdges[0].assign ( allEdges.begin(), allEdges.begin() + nbOutEdges ); theSinuEdges[1].assign ( allEdges.begin() + nbOutEdges, allEdges.end() ); + theSinuFace._sinuEdges = allEdges; return true; } if ( theSinuFace._nbWires > 2 ) @@ -435,8 +553,8 @@ namespace allEdges, theShortEdges[ nbBranchPoints > 0 ] )) return false; - for ( size_t iS = 0; iS < theShortEdges[ nbBranchPoints ].size(); ++iS ) - shortMap.Add( theShortEdges[ nbBranchPoints ][ iS ]); + for ( size_t iS = 0; iS < theShortEdges[ nbBranchPoints > 0 ].size(); ++iS ) + shortMap.Add( theShortEdges[ nbBranchPoints > 0 ][ iS ]); ++nbBranchPoints; } @@ -485,7 +603,7 @@ namespace theSinuEdges [0].size() > 0 && theSinuEdges [1].size() > 0 ); // the sinuous EDGEs can be composite and C0 continuous, - // therefor we use a complex criterion to find TWO short non-sinuous EDGEs + // therefore we use a complex criterion to find TWO short non-sinuous EDGEs // and the rest EDGEs will be treated as sinuous. // A short edge should have the following features: // a) straight @@ -591,27 +709,45 @@ namespace //================================================================================ TopoDS_Edge makeEdgeFromMA( SMESH_MesherHelper& theHelper, - const SMESH_MAT2d::MedialAxis& theMA ) + const SMESH_MAT2d::MedialAxis& theMA, + const double theMinSegLen) { - 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(); TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() ); Handle(Geom_Surface) surface = BRep_Tool::Surface( face ); - // cout << "from salome.geom import geomBuilder" << endl; - // cout << "geompy = geomBuilder.New(salome.myStudy)" << endl; - Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, uv.size()); - for ( size_t i = 0; i < uv.size(); ++i ) + vector< gp_Pnt > pnt; + pnt.reserve( uv.size() * 2 ); + pnt.push_back( surface->Value( uv[0].X(), uv[0].Y() )); + for ( size_t i = 1; i < uv.size(); ++i ) { gp_Pnt p = surface->Value( uv[i].X(), uv[i].Y() ); + int nbDiv = int( p.Distance( pnt.back() ) / theMinSegLen ); + for ( int iD = 1; iD < nbDiv; ++iD ) + { + double R = iD / double( nbDiv ); + gp_XY uvR = uv[i-1] * (1 - R) + uv[i] * R; + pnt.push_back( surface->Value( uvR.X(), uvR.Y() )); + } + pnt.push_back( p ); + } + + // cout << "from salome.geom import geomBuilder" << endl; + // cout << "geompy = geomBuilder.New()" << endl; + Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, pnt.size()); + for ( size_t i = 0; i < pnt.size(); ++i ) + { + gp_Pnt& p = pnt[i]; points->SetValue( i+1, p ); - //cout << "geompy.MakeVertex( "<< p.X()<<", " << p.Y()<<", " << p.Z()<<" )" << endl; + // cout << "geompy.MakeVertex( "<< p.X()<<", " << p.Y()<<", " << p.Z() + // <<" theName = 'p_" << i << "')" << endl; } GeomAPI_Interpolate interpol( points, /*isClosed=*/false, gp::Resolution()); @@ -657,23 +793,45 @@ namespace const SMESH_MAT2d::MedialAxis& theMA, const SinuousFace& theSinuFace, SMESH_Algo* the1dAlgo, + const double theMinSegLen, vector& theMAParams ) { - // check if all EDGEs of one size are meshed, then MA discretization is not needed + // Check if all EDGEs of one size are meshed, then MA discretization is not needed SMESH_Mesh* mesh = theHelper.GetMesh(); size_t nbComputedEdges[2] = { 0, 0 }; for ( size_t iS = 0; iS < 2; ++iS ) for ( size_t i = 0; i < theSinuFace._sinuSide[iS].size(); ++i ) { - bool isComputed = ( ! mesh->GetSubMesh( theSinuFace._sinuSide[iS][i] )->IsEmpty() ); + const TopoDS_Edge& sinuEdge = theSinuFace._sinuSide[iS][i]; + SMESH_subMesh* sm = mesh->GetSubMesh( sinuEdge ); + bool isComputed = ( !sm->IsEmpty() ); + if ( isComputed ) + { + TopAbs_ShapeEnum shape = getHypShape( mesh, sinuEdge ); + if ( shape == TopAbs_SHAPE || shape <= TopAbs_FACE ) + { + // EDGE computed using global hypothesis -> clear it + bool hasComputedFace = false; + PShapeIteratorPtr faceIt = theHelper.GetAncestors( sinuEdge, *mesh, TopAbs_FACE ); + while ( const TopoDS_Shape* face = faceIt->next() ) + if (( !face->IsSame( theSinuFace.Face() )) && + ( hasComputedFace = !mesh->GetSubMesh( *face )->IsEmpty() )) + break; + if ( !hasComputedFace ) + { + sm->ComputeStateEngine( SMESH_subMesh::CLEAN ); + isComputed = false; + } + } + } nbComputedEdges[ iS ] += isComputed; } if ( nbComputedEdges[0] == theSinuFace._sinuSide[0].size() || nbComputedEdges[1] == theSinuFace._sinuSide[1].size() ) return true; // discretization is not needed - - TopoDS_Edge branchEdge = makeEdgeFromMA( theHelper, theMA ); + // Make MA EDGE + TopoDS_Edge branchEdge = makeEdgeFromMA( theHelper, theMA, theMinSegLen ); if ( branchEdge.IsNull() ) return false; @@ -681,18 +839,24 @@ namespace // BRepTools::Write( branchEdge, file); // cout << "Write " << file << endl; - // look for a most local hyps assigned to theSinuEdges - TopoDS_Edge edge = theSinuFace._sinuEdges[0]; - int mostSimpleShape = (int) getHypShape( mesh, edge ); - for ( size_t i = 1; i < theSinuFace._sinuEdges.size(); ++i ) + + // Find 1D algo to mesh branchEdge + + // look for a most local 1D hyp assigned to the FACE + int mostSimpleShape = -1, maxShape = TopAbs_EDGE; + TopoDS_Edge edge; + for ( size_t i = 0; i < theSinuFace._sinuEdges.size(); ++i ) { - int shapeType = (int) getHypShape( mesh, theSinuFace._sinuEdges[i] ); - if ( shapeType > mostSimpleShape ) + TopAbs_ShapeEnum shapeType = getHypShape( mesh, theSinuFace._sinuEdges[i] ); + if ( mostSimpleShape < shapeType && shapeType < maxShape ) + { edge = theSinuFace._sinuEdges[i]; + mostSimpleShape = shapeType; + } } SMESH_Algo* algo = the1dAlgo; - if ( mostSimpleShape != TopAbs_SHAPE ) + if ( mostSimpleShape > -1 ) { algo = mesh->GetSubMesh( edge )->GetAlgo(); SMESH_Hypothesis::Hypothesis_Status status; @@ -703,7 +867,7 @@ namespace TmpMesh tmpMesh; tmpMesh.ShapeToMesh( branchEdge ); try { - mesh->GetGen()->Compute( tmpMesh, branchEdge, true, true ); // make nodes on VERTEXes + mesh->GetGen()->Compute( tmpMesh, branchEdge, SMESH_Gen::SHAPE_ONLY_UPWARD ); // make nodes on VERTEXes if ( !algo->Compute( tmpMesh, branchEdge )) return false; } @@ -787,7 +951,7 @@ namespace { const SMDS_MeshNode* _node; double _u; - int _edgeInd; // index in theSinuEdges vector + size_t _edgeInd; // index in theSinuEdges vector NodePoint(): _node(0), _u(0), _edgeInd(-1) {} NodePoint(const SMDS_MeshNode* n, double u, size_t iEdge ): _node(n), _u(u), _edgeInd(iEdge) {} @@ -795,9 +959,10 @@ namespace NodePoint(const SMESH_MAT2d::BoundaryPoint& p) : _node(0), _u(p._param), _edgeInd(p._edgeIndex) {} gp_Pnt Point(const vector< Handle(Geom_Curve) >& curves) const { - return curves[ _edgeInd ]->Value( _u ); + return _node ? SMESH_TNodeXYZ(_node) : curves[ _edgeInd ]->Value( _u ); } }; + typedef multimap< double, pair< NodePoint, NodePoint > > TMAPar2NPoints; //================================================================================ /*! @@ -810,9 +975,11 @@ namespace */ //================================================================================ - bool findVertex( NodePoint& theNodePnt, - const vector& theSinuEdges, - SMESHDS_Mesh* theMeshDS) + bool findVertexAndNode( NodePoint& theNodePnt, + const vector& theSinuEdges, + SMESHDS_Mesh* theMeshDS = 0, + size_t theEdgeIndPrev = 0, + size_t theEdgeIndNext = 0) { if ( theNodePnt._edgeInd >= theSinuEdges.size() ) return false; @@ -826,8 +993,10 @@ namespace V = SMESH_MesherHelper::IthVertex( 0, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false); else if ( Abs( l - theNodePnt._u ) < tol ) V = SMESH_MesherHelper::IthVertex( 1, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false); + else if ( theEdgeIndPrev != theEdgeIndNext ) + TopExp::CommonVertex( theSinuEdges[theEdgeIndPrev], theSinuEdges[theEdgeIndNext], V ); - if ( !V.IsNull() ) + if ( !V.IsNull() && theMeshDS ) { theNodePnt._node = SMESH_Algo::VertexNode( V, theMeshDS ); if ( !theNodePnt._node ) @@ -836,9 +1005,8 @@ namespace theNodePnt._node = theMeshDS->AddNode( p.X(), p.Y(), p.Z() ); theMeshDS->SetNodeOnVertex( theNodePnt._node, V ); } - return true; } - return false; + return !V.IsNull(); } //================================================================================ @@ -850,31 +1018,75 @@ namespace * \param [in] theDivPoints - projections of VERTEXes to MA * \param [in] theSinuEdges - the sinuous EDGEs * \param [in] theSideEdgeIDs - indices of sinuous EDGEs per side - * \param [in] theIsEdgeComputed - is sinuous EGDE is meshed + * \param [in] theIsEdgeComputed - is sinuous EDGE is meshed * \param [in,out] thePointsOnE - the map to fill * \param [out] theNodes2Merge - the map of nodes to merge */ //================================================================================ - bool projectVertices( SMESH_MesherHelper& theHelper, - //const double theMinSegLen, - const SMESH_MAT2d::MedialAxis& theMA, - const vector< SMESH_MAT2d::BranchPoint >& theDivPoints, - const vector& theSinuEdges, - const vector< Handle(Geom_Curve) >& theCurves, - const vector< bool >& theIsEdgeComputed, - map< double, pair< NodePoint, NodePoint > > & thePointsOnE, - TMergeMap& theNodes2Merge) + bool projectVertices( SMESH_MesherHelper& theHelper, + const SMESH_MAT2d::MedialAxis& theMA, + vector< SMESH_MAT2d::BranchPoint >& theDivPoints, + const vector< std::size_t > & theEdgeIDs1, + const vector< std::size_t > & theEdgeIDs2, + const vector< bool >& theIsEdgeComputed, + TMAPar2NPoints & thePointsOnE, + SinuousFace& theSinuFace) { if ( theDivPoints.empty() ) return true; SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); + const vector< TopoDS_Edge >& theSinuEdges = theSinuFace._sinuEdges; + const vector< Handle(Geom_Curve) >& theCurves = theSinuFace._sinuCurves; double uMA; - SMESH_MAT2d::BoundaryPoint bp[2]; - const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0]; + SMESH_MAT2d::BoundaryPoint bp[2]; // 2 sinuous sides + const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); + { + // add to thePointsOnE NodePoint's of ends of theSinuEdges + if ( !branch.getBoundaryPoints( 0., bp[0], bp[1] ) || + !theMA.getBoundary().moveToClosestEdgeEnd( bp[0] )) return false; + if ( !theSinuFace.IsRing() && + !theMA.getBoundary().moveToClosestEdgeEnd( bp[1] )) return false; + NodePoint np0( bp[0] ), np1( bp[1] ); + findVertexAndNode( np0, theSinuEdges, meshDS ); + findVertexAndNode( np1, theSinuEdges, meshDS ); + thePointsOnE.insert( make_pair( -0.1, make_pair( np0, np1 ))); + } + if ( !theSinuFace.IsRing() ) + { + if ( !branch.getBoundaryPoints( 1., bp[0], bp[1] ) || + !theMA.getBoundary().moveToClosestEdgeEnd( bp[0] ) || + !theMA.getBoundary().moveToClosestEdgeEnd( bp[1] )) return false; + NodePoint np0( bp[0] ), np1( bp[1] ); + findVertexAndNode( np0, theSinuEdges, meshDS ); + findVertexAndNode( np1, theSinuEdges, meshDS ); + thePointsOnE.insert( make_pair( 1.1, make_pair( np0, np1))); + } + else + { + // project a VERTEX of outer sinuous side corresponding to branch(0.) + // which is not included into theDivPoints + if ( ! ( theDivPoints[0]._iEdge == 0 && + theDivPoints[0]._edgeParam == 0. )) // recursive call + { + SMESH_MAT2d::BranchPoint brp( &branch, 0, 0. ); + vector< SMESH_MAT2d::BranchPoint > divPoint( 1, brp ); + vector< std::size_t > edgeIDs1(2), edgeIDs2(2); + edgeIDs1[0] = theEdgeIDs1.back(); + edgeIDs1[1] = theEdgeIDs1[0]; + edgeIDs2[0] = theEdgeIDs2.back(); + edgeIDs2[1] = theEdgeIDs2[0]; + projectVertices( theHelper, theMA, divPoint, edgeIDs1, edgeIDs2, + theIsEdgeComputed, thePointsOnE, theSinuFace ); + } + } + // project theDivPoints and keep projections to merge + + TMAPar2NPoints::iterator u2NP; + vector< TMAPar2NPoints::iterator > projToMerge; for ( size_t i = 0; i < theDivPoints.size(); ++i ) { if ( !branch.getParameter( theDivPoints[i], uMA )) @@ -882,23 +1094,53 @@ namespace if ( !branch.getBoundaryPoints( theDivPoints[i], bp[0], bp[1] )) return false; - NodePoint np[2] = { NodePoint( bp[0] ), - NodePoint( bp[1] )}; - bool isVertex[2] = { findVertex( np[0], theSinuEdges, meshDS ), - findVertex( np[1], theSinuEdges, meshDS )}; + NodePoint np[2] = { + NodePoint( bp[0] ), + NodePoint( bp[1] ) + }; + bool isVertex[2] = { + findVertexAndNode( np[0], theSinuEdges, meshDS, theEdgeIDs1[i], theEdgeIDs1[i+1] ), + findVertexAndNode( np[1], theSinuEdges, meshDS, theEdgeIDs2[i], theEdgeIDs2[i+1] ) + }; + const size_t iVert = isVertex[0] ? 0 : 1; // side with a VERTEX + const size_t iNode = 1 - iVert; // opposite (meshed?) side + + if ( isVertex[0] != isVertex[1] ) // try to find an opposite VERTEX + { + theMA.getBoundary().moveToClosestEdgeEnd( bp[iNode] ); // EDGE -> VERTEX + SMESH_MAT2d::BranchPoint brp; + theMA.getBoundary().getBranchPoint( bp[iNode], brp ); // WIRE -> MA + SMESH_MAT2d::BoundaryPoint bp2[2]; + branch.getBoundaryPoints( brp, bp2[0], bp2[1] ); // MA -> WIRE + NodePoint np2[2] = { NodePoint( bp2[0]), NodePoint( bp2[1]) }; + findVertexAndNode( np2[0], theSinuEdges, meshDS ); + findVertexAndNode( np2[1], theSinuEdges, meshDS ); + if ( np2[ iVert ]._node == np[ iVert ]._node && + np2[ iNode ]._node) + { + np[ iNode ] = np2[ iNode ]; + isVertex[ iNode ] = true; + } + } - map< double, pair< NodePoint, NodePoint > >::iterator u2NP = - thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1]))).first; + u2NP = thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1]))); if ( !isVertex[0] && !isVertex[1] ) return false; // error if ( isVertex[0] && isVertex[1] ) continue; - const size_t iVertex = isVertex[0] ? 0 : 1; - const size_t iNode = 1 - iVertex; - bool isOppComputed = theIsEdgeComputed[ np[ iNode ]._edgeInd ]; - if ( !isOppComputed ) - continue; + // bool isOppComputed = theIsEdgeComputed[ np[ iNode ]._edgeInd ]; + // if ( isOppComputed ) + projToMerge.push_back( u2NP ); + } + + // merge projections + + for ( size_t i = 0; i < projToMerge.size(); ++i ) + { + u2NP = projToMerge[i]; + const size_t iVert = get( u2NP->second, 0 )._node ? 0 : 1; // side with a VERTEX + const size_t iNode = 1 - iVert; // opposite (meshed?) side // a VERTEX is projected on a meshed EDGE; there are two options: // 1) a projected point is joined with a closet node if a strip between this and neighbor @@ -908,10 +1150,12 @@ namespace // projection is set to the BoundaryPoint of this projection // evaluate distance to neighbor projections - const double rShort = 0.2; - bool isShortPrev[2], isShortNext[2]; - map< double, pair< NodePoint, NodePoint > >::iterator u2NPPrev = u2NP, u2NPNext = u2NP; + const double rShort = 0.33; + bool isShortPrev[2], isShortNext[2], isPrevCloser[2]; + TMAPar2NPoints::iterator u2NPPrev = u2NP, u2NPNext = u2NP; --u2NPPrev; ++u2NPNext; + if ( u2NPNext == thePointsOnE.end() ) + u2NPNext = thePointsOnE.begin(); // hope theSinuFace.IsRing() for ( int iS = 0; iS < 2; ++iS ) // side with Vertex and side with Nodes { NodePoint np = get( u2NP->second, iS ); @@ -923,58 +1167,490 @@ namespace double distPrev = p.Distance( pPrev ); double distNext = p.Distance( pNext ); double r = distPrev / ( distPrev + distNext ); - isShortPrev[iS] = ( r < rShort ); - isShortNext[iS] = (( 1 - r ) > ( 1 - rShort )); + isShortPrev [iS] = ( r < rShort ); + isShortNext [iS] = (( 1 - r ) < rShort ); + isPrevCloser[iS] = (( r < 0.5 ) && ( theSinuFace.IsRing() || u2NPPrev->first > 0 )); } - map< double, pair< NodePoint, NodePoint > >::iterator u2NPClose; + TMAPar2NPoints::iterator u2NPClose; if (( isShortPrev[0] && isShortPrev[1] ) || // option 2) -> remove a too close projection ( isShortNext[0] && isShortNext[1] )) { - u2NPClose = isShortPrev[0] ? u2NPPrev : u2NPNext; - NodePoint& npProj = get( u2NP->second, iNode ); // NP of VERTEX projection - NodePoint npCloseN = get( u2NPClose->second, iNode); // NP close to npProj - NodePoint npCloseV = get( u2NPClose->second, iVertex); // NP close to VERTEX - if ( !npCloseV._node ) + u2NPClose = isPrevCloser[0] ? u2NPPrev : u2NPNext; + NodePoint& npProj = get( u2NP->second, iNode ); // NP of VERTEX projection + NodePoint& npVert = get( u2NP->second, iVert ); // NP of VERTEX + NodePoint npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj + NodePoint npCloseV = get( u2NPClose->second, iVert ); // NP close to npVert + if ( !npCloseV._node || npCloseV._node == npVert._node ) { npProj = npCloseN; - thePointsOnE.erase( isShortPrev[0] ? u2NPPrev : u2NPNext ); + thePointsOnE.erase( u2NPClose ); continue; } else { - // can't remove the neighbor projection as it is also from VERTEX, -> option 1) + // can't remove the neighbor projection as it is also from VERTEX -> option 1) } } - // else option 1) - wide enough -> "duplicate" existing node + // else: option 1) - wide enough -> "duplicate" existing node { - u2NPClose = isShortPrev[ iNode ] ? u2NPPrev : u2NPNext; - NodePoint& npProj = get( u2NP->second, iNode ); // NP of VERTEX projection - NodePoint& npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj - // npProj._edgeInd = npCloseN._edgeInd; + u2NPClose = isPrevCloser[ iNode ] ? u2NPPrev : u2NPNext; + NodePoint& npProj = get( u2NP->second, iNode ); // NP of VERTEX projection + NodePoint& npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj + npProj = npCloseN; + npProj._node = 0; + //npProj._edgeInd = npCloseN._edgeInd; // npProj._u = npCloseN._u + 1e-3 * Abs( get( u2NPPrev->second, iNode )._u - // get( u2NPNext->second, iNode )._u ); - gp_Pnt p = npProj.Point( theCurves ); - npProj._node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); - meshDS->SetNodeOnEdge( npProj._node, theSinuEdges[ npProj._edgeInd ], npProj._u ); + // gp_Pnt p = npProj.Point( theCurves ); + // npProj._node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); + // meshDS->SetNodeOnEdge( npProj._node, theSinuEdges[ npProj._edgeInd ], npProj._u ); - theNodes2Merge[ npCloseN._node ].push_back( npProj._node ); + //theNodes2Merge[ npCloseN._node ].push_back( npProj._node ); } } + + // remove auxiliary NodePoint's of ends of theSinuEdges + for ( u2NP = thePointsOnE.begin(); u2NP->first < 0; ) + thePointsOnE.erase( u2NP++ ); + thePointsOnE.erase( 1.1 ); + return true; } + double getUOnEdgeByPoint( const size_t iEdge, + const NodePoint* point, + SinuousFace& sinuFace ) + { + if ( point->_edgeInd == iEdge ) + return point->_u; + + TopoDS_Vertex V0 = TopExp::FirstVertex( sinuFace._sinuEdges[ iEdge ]); + TopoDS_Vertex V1 = TopExp::LastVertex ( sinuFace._sinuEdges[ iEdge ]); + gp_Pnt p0 = BRep_Tool::Pnt( V0 ); + gp_Pnt p1 = BRep_Tool::Pnt( V1 ); + gp_Pnt p = point->Point( sinuFace._sinuCurves ); + + double f,l; + BRep_Tool::Range( sinuFace._sinuEdges[ iEdge ], f,l ); + return p.SquareDistance( p0 ) < p.SquareDistance( p1 ) ? f : l; + } + + //================================================================================ + /*! + * \brief Move coincident nodes to make node params on EDGE unique + * \param [in] theHelper - the helper + * \param [in] thePointsOnE - nodes on two opposite river sides + * \param [in] theSinuFace - the sinuous FACE + * \param [out] theNodes2Merge - the map of nodes to merge + */ + //================================================================================ + + void separateNodes( SMESH_MesherHelper& theHelper, + const SMESH_MAT2d::MedialAxis& theMA, + TMAPar2NPoints & thePointsOnE, + SinuousFace& theSinuFace, + const vector< bool >& theIsComputedEdge) + { + if ( thePointsOnE.size() < 2 ) + return; + + SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); + const vector& theSinuEdges = theSinuFace._sinuEdges; + const vector< Handle(Geom_Curve) >& curves = theSinuFace._sinuCurves; + + //SMESH_MAT2d::BoundaryPoint bp[2]; + //const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); + + typedef TMAPar2NPoints::iterator TIterator; + + for ( int iSide = 0; iSide < 2; ++iSide ) // loop on two sinuous sides + { + // get a tolerance to compare points + double tol = Precision::Confusion(); + for ( size_t i = 0; i < theSinuFace._sinuSide[ iSide ].size(); ++i ) + tol = Max( tol , BRep_Tool::Tolerance( theSinuFace._sinuSide[ iSide ][ i ])); + + // find coincident points + TIterator u2NP = thePointsOnE.begin(); + vector< TIterator > sameU2NP( 1, u2NP++ ); + while ( u2NP != thePointsOnE.end() ) + { + for ( ; u2NP != thePointsOnE.end(); ++u2NP ) + { + NodePoint& np1 = get( sameU2NP.back()->second, iSide ); + NodePoint& np2 = get( u2NP ->second, iSide ); + + if (( !np1._node || !np2._node ) && + ( np1.Point( curves ).SquareDistance( np2.Point( curves )) < tol*tol )) + { + sameU2NP.push_back( u2NP ); + } + else if ( sameU2NP.size() == 1 ) + { + sameU2NP[ 0 ] = u2NP; + } + else + { + break; + } + } + + if ( sameU2NP.size() > 1 ) + { + // find an existing node on VERTEX among sameU2NP and get underlying EDGEs + const SMDS_MeshNode* existingNode = 0; + set< size_t > edgeInds; + NodePoint* np; + for ( size_t i = 0; i < sameU2NP.size(); ++i ) + { + np = &get( sameU2NP[i]->second, iSide ); + if ( np->_node ) + if ( !existingNode || np->_node->GetPosition()->GetDim() == 0 ) + existingNode = np->_node; + edgeInds.insert( np->_edgeInd ); + } + list< const SMDS_MeshNode* >& mergeNodes = theSinuFace._nodesToMerge[ existingNode ]; + + TIterator u2NPprev = sameU2NP.front(); + TIterator u2NPnext = sameU2NP.back() ; + if ( u2NPprev->first < 0. ) ++u2NPprev; + if ( u2NPnext->first > 1. ) --u2NPnext; + + set< size_t >::iterator edgeID = edgeInds.begin(); + for ( ; edgeID != edgeInds.end(); ++edgeID ) + { + // get U range on iEdge within which the equal points will be distributed + double u0, u1; + np = &get( u2NPprev->second, iSide ); + u0 = getUOnEdgeByPoint( *edgeID, np, theSinuFace ); + + np = &get( u2NPnext->second, iSide ); + u1 = getUOnEdgeByPoint( *edgeID, np, theSinuFace ); + + if ( u0 == u1 ) + { + if ( u2NPprev != thePointsOnE.begin() ) --u2NPprev; + if ( u2NPnext != --thePointsOnE.end() ) ++u2NPnext; + np = &get( u2NPprev->second, iSide ); + u0 = getUOnEdgeByPoint( *edgeID, np, theSinuFace ); + np = &get( u2NPnext->second, iSide ); + u1 = getUOnEdgeByPoint( *edgeID, np, theSinuFace ); + } + + // distribute points and create nodes + double du = ( u1 - u0 ) / ( sameU2NP.size() + 1 /*!existingNode*/ ); + double u = u0 + du; + for ( size_t i = 0; i < sameU2NP.size(); ++i ) + { + np = &get( sameU2NP[i]->second, iSide ); + if ( !np->_node && *edgeID == np->_edgeInd ) + { + np->_u = u; + u += du; + gp_Pnt p = np->Point( curves ); + np->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); + meshDS->SetNodeOnEdge( np->_node, theSinuEdges[ *edgeID ], np->_u ); + + if ( theIsComputedEdge[ *edgeID ]) + mergeNodes.push_back( np->_node ); + } + } + } + + sameU2NP.resize( 1 ); + u2NP = ++sameU2NP.back(); + sameU2NP[ 0 ] = u2NP; + + } // if ( sameU2NP.size() > 1 ) + } // while ( u2NP != thePointsOnE.end() ) + } // for ( int iSide = 0; iSide < 2; ++iSide ) + + return; + } // separateNodes() + + + //================================================================================ + /*! + * \brief Find association of nodes existing on the sinuous sides of a ring + * + * TMAPar2NPoints filled here is used in setQuadSides() only if theSinuFace.IsRing() + * to find most distant nodes of the inner and the outer wires + */ + //================================================================================ + + void assocNodes( SMESH_MesherHelper& theHelper, + SinuousFace& theSinuFace, + const SMESH_MAT2d::MedialAxis& theMA, + TMAPar2NPoints & thePointsOnE ) + { + SMESH_Mesh* mesh = theHelper.GetMesh(); + SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); + + list< TopoDS_Edge > ee1( theSinuFace._sinuSide [0].begin(), theSinuFace._sinuSide [0].end() ); + list< TopoDS_Edge > ee2( theSinuFace._sinuSide [1].begin(), theSinuFace._sinuSide [1].end() ); + StdMeshers_FaceSide sideOut( theSinuFace.Face(), ee1, mesh, true, true, &theHelper ); + StdMeshers_FaceSide sideIn ( theSinuFace.Face(), ee2, mesh, true, true, &theHelper ); + const UVPtStructVec& uvsOut = sideOut.GetUVPtStruct(); + const UVPtStructVec& uvsIn = sideIn.GetUVPtStruct(); + // if ( uvs1.size() != uvs2.size() ) + // return; + + const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); + SMESH_MAT2d::BoundaryPoint bp[2]; + SMESH_MAT2d::BranchPoint brp; + + map< double, const SMDS_MeshNode* > nodeParams; // params of existing nodes + map< double, const SMDS_MeshNode* >::iterator u2n; + + // find a node of sideOut most distant from sideIn + + vector< BRepAdaptor_Curve > curvesIn( theSinuFace._sinuSide[1].size() ); + for ( size_t iE = 0; iE < theSinuFace._sinuSide[1].size(); ++iE ) + curvesIn[ iE ].Initialize( theSinuFace._sinuSide[1][iE] ); + + double maxDist = 0; + SMESH_MAT2d::BoundaryPoint bpIn; // closest IN point + const SMDS_MeshNode* nOut = 0; + const size_t nbEOut = theSinuFace._sinuSide[0].size(); + for ( size_t iE = 0; iE < nbEOut; ++iE ) + { + const TopoDS_Edge& E = theSinuFace._sinuSide[0][iE]; + + if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, E, /*skipMedium=*/true, nodeParams )) + return; + for ( u2n = nodeParams.begin(); u2n != nodeParams.end(); ++u2n ) + { + // point on EDGE (u2n) --> MA point (brp) + if ( !theMA.getBoundary().getBranchPoint( iE, u2n->first, brp ) || + !branch.getBoundaryPoints( brp, bp[0], bp[1] )) + return; + gp_Pnt pOut = SMESH_TNodeXYZ( u2n->second ); + gp_Pnt pIn = curvesIn[ bp[1]._edgeIndex - nbEOut ].Value( bp[1]._param ); + double dist = pOut.SquareDistance( pIn ); + if ( dist > maxDist ) + { + maxDist = dist; + nOut = u2n->second; + bpIn = bp[1]; + } + } + } + const SMDS_MeshNode* nIn = 0; + if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, + theSinuFace._sinuEdges[ bpIn._edgeIndex ], + /*skipMedium=*/true, + nodeParams )) + return; + u2n = nodeParams.lower_bound( bpIn._param ); + if ( u2n == nodeParams.end() ) + nIn = nodeParams.rbegin()->second; + else + nIn = u2n->second; + + // find position of distant nodes in uvsOut and uvsIn + size_t iDistOut, iDistIn; + for ( iDistOut = 0; iDistOut < uvsOut.size(); ++iDistOut ) + { + if ( uvsOut[iDistOut].node == nOut ) + break; + } + for ( iDistIn = 0; iDistIn < uvsIn.size(); ++iDistIn ) + { + if ( uvsIn[iDistIn].node == nIn ) + break; + } + if ( iDistOut == uvsOut.size() || iDistIn == uvsIn.size() ) + return; + + // store opposite nodes in thePointsOnE (param and EDGE have no sense) + pair< NodePoint, NodePoint > oppNodes( NodePoint( nOut, 0, 0 ), NodePoint( nIn, 0, 0)); + thePointsOnE.insert( make_pair( uvsOut[ iDistOut ].normParam, oppNodes )); + int iOut = iDistOut, iIn = iDistIn; + int i, nbNodes = std::min( uvsOut.size(), uvsIn.size() ); + if ( nbNodes > 5 ) nbNodes = 5; + for ( i = 0, ++iOut, --iIn; i < nbNodes; ++iOut, --iIn, ++i ) + { + iOut = theHelper.WrapIndex( iOut, uvsOut.size() ); + iIn = theHelper.WrapIndex( iIn, uvsIn.size() ); + oppNodes.first._node = uvsOut[ iOut ].node; + oppNodes.second._node = uvsIn[ iIn ].node; + thePointsOnE.insert( make_pair( uvsOut[ iOut ].normParam, oppNodes )); + } + + return; + } // assocNodes() + + //================================================================================ + /*! + * \brief Setup sides of SinuousFace::_quad + * \param [in] theHelper - helper + * \param [in] thePointsOnEdges - NodePoint's on sinuous sides + * \param [in,out] theSinuFace - the FACE + * \param [in] the1dAlgo - algorithm to use for radial discretization of a ring FACE + * \return bool - is OK + */ + //================================================================================ + + bool setQuadSides(SMESH_MesherHelper& theHelper, + const TMAPar2NPoints& thePointsOnEdges, + SinuousFace& theFace, + SMESH_Algo* the1dAlgo) + { + SMESH_Mesh* mesh = theHelper.GetMesh(); + const TopoDS_Face& face = theFace._quad->face; + SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( *mesh, face ); + if ( !proxyMesh ) + return false; + + list< TopoDS_Edge > side[4]; + side[0].insert( side[0].end(), theFace._shortSide[0].begin(), theFace._shortSide[0].end() ); + side[1].insert( side[1].end(), theFace._sinuSide [1].begin(), theFace._sinuSide [1].end() ); + side[2].insert( side[2].end(), theFace._shortSide[1].begin(), theFace._shortSide[1].end() ); + side[3].insert( side[3].end(), theFace._sinuSide [0].begin(), theFace._sinuSide [0].end() ); + + for ( int i = 0; i < 4; ++i ) + { + theFace._quad->side[i] = StdMeshers_FaceSide::New( face, side[i], mesh, i < QUAD_TOP_SIDE, + /*skipMediumNodes=*/true, + &theHelper, proxyMesh ); + } + + if ( theFace.IsRing() ) + { + // -------------------------------------- + // Discretize a ring in radial direction + // -------------------------------------- + + if ( thePointsOnEdges.size() < 4 ) + return false; + + int nbOut = theFace._quad->side[ 1 ].GetUVPtStruct().size(); + int nbIn = theFace._quad->side[ 3 ].GetUVPtStruct().size(); + if ( nbOut == 0 || nbIn == 0 ) + return false; + + // find most distant opposite nodes + double maxDist = 0, dist; + TMAPar2NPoints::const_iterator u2NPdist, u2NP = thePointsOnEdges.begin(); + for ( ; u2NP != thePointsOnEdges.end(); ++u2NP ) + { + SMESH_TNodeXYZ xyz( u2NP->second.first._node ); // node out + dist = xyz.SquareDistance( u2NP->second.second._node );// node in + if ( dist > maxDist ) + { + u2NPdist = u2NP; + maxDist = dist; + } + } + // compute distribution of radial nodes + list< double > params; // normalized params + static_cast< StdMeshers_QuadFromMedialAxis_1D2D::Algo1D* > + ( the1dAlgo )->ComputeDistribution( theHelper, + SMESH_TNodeXYZ( u2NPdist->second.first._node ), + SMESH_TNodeXYZ( u2NPdist->second.second._node ), + params ); + + // add a radial quad side + + theHelper.SetElementsOnShape( true ); + u2NP = thePointsOnEdges.begin(); + const SMDS_MeshNode* nOut = u2NP->second.first._node; + const SMDS_MeshNode* nIn = u2NP->second.second._node; + nOut = proxyMesh->GetProxyNode( nOut ); + nIn = proxyMesh->GetProxyNode( nIn ); + gp_XY uvOut = theHelper.GetNodeUV( face, nOut ); + gp_XY uvIn = theHelper.GetNodeUV( face, nIn ); + Handle(Geom_Surface) surface = BRep_Tool::Surface( face ); + UVPtStructVec uvsNew; UVPtStruct uvPt; + uvPt.node = nOut; + uvPt.u = uvOut.X(); + uvPt.v = uvOut.Y(); + uvsNew.push_back( uvPt ); + for (list::iterator itU = params.begin(); itU != params.end(); ++itU ) + { + gp_XY uv = ( 1 - *itU ) * uvOut + *itU * uvIn; // applied in direction Out -> In + gp_Pnt p = surface->Value( uv.X(), uv.Y() ); + uvPt.node = theHelper.AddNode( p.X(), p.Y(), p.Z(), /*id=*/0, uv.X(), uv.Y() ); + uvPt.u = uv.X(); + uvPt.v = uv.Y(); + uvsNew.push_back( uvPt ); + } + uvPt.node = nIn; + uvPt.u = uvIn.X(); + uvPt.v = uvIn.Y(); + uvsNew.push_back( uvPt ); + + theFace._quad->side[ 0 ] = StdMeshers_FaceSide::New( uvsNew ); + theFace._quad->side[ 2 ] = theFace._quad->side[ 0 ]; + if ( nbIn != nbOut ) + theFace._quad->side[ 2 ] = StdMeshers_FaceSide::New( uvsNew ); + + // assure that the outer sinuous side starts at nOut + { + const UVPtStructVec& uvsOut = theFace._quad->side[ 3 ].GetUVPtStruct(); // _sinuSide[0] + size_t i; // find UVPtStruct holding nOut + for ( i = 0; i < uvsOut.size(); ++i ) + if ( nOut == uvsOut[i].node ) + break; + if ( i == uvsOut.size() ) + return false; + + if ( i != 0 && i != uvsOut.size()-1 ) + { + // create a new OUT quad side + uvsNew.clear(); + uvsNew.reserve( uvsOut.size() ); + uvsNew.insert( uvsNew.end(), uvsOut.begin() + i, uvsOut.end() ); + uvsNew.insert( uvsNew.end(), uvsOut.begin() + 1, uvsOut.begin() + i + 1); + theFace._quad->side[ 3 ] = StdMeshers_FaceSide::New( uvsNew ); + } + } + + // rotate the IN side if opposite nodes of IN and OUT sides don't match + + const SMDS_MeshNode * nIn0 = theFace._quad->side[ 1 ].First().node; + if ( nIn0 != nIn ) + { + nIn = proxyMesh->GetProxyNode( nIn ); + const UVPtStructVec& uvsIn = theFace._quad->side[ 1 ].GetUVPtStruct(); // _sinuSide[1] + size_t i; // find UVPtStruct holding nIn + for ( i = 0; i < uvsIn.size(); ++i ) + if ( nIn == uvsIn[i].node ) + break; + if ( i == uvsIn.size() ) + return false; + + // create a new IN quad side + uvsNew.clear(); + uvsNew.reserve( uvsIn.size() ); + uvsNew.insert( uvsNew.end(), uvsIn.begin() + i, uvsIn.end() ); + uvsNew.insert( uvsNew.end(), uvsIn.begin() + 1, uvsIn.begin() + i + 1); + theFace._quad->side[ 1 ] = StdMeshers_FaceSide::New( uvsNew ); + } + + if ( theFace._quad->side[ 1 ].GetUVPtStruct().empty() || + theFace._quad->side[ 3 ].GetUVPtStruct().empty() ) + return false; + + } // if ( theFace.IsRing() ) + + return true; + + } // setQuadSides() + //================================================================================ /*! * \brief Divide the sinuous EDGEs by projecting the division point of Medial - * Axis to the EGDEs + * Axis to the EDGEs * \param [in] theHelper - the helper * \param [in] theMinSegLen - minimal segment length * \param [in] theMA - the Medial Axis * \param [in] theMAParams - parameters of division points of \a theMA * \param [in] theSinuEdges - the EDGEs to make nodes on * \param [in] theSinuSide0Size - the number of EDGEs in the 1st sinuous side + * \param [in] the1dAlgo - algorithm to use for radial discretization of a ring FACE * \return bool - is OK or not */ //================================================================================ @@ -983,9 +1659,10 @@ namespace double /*theMinSegLen*/, SMESH_MAT2d::MedialAxis& theMA, vector& theMAParams, - SinuousFace& theSinuFace) + SinuousFace& theSinuFace, + SMESH_Algo* the1dAlgo) { - if ( theMA.getBranches().size() != 1 ) + if ( theMA.nbBranches() != 1 ) return false; // normalize theMAParams @@ -997,11 +1674,13 @@ namespace SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); double f,l; + // get data of sinuous EDGEs and remove unnecessary nodes const vector< TopoDS_Edge >& theSinuEdges = theSinuFace._sinuEdges; - vector< Handle(Geom_Curve) > curves ( theSinuEdges.size() ); - vector< int > edgeIDs( theSinuEdges.size() ); - vector< bool > isComputed( theSinuEdges.size() ); - //bool hasComputed = false; + vector< Handle(Geom_Curve) >& curves = theSinuFace._sinuCurves; + vector< int > edgeIDs ( theSinuEdges.size() ); // IDs in the main shape + vector< bool > isComputed( theSinuEdges.size() ); + curves.resize( theSinuEdges.size(), 0 ); + bool allComputed = true; for ( size_t i = 0; i < theSinuEdges.size(); ++i ) { curves[i] = BRep_Tool::Curve( theSinuEdges[i], f,l ); @@ -1010,164 +1689,235 @@ namespace SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] ); edgeIDs [i] = sm->GetId(); isComputed[i] = ( !sm->IsEmpty() ); - if ( isComputed[i] ) - { - TopAbs_ShapeEnum shape = getHypShape( mesh, theSinuEdges[i] ); - if ( shape == TopAbs_SHAPE || shape <= TopAbs_FACE ) - { - // EDGE computed using global hypothesis -> clear it - bool hasComputedFace = false; - PShapeIteratorPtr faceIt = theHelper.GetAncestors( theSinuEdges[i], *mesh, TopAbs_FACE ); - while ( const TopoDS_Shape* face = faceIt->next() ) - if (( !face->IsSame( theSinuFace.Face())) && - ( hasComputedFace = !mesh->GetSubMesh( *face )->IsEmpty() )) - break; - if ( !hasComputedFace ) - sm->ComputeStateEngine( SMESH_subMesh::CLEAN ); - isComputed[i] = false; - } - } + if ( !isComputed[i] ) + allComputed = false; } - 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; - vector< SMESH_MAT2d::BranchPoint > divPoints; - branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints ); - for ( size_t i = 0; i < edgeIDs1.size(); ++i ) - if ( isComputed[ edgeIDs1[i]] && - isComputed[ edgeIDs2[i]]) - return false; - - // map param on MA to parameters of nodes on a pair of theSinuEdges - typedef map< double, pair< NodePoint, NodePoint > > TMAPar2NPoints; TMAPar2NPoints pointsOnE; - vector maParams; - - // compute params of nodes on EDGEs by projecting division points from MA - //const double tol = 1e-5 * theMAParams.back(); - size_t iEdgePair = 0; - while ( iEdgePair < edgeIDs1.size() ) + // check that computed EDGEs are opposite and equally meshed + if ( allComputed ) { - if ( isComputed[ edgeIDs1[ iEdgePair ]] || - isComputed[ edgeIDs2[ iEdgePair ]]) - { - // "projection" from one side to the other + // int nbNodes[2] = { 0, 0 }; + // for ( int iSide = 0; iSide < 2; ++iSide ) // loop on two sinuous sides + // nbNodes[ iSide ] += meshDS->MeshElements( theSinuFace._sinuSide[ iSide ])->NbNodes() - 1; - size_t iEdgeComputed = edgeIDs1[iEdgePair], iSideComputed = 0; - if ( !isComputed[ iEdgeComputed ]) - ++iSideComputed, iEdgeComputed = edgeIDs2[iEdgePair]; + // if ( nbNodes[0] != nbNodes[1] ) + // return false; - map< double, const SMDS_MeshNode* > nodeParams; // params of existing nodes - if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iEdgeComputed ], /*skipMedium=*/true, nodeParams )) - return false; - - SMESH_MAT2d::BoundaryPoint& bndPnt = bp[ 1-iSideComputed ]; - SMESH_MAT2d::BranchPoint brp; - NodePoint npN, npB; - NodePoint& np0 = iSideComputed ? npB : npN; - NodePoint& np1 = iSideComputed ? npN : npB; + if ( theSinuFace.IsRing() ) + assocNodes( theHelper, theSinuFace, theMA, pointsOnE ); + } + else + { + vector< std::size_t > edgeIDs1, edgeIDs2; // indices in theSinuEdges + vector< SMESH_MAT2d::BranchPoint > divPoints; + branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints ); - double maParam1st, maParamLast, maParam; - if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.begin()->first, brp )) + for ( size_t i = 0; i < edgeIDs1.size(); ++i ) + if ( isComputed[ edgeIDs1[i]] && + isComputed[ edgeIDs2[i]] ) + { + int nbNodes1 = meshDS->MeshElements(edgeIDs[ edgeIDs1[i]] )->NbNodes(); + int nbNodes2 = meshDS->MeshElements(edgeIDs[ edgeIDs2[i]] )->NbNodes(); + if ( nbNodes1 != nbNodes2 ) + return false; + if (( int(i)-1 >= 0 ) && + ( edgeIDs1[i-1] == edgeIDs1[i] || + edgeIDs2[i-1] == edgeIDs2[i] )) return false; - branch.getParameter( brp, maParam1st ); - if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.rbegin()->first, brp )) + if (( i+1 < edgeIDs1.size() ) && + ( edgeIDs1[i+1] == edgeIDs1[i] || + edgeIDs2[i+1] == edgeIDs2[i] )) return false; - branch.getParameter( brp, maParamLast ); + } + + // map (param on MA) to (parameters of nodes on a pair of theSinuEdges) + vector maParams; + set projectedEdges; // treated EDGEs which 'isComputed' + + // compute params of nodes on EDGEs by projecting division points from MA - map< double, const SMDS_MeshNode* >::iterator u2n = nodeParams.begin(), u2nEnd = --nodeParams.end(); - TMAPar2NPoints::iterator end = pointsOnE.end(), pos = end; - TMAPar2NPoints::iterator & hint = (maParamLast > maParam1st) ? end : pos; - for ( ++u2n; u2n != u2nEnd; ++u2n ) + for ( size_t iEdgePair = 0; iEdgePair < edgeIDs1.size(); ++iEdgePair ) + // loop on pairs of opposite EDGEs + { + if ( projectedEdges.count( edgeIDs1[ iEdgePair ]) || + projectedEdges.count( edgeIDs2[ iEdgePair ]) ) + continue; + + // -------------------------------------------------------------------------------- + if ( isComputed[ edgeIDs1[ iEdgePair ]] != // one EDGE is meshed + isComputed[ edgeIDs2[ iEdgePair ]]) { - if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, u2n->first, brp )) + // "projection" from one side to the other + + size_t iEdgeComputed = edgeIDs1[iEdgePair], iSideComputed = 0; + if ( !isComputed[ iEdgeComputed ]) + ++iSideComputed, iEdgeComputed = edgeIDs2[iEdgePair]; + + map< double, const SMDS_MeshNode* > nodeParams; // params of existing nodes + if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iEdgeComputed ], /*skipMedium=*/true, nodeParams )) return false; - if ( !branch.getBoundaryPoints( brp, bp[0], bp[1] )) + + projectedEdges.insert( iEdgeComputed ); + + SMESH_MAT2d::BoundaryPoint& bndPnt = bp[ 1-iSideComputed ]; + SMESH_MAT2d::BranchPoint brp; + NodePoint npN, npB; // NodePoint's initialized by node and BoundaryPoint + NodePoint& np0 = iSideComputed ? npB : npN; + NodePoint& np1 = iSideComputed ? npN : npB; + + double maParam1st, maParamLast, maParam; + if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.begin()->first, brp )) return false; - if ( !branch.getParameter( brp, maParam )) + branch.getParameter( brp, maParam1st ); + if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.rbegin()->first, brp )) return false; + branch.getParameter( brp, maParamLast ); - npN = NodePoint( u2n->second, u2n->first, iEdgeComputed ); - npB = NodePoint( bndPnt ); - pos = pointsOnE.insert( hint, make_pair( maParam, make_pair( np0, np1 ))); + map< double, const SMDS_MeshNode* >::iterator u2n = nodeParams.begin(), u2nEnd = nodeParams.end(); + TMAPar2NPoints::iterator end = pointsOnE.end(), pos = end; + TMAPar2NPoints::iterator & hint = (maParamLast > maParam1st) ? end : pos; + for ( ++u2n, --u2nEnd; u2n != u2nEnd; ++u2n ) + { + // point on EDGE (u2n) --> MA point (brp) + if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, u2n->first, brp )) + return false; + // MA point --> points on 2 EDGEs (bp) + if ( !branch.getBoundaryPoints( brp, bp[0], bp[1] ) || + !branch.getParameter( brp, maParam )) + return false; + + npN = NodePoint( u2n->second, u2n->first, iEdgeComputed ); + npB = NodePoint( bndPnt ); + pos = pointsOnE.insert( hint, make_pair( maParam, make_pair( np0, np1 ))); + } } + // -------------------------------------------------------------------------------- + else if ( !isComputed[ edgeIDs1[ iEdgePair ]] && // none of EDGEs is meshed + !isComputed[ edgeIDs2[ iEdgePair ]]) + { + // "projection" from MA + maParams.clear(); + if ( !getParamsForEdgePair( iEdgePair, divPoints, theMAParams, maParams )) + return false; - // move iEdgePair forward - while ( iEdgePair < edgeIDs1.size() ) - if ( edgeIDs1[ iEdgePair ] == bp[0]._edgeIndex && - edgeIDs2[ iEdgePair ] == bp[1]._edgeIndex ) - break; - else - ++iEdgePair; - } - else - { - // projection from MA - maParams.clear(); - if ( !getParamsForEdgePair( iEdgePair, divPoints, theMAParams, maParams )) - return false; + for ( size_t i = 1; i < maParams.size()-1; ++i ) + { + if ( !branch.getBoundaryPoints( maParams[i], bp[0], bp[1] )) + return false; - for ( size_t i = 1; i < maParams.size()-1; ++i ) + pointsOnE.insert( pointsOnE.end(), make_pair( maParams[i], make_pair( NodePoint(bp[0]), + NodePoint(bp[1])))); + } + } + // -------------------------------------------------------------------------------- + else if ( isComputed[ edgeIDs1[ iEdgePair ]] && // equally meshed EDGES + isComputed[ edgeIDs2[ iEdgePair ]]) { - if ( !branch.getBoundaryPoints( maParams[i], bp[0], bp[1] )) + // add existing nodes + + size_t iE0 = edgeIDs1[ iEdgePair ]; + size_t iE1 = edgeIDs2[ iEdgePair ]; + map< double, const SMDS_MeshNode* > nodeParams[2]; // params of existing nodes + if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iE0 ], + /*skipMedium=*/false, nodeParams[0] ) || + !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iE1 ], + /*skipMedium=*/false, nodeParams[1] ) || + nodeParams[0].size() != nodeParams[1].size() ) return false; - pointsOnE.insert( pointsOnE.end(), make_pair( maParams[i], make_pair( NodePoint(bp[0]), - NodePoint(bp[1])))); + if ( nodeParams[0].size() <= 2 ) + continue; // nodes on VERTEXes only + + bool reverse = ( theSinuEdges[0].Orientation() == theSinuEdges[1].Orientation() ); + double maParam; + SMESH_MAT2d::BranchPoint brp; + std::pair< NodePoint, NodePoint > npPair; + + map< double, const SMDS_MeshNode* >::iterator + u2n0F = ++nodeParams[0].begin(), + u2n1F = ++nodeParams[1].begin(); + map< double, const SMDS_MeshNode* >::reverse_iterator + u2n1R = ++nodeParams[1].rbegin(); + for ( ; u2n0F != nodeParams[0].end(); ++u2n0F ) + { + if ( !theMA.getBoundary().getBranchPoint( iE0, u2n0F->first, brp ) || + !branch.getParameter( brp, maParam )) + return false; + + npPair.first = NodePoint( u2n0F->second, u2n0F->first, iE0 ); + if ( reverse ) + { + npPair.second = NodePoint( u2n1R->second, u2n1R->first, iE1 ); + ++u2n1R; + } + else + { + npPair.second = NodePoint( u2n1F->second, u2n1F->first, iE1 ); + ++u2n1F; + } + pointsOnE.insert( make_pair( maParam, npPair )); + } } - } - ++iEdgePair; - } + } // loop on pairs of opposite EDGEs - if ( !projectVertices( theHelper, theMA, divPoints, theSinuEdges, curves, - isComputed, pointsOnE, theSinuFace._nodesToMerge )) - return false; + if ( !projectVertices( theHelper, theMA, divPoints, edgeIDs1, edgeIDs2, + isComputed, pointsOnE, theSinuFace )) + return false; - // create nodes - TMAPar2NPoints::iterator u2np = pointsOnE.begin(); - for ( ; u2np != pointsOnE.end(); ++u2np ) - { - NodePoint* np[2] = { & u2np->second.first, & u2np->second.second }; - for ( int iSide = 0; iSide < 2; ++iSide ) + separateNodes( theHelper, theMA, pointsOnE, theSinuFace, isComputed ); + + // create nodes + TMAPar2NPoints::iterator u2np = pointsOnE.begin(); + for ( ; u2np != pointsOnE.end(); ++u2np ) { - if ( np[ iSide ]->_node ) continue; - size_t iEdge = np[ iSide ]->_edgeInd; - double u = np[ iSide ]->_u; - gp_Pnt p = curves[ iEdge ]->Value( u ); - np[ iSide ]->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); - meshDS->SetNodeOnEdge( np[ iSide ]->_node, edgeIDs[ iEdge ], u ); + NodePoint* np[2] = { & u2np->second.first, & u2np->second.second }; + for ( int iSide = 0; iSide < 2; ++iSide ) + { + if ( np[ iSide ]->_node ) continue; + size_t iEdge = np[ iSide ]->_edgeInd; + double u = np[ iSide ]->_u; + gp_Pnt p = curves[ iEdge ]->Value( u ); + np[ iSide ]->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); + meshDS->SetNodeOnEdge( np[ iSide ]->_node, edgeIDs[ iEdge ], u ); + } } - } - // create mesh segments on EDGEs - theHelper.SetElementsOnShape( false ); - TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() ); - for ( size_t i = 0; i < theSinuEdges.size(); ++i ) - { - SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] ); - if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 ) - continue; + // create mesh segments on EDGEs + theHelper.SetElementsOnShape( false ); + TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() ); + for ( size_t i = 0; i < theSinuEdges.size(); ++i ) + { + SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] ); + if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 ) + continue; + + StdMeshers_FaceSide side( face, theSinuEdges[i], mesh, + /*isFwd=*/true, /*skipMediumNodes=*/true, &theHelper ); + vector nodes = side.GetOrderedNodes(); + for ( size_t in = 1; in < nodes.size(); ++in ) + { + const SMDS_MeshElement* seg = theHelper.AddEdge( nodes[in-1], nodes[in], 0, false ); + meshDS->SetMeshElementOnShape( seg, edgeIDs[ i ] ); + } + } - StdMeshers_FaceSide side( face, theSinuEdges[i], mesh, - /*isFwd=*/true, /*skipMediumNodes=*/true ); - vector nodes = side.GetOrderedNodes(); - for ( size_t in = 1; in < nodes.size(); ++in ) + // update sub-meshes on VERTEXes + for ( size_t i = 0; i < theSinuEdges.size(); ++i ) { - const SMDS_MeshElement* seg = theHelper.AddEdge( nodes[in-1], nodes[in], 0, false ); - meshDS->SetMeshElementOnShape( seg, edgeIDs[ i ] ); + mesh->GetSubMesh( theHelper.IthVertex( 0, theSinuEdges[i] )) + ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + mesh->GetSubMesh( theHelper.IthVertex( 1, theSinuEdges[i] )) + ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); } } - // update sub-meshes on VERTEXes - for ( size_t i = 0; i < theSinuEdges.size(); ++i ) - { - mesh->GetSubMesh( theHelper.IthVertex( 0, theSinuEdges[i] )) - ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); - mesh->GetSubMesh( theHelper.IthVertex( 1, theSinuEdges[i] )) - ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); - } + // Setup sides of a quadrangle + if ( !setQuadSides( theHelper, pointsOnE, theSinuFace, the1dAlgo )) + return false; return true; } @@ -1180,16 +1930,30 @@ namespace bool computeShortEdges( SMESH_MesherHelper& theHelper, const vector& theShortEdges, - SMESH_Algo* the1dAlgo ) + SMESH_Algo* the1dAlgo, + const bool theHasRadialHyp, + const bool theIs2nd) { + SMESH_Hypothesis::Hypothesis_Status aStatus; for ( size_t i = 0; i < theShortEdges.size(); ++i ) { - theHelper.GetGen()->Compute( *theHelper.GetMesh(), theShortEdges[i], true, true ); + if ( !theHasRadialHyp ) + // use global hyps + theHelper.GetGen()->Compute( *theHelper.GetMesh(), theShortEdges[i], + SMESH_Gen::SHAPE_ONLY_UPWARD ); SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh(theShortEdges[i] ); if ( sm->IsEmpty() ) { + // use 2D hyp or minSegLen try { + // compute VERTEXes + SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false); + while ( smIt->more() ) + smIt->next()->ComputeStateEngine( SMESH_subMesh::COMPUTE ); + + // compute EDGE + the1dAlgo->CheckHypothesis( *theHelper.GetMesh(), theShortEdges[i], aStatus ); if ( !the1dAlgo->Compute( *theHelper.GetMesh(), theShortEdges[i] )) return false; } @@ -1223,12 +1987,12 @@ namespace const double dksi = 0.5, deta = 0.5; const double dksi2 = dksi*dksi, deta2 = deta*deta; double err = 0., g11, g22, g12; - int nbErr = 0; + //int nbErr = 0; FaceQuadStruct& q = *quad; UVPtStruct pNew; - double refArea = area( q.UVPt(0,0), q.UVPt(1,0), q.UVPt(1,1) ); + //double refArea = area( q.UVPt(0,0), q.UVPt(1,0), q.UVPt(1,1) ); for ( int iLoop = 0; iLoop < nbLoops; ++iLoop ) { @@ -1302,6 +2066,7 @@ namespace TMergeMap::iterator n2nn = theSinuFace._nodesToMerge.begin(); for ( ; n2nn != theSinuFace._nodesToMerge.end(); ++n2nn ) { + if ( !n2nn->first ) continue; nodesGroups.push_back( list< const SMDS_MeshNode* >() ); list< const SMDS_MeshNode* > & group = nodesGroups.back(); @@ -1313,6 +2078,17 @@ namespace } // namespace +//================================================================================ +/*! + * \brief Sets event listener which removes mesh from EDGEs when 2D hyps change + */ +//================================================================================ + +void StdMeshers_QuadFromMedialAxis_1D2D::SetEventListener(SMESH_subMesh* faceSubMesh) +{ + faceSubMesh->SetEventListener( new EdgeCleaner, 0, faceSubMesh ); +} + //================================================================================ /*! * \brief Create quadrangle elements @@ -1324,57 +2100,41 @@ namespace */ //================================================================================ -bool StdMeshers_QuadFromMedialAxis_1D2D::computeQuads( SMESH_MesherHelper& theHelper, - const TopoDS_Face& theFace, - const vector theSinuEdges[2], - const vector theShortEdges[2]) +bool StdMeshers_QuadFromMedialAxis_1D2D::computeQuads( SMESH_MesherHelper& theHelper, + FaceQuadStruct::Ptr theQuad) { - SMESH_Mesh* mesh = theHelper.GetMesh(); - SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( *mesh, theFace ); - if ( !proxyMesh ) - return false; - - StdMeshers_Quadrangle_2D::myProxyMesh = proxyMesh; StdMeshers_Quadrangle_2D::myHelper = &theHelper; StdMeshers_Quadrangle_2D::myNeedSmooth = false; StdMeshers_Quadrangle_2D::myCheckOri = false; StdMeshers_Quadrangle_2D::myQuadList.clear(); - // fill FaceQuadStruct - - list< TopoDS_Edge > side[4]; - side[0].insert( side[0].end(), theShortEdges[0].begin(), theShortEdges[0].end() ); - side[1].insert( side[1].end(), theSinuEdges[1].begin(), theSinuEdges[1].end() ); - side[2].insert( side[2].end(), theShortEdges[1].begin(), theShortEdges[1].end() ); - side[3].insert( side[3].end(), theSinuEdges[0].begin(), theSinuEdges[0].end() ); - - FaceQuadStruct::Ptr quad( new FaceQuadStruct ); - quad->side.resize( 4 ); - quad->face = theFace; - for ( int i = 0; i < 4; ++i ) - { - quad->side[i] = StdMeshers_FaceSide::New( theFace, side[i], mesh, i < QUAD_TOP_SIDE, - /*skipMediumNodes=*/true, proxyMesh ); - } - int nbNodesShort0 = quad->side[0].NbPoints(); - int nbNodesShort1 = quad->side[2].NbPoints(); + int nbNodesShort0 = theQuad->side[0].NbPoints(); + int nbNodesShort1 = theQuad->side[2].NbPoints(); + int nbNodesSinu0 = theQuad->side[1].NbPoints(); + int nbNodesSinu1 = theQuad->side[3].NbPoints(); // compute UV of internal points - myQuadList.push_back( quad ); - if ( !StdMeshers_Quadrangle_2D::setNormalizedGrid( quad )) - return false; + myQuadList.push_back( theQuad ); + // if ( !StdMeshers_Quadrangle_2D::setNormalizedGrid( theQuad )) + // return false; // elliptic smooth of internal points to get boundary cell normal to the boundary - ellipticSmooth( quad, 1 ); - + bool isRing = theQuad->side[0].grid->Edge(0).IsNull(); + if ( !isRing ) { + if ( !StdMeshers_Quadrangle_2D::setNormalizedGrid( theQuad )) + return false; + ellipticSmooth( theQuad, 1 ); + } // create quadrangles bool ok; - if ( nbNodesShort0 == nbNodesShort1 ) - ok = StdMeshers_Quadrangle_2D::computeQuadDominant( *mesh, theFace, quad ); + theHelper.SetElementsOnShape( true ); + if ( nbNodesShort0 == nbNodesShort1 && nbNodesSinu0 == nbNodesSinu1 ) + ok = StdMeshers_Quadrangle_2D::computeQuadDominant( *theHelper.GetMesh(), + theQuad->face, theQuad ); else - ok = StdMeshers_Quadrangle_2D::computeTriangles( *mesh, theFace, quad ); + ok = StdMeshers_Quadrangle_2D::computeTriangles( *theHelper.GetMesh(), + theQuad->face, theQuad ); - StdMeshers_Quadrangle_2D::myProxyMesh.reset(); StdMeshers_Quadrangle_2D::myHelper = 0; return ok; @@ -1401,36 +2161,35 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::Compute(SMESH_Mesh& theMesh, if ( getSinuousEdges( helper, sinuFace )) { - _progress = 0.2; - - // if ( sinuFace._sinuEdges.size() > 2 ) - // return error(COMPERR_BAD_SHAPE, "Not yet supported case" ); + _progress = 0.4; double minSegLen = getMinSegLen( helper, sinuFace._sinuEdges ); SMESH_MAT2d::MedialAxis ma( F, sinuFace._sinuEdges, minSegLen, /*ignoreCorners=*/true ); if ( !_regular1D ) - _regular1D = new Algo1D( _studyId, _gen ); + _regular1D = new Algo1D( _gen ); _regular1D->SetSegmentLength( minSegLen ); vector maParams; - if ( ! divideMA( helper, ma, sinuFace, _regular1D, maParams )) + if ( ! divideMA( helper, ma, sinuFace, _regular1D, minSegLen, maParams )) return error(COMPERR_BAD_SHAPE); - _progress = 0.4; + _progress = 0.8; + if ( _hyp2D ) + _regular1D->SetRadialDistribution( _hyp2D ); - if ( !computeShortEdges( helper, sinuFace._shortSide[0], _regular1D ) || - !computeShortEdges( helper, sinuFace._shortSide[1], _regular1D )) + if ( !computeShortEdges( helper, sinuFace._shortSide[0], _regular1D, _hyp2D, 0 ) || + !computeShortEdges( helper, sinuFace._shortSide[1], _regular1D, _hyp2D, 1 )) return error("Failed to mesh short edges"); - _progress = 0.6; + _progress = 0.85; - if ( !computeSinuEdges( helper, minSegLen, ma, maParams, sinuFace )) + if ( !computeSinuEdges( helper, minSegLen, ma, maParams, sinuFace, _regular1D )) return error("Failed to mesh sinuous edges"); - _progress = 0.8; + _progress = 0.9; - bool ok = computeQuads( helper, F, sinuFace._sinuSide, sinuFace._shortSide ); + bool ok = computeQuads( helper, sinuFace._quad ); if ( ok ) mergeNodes( helper, sinuFace ); @@ -1456,3 +2215,31 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::Evaluate(SMESH_Mesh & theMesh, return StdMeshers_Quadrangle_2D::Evaluate(theMesh,theShape,theResMap); } +//================================================================================ +/*! + * \brief Return true if the algorithm can mesh this shape + * \param [in] aShape - shape to check + * \param [in] toCheckAll - if true, this check returns OK if all shapes are OK, + * else, returns OK if at least one shape is OK + */ +//================================================================================ + +bool StdMeshers_QuadFromMedialAxis_1D2D::IsApplicable( const TopoDS_Shape & aShape, + bool toCheckAll ) +{ + TmpMesh tmpMesh; + SMESH_MesherHelper helper( tmpMesh ); + + int nbFoundFaces = 0; + for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces ) + { + const TopoDS_Face& face = TopoDS::Face( exp.Current() ); + SinuousFace sinuFace( face ); + bool isApplicable = getSinuousEdges( helper, sinuFace ); + + if ( toCheckAll && !isApplicable ) return false; + if ( !toCheckAll && isApplicable ) return true; + } + return ( toCheckAll && nbFoundFaces != 0 ); +} +