X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_RadialQuadrangle_1D2D.cxx;h=d03470fb341870b223f48e6922217376695d5753;hp=95423e861795c5e0a07542c5d66bf2aa99a26033;hb=aa67cc96d730566d202d1014d97b7f0b3a4d71f4;hpb=3add66f2196c5b8e04e69d87c9f6a28e9a9a2eb2 diff --git a/src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx b/src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx index 95423e861..d03470fb3 100644 --- a/src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx +++ b/src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx @@ -1,7 +1,4 @@ -// Copyright (C) 2007-2008 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 +// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -19,6 +16,7 @@ // // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // + // SMESH SMESH : implementaion of SMESH idl descriptions // File : StdMeshers_RadialQuadrangle_1D2D.cxx // Module : SMESH @@ -27,15 +25,15 @@ // #include "StdMeshers_RadialQuadrangle_1D2D.hxx" -//#include "StdMeshers_ProjectionUtils.hxx" #include "StdMeshers_NumberOfLayers.hxx" #include "StdMeshers_LayerDistribution.hxx" -//#include "StdMeshers_Prism_3D.hxx" #include "StdMeshers_Regular_1D.hxx" +#include "StdMeshers_NumberOfSegments.hxx" #include "SMDS_MeshNode.hxx" #include "SMESHDS_SubMesh.hxx" #include "SMESH_Gen.hxx" +#include "SMESH_HypoFilter.hxx" #include "SMESH_Mesh.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_subMesh.hxx" @@ -45,23 +43,16 @@ #include #include -//#include #include -#include -#include -//#include -//#include -//#include -//#include -//#include - - -#include +#include #include #include +#include #include #include -#include +#include +#include +#include using namespace std; @@ -69,8 +60,6 @@ using namespace std; #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; } #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z()) -//typedef StdMeshers_ProjectionUtils TAssocTool; - //======================================================================= //function : StdMeshers_RadialQuadrangle_1D2D @@ -90,6 +79,7 @@ StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int hypId, myNbLayerHypo = 0; myDistributionHypo = 0; _requireDescretBoundary = false; + _supportSubmeshes = true; } @@ -121,8 +111,8 @@ bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis const list &hyps = GetUsedHypothesis(aMesh, aShape); if ( hyps.size() == 0 ) { - aStatus = SMESH_Hypothesis::HYP_MISSING; - return false; // can't work with no hypothesis + aStatus = SMESH_Hypothesis::HYP_OK; + return true; // can work with no hypothesis } if ( hyps.size() > 1 ) { @@ -148,13 +138,292 @@ bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis return true; } +namespace +{ + // ------------------------------------------------------------------------------ + /*! + * \brief Listener used to mark edges meshed by StdMeshers_RadialQuadrangle_1D2D + */ + class TEdgeMarker : public SMESH_subMeshEventListener + { + TEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false) {} + public: + //!< Return static listener + static SMESH_subMeshEventListener* getListener() + { + static TEdgeMarker theEdgeMarker; + return &theEdgeMarker; + } + //! Clear face sumbesh if something happens on edges + void ProcessEvent(const int event, + const int eventType, + SMESH_subMesh* edgeSubMesh, + EventListenerData* data, + const SMESH_Hypothesis* /*hyp*/) + { + if ( data && !data->mySubMeshes.empty() && eventType == SMESH_subMesh::ALGO_EVENT) + { + ASSERT( data->mySubMeshes.front() != edgeSubMesh ); + SMESH_subMesh* faceSubMesh = data->mySubMeshes.front(); + faceSubMesh->ComputeStateEngine( SMESH_subMesh::CLEAN ); + } + } + }; + + // ------------------------------------------------------------------------------ + /*! + * \brief Mark an edge as computed by StdMeshers_RadialQuadrangle_1D2D + */ + void markEdgeAsComputedByMe(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh) + { + if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge )) + { + if ( !edgeSM->GetEventListenerData( TEdgeMarker::getListener() )) + faceSubMesh->SetEventListener( TEdgeMarker::getListener(), + SMESH_subMeshEventListenerData::MakeData(faceSubMesh), + edgeSM); + } + } + // ------------------------------------------------------------------------------ + /*! + * \brief Return true if a radial edge was meshed with StdMeshers_RadialQuadrangle_1D2D with + * the same radial distribution + */ +// bool isEdgeCompatiballyMeshed(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh) +// { +// if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge )) +// { +// if ( SMESH_subMeshEventListenerData* otherFaceData = +// edgeSM->GetEventListenerData( TEdgeMarker::getListener() )) +// { +// // compare hypothesis aplied to two disk faces sharing radial edges +// SMESH_Mesh& mesh = *faceSubMesh->GetFather(); +// SMESH_Algo* radialQuadAlgo = mesh.GetGen()->GetAlgo(mesh, faceSubMesh->GetSubShape() ); +// SMESH_subMesh* otherFaceSubMesh = otherFaceData->mySubMeshes.front(); +// list hyps1 = +// radialQuadAlgo->GetUsedHypothesis( mesh, faceSubMesh->GetSubShape()); +// list hyps2 = +// radialQuadAlgo->GetUsedHypothesis( mesh, otherFaceSubMesh->GetSubShape()); +// if( hyps1.empty() && hyps2.empty() ) +// return true; // defaul hyps +// if ( hyps1.size() != hyps2.size() ) +// return false; +// return *hyps1.front() == *hyps2.front(); +// } +// } +// return false; +// } + + //================================================================================ + /*! + * \brief Return base curve of the edge and extremum parameters + */ + //================================================================================ + + Handle(Geom_Curve) getCurve(const TopoDS_Edge& edge, double* f=0, double* l=0) + { + Handle(Geom_Curve) C; + if ( !edge.IsNull() ) + { + double first = 0., last = 0.; + C = BRep_Tool::Curve(edge, first, last); + if ( !C.IsNull() ) + { + Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C); + while( !tc.IsNull() ) { + C = tc->BasisCurve(); + tc = Handle(Geom_TrimmedCurve)::DownCast(C); + } + if ( f ) *f = first; + if ( l ) *l = last; + } + } + return C; + } + + //================================================================================ + /*! + * \brief Return edges of the face + * \retval int - nb of edges + */ + //================================================================================ + + int analyseFace(const TopoDS_Shape& face, + TopoDS_Edge& CircEdge, + TopoDS_Edge& LinEdge1, + TopoDS_Edge& LinEdge2) + { + CircEdge.Nullify(); LinEdge1.Nullify(); LinEdge2.Nullify(); + int nbe = 0; + + for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next(), ++nbe ) + { + const TopoDS_Edge& E = TopoDS::Edge( exp.Current() ); + double f,l; + Handle(Geom_Curve) C = getCurve(E,&f,&l); + if ( !C.IsNull() ) + { + if ( C->IsKind( STANDARD_TYPE(Geom_Circle))) + { + if ( CircEdge.IsNull() ) + CircEdge = E; + else + return 0; + } + else if ( LinEdge1.IsNull() ) + LinEdge1 = E; + else + LinEdge2 = E; + } + } + return nbe; + } + +//================================================================================ +//================================================================================ +/*! + * \brief Class computing layers distribution using data of + * StdMeshers_LayerDistribution hypothesis + */ +//================================================================================ +//================================================================================ + +class TNodeDistributor: public StdMeshers_Regular_1D +{ + list myUsedHyps; +public: + // ----------------------------------------------------------------------------- + static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh) + { + const int myID = -1000; + map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo; + map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID ); + if ( id_algo == algoMap.end() ) + return new TNodeDistributor( myID, 0, aMesh.GetGen() ); + return static_cast< TNodeDistributor* >( id_algo->second ); + } + // ----------------------------------------------------------------------------- + //! Computes distribution of nodes on a straight line ending at pIn and pOut + bool Compute( vector< double > & positions, + gp_Pnt pIn, + gp_Pnt pOut, + SMESH_Mesh& aMesh, + const SMESH_Hypothesis* hyp1d) + { + if ( !hyp1d ) return error( "Invalid LayerDistribution hypothesis"); + + double len = pIn.Distance( pOut ); + if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells"); + + myUsedHyps.clear(); + myUsedHyps.push_back( hyp1d ); + + TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut ); + SMESH_Hypothesis::Hypothesis_Status aStatus; + if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus )) + return error( "StdMeshers_Regular_1D::CheckHypothesis() failed " + "with LayerDistribution hypothesis"); + + BRepAdaptor_Curve C3D(edge); + double f = C3D.FirstParameter(), l = C3D.LastParameter(); + list< double > params; + if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false )) + return error("StdMeshers_Regular_1D failed to compute layers distribution"); + + positions.clear(); + positions.reserve( params.size() ); + for (list::iterator itU = params.begin(); itU != params.end(); itU++) + positions.push_back( *itU / len ); + return true; + } + // ----------------------------------------------------------------------------- + //! Make mesh on an adge using assigned 1d hyp or defaut nb of segments + bool ComputeCircularEdge(SMESH_Mesh& aMesh, + const TopoDS_Edge& anEdge) + { + _gen->Compute( aMesh, anEdge); + SMESH_subMesh *sm = aMesh.GetSubMesh(anEdge); + if ( sm->GetComputeState() != SMESH_subMesh::COMPUTE_OK) + { + // find any 1d hyp assigned (there can be a hyp w/o algo) + myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true); + Hypothesis_Status aStatus; + if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus )) + { + // no valid 1d hyp assigned, use default nb of segments + _hypType = NB_SEGMENTS; + _ivalue[ DISTR_TYPE_IND ] = StdMeshers_NumberOfSegments::DT_Regular; + _ivalue[ NB_SEGMENTS_IND ] = _gen->GetDefaultNbSegments(); + } + return StdMeshers_Regular_1D::Compute( aMesh, anEdge ); + } + return true; + } + // ----------------------------------------------------------------------------- + //! Make mesh on an adge using assigned 1d hyp or defaut nb of segments + bool EvaluateCircularEdge(SMESH_Mesh& aMesh, + const TopoDS_Edge& anEdge, + MapShapeNbElems& aResMap) + { + _gen->Evaluate( aMesh, anEdge, aResMap ); + if ( aResMap.count( aMesh.GetSubMesh( anEdge ))) + return true; + + // find any 1d hyp assigned + myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true); + Hypothesis_Status aStatus; + if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus )) + { + // no valid 1d hyp assigned, use default nb of segments + _hypType = NB_SEGMENTS; + _ivalue[ DISTR_TYPE_IND ] = StdMeshers_NumberOfSegments::DT_Regular; + _ivalue[ NB_SEGMENTS_IND ] = _gen->GetDefaultNbSegments(); + } + return StdMeshers_Regular_1D::Evaluate( aMesh, anEdge, aResMap ); + } +protected: + // ----------------------------------------------------------------------------- + TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen) + : StdMeshers_Regular_1D( hypId, studyId, gen) + { + } + // ----------------------------------------------------------------------------- + virtual const list & + GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool) + { + return myUsedHyps; + } + // ----------------------------------------------------------------------------- +}; +} + +//======================================================================= +/*! + * \brief Allow algo to do something after persistent restoration + * \param subMesh - restored submesh + * + * call markEdgeAsComputedByMe() + */ +//======================================================================= + +void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh) +{ + if ( !faceSubMesh->IsEmpty() ) + { + TopoDS_Edge CircEdge, LinEdge1, LinEdge2; + analyseFace( faceSubMesh->GetSubShape(), CircEdge, LinEdge1, LinEdge2 ); + if ( !CircEdge.IsNull() ) markEdgeAsComputedByMe( CircEdge, faceSubMesh ); + if ( !LinEdge1.IsNull() ) markEdgeAsComputedByMe( LinEdge1, faceSubMesh ); + if ( !LinEdge2.IsNull() ) markEdgeAsComputedByMe( LinEdge2, faceSubMesh ); + } +} //======================================================================= //function : Compute //purpose : //======================================================================= -bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh, +bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) { TopExp_Explorer exp; @@ -162,82 +431,63 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh, myHelper = new SMESH_MesherHelper( aMesh ); myHelper->IsQuadraticSubMesh( aShape ); + // to delete helper at exit from Compute() + auto_ptr helperDeleter( myHelper ); - myLayerPositions.clear(); - - TopoDS_Edge E1,E2,E3; - Handle(Geom_Curve) C1,C2,C3; - double f1,l1,f2,l2,f3,l3; - int nbe = 0; - for ( exp.Init( aShape, TopAbs_EDGE ); exp.More(); exp.Next() ) { - nbe++; - TopoDS_Edge E = TopoDS::Edge( exp.Current() ); - if(nbe==1) { - E1 = E; - C1 = BRep_Tool::Curve(E,f1,l1); - } - else if(nbe==2) { - E2 = E; - C2 = BRep_Tool::Curve(E,f2,l2); - } - else if(nbe==3) { - E3 = E; - C3 = BRep_Tool::Curve(E,f3,l3); - } - } - - if(nbe>3) - return error(COMPERR_BAD_SHAPE); + TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh); - gp_Pnt P0,P1; + TopoDS_Edge CircEdge, LinEdge1, LinEdge2; + int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 ); + Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge )); + if( nbe>3 || nbe < 1 || aCirc.IsNull() ) + return error("The face must be a full circle or a part of circle (i.e. the number of edges is less or equal to 3 and one of them is a circle curve)"); + + gp_Pnt P0, P1; // points for rotation TColgp_SequenceOfPnt Points; // angles for rotation TColStd_SequenceOfReal Angles; // Nodes1 and Nodes2 - nodes along radiuses // CNodes - nodes on circle edge - std::vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes; + vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes; SMDS_MeshNode * NC; // parameters edge nodes on face - TColgp_SequenceOfPnt2d Pnts2d1, Pnts2d2; + TColgp_SequenceOfPnt2d Pnts2d1; gp_Pnt2d PC; int faceID = meshDS->ShapeToIndex(aShape); TopoDS_Face F = TopoDS::Face(aShape); Handle(Geom_Surface) S = BRep_Tool::Surface(F); - //cout<<"RadialQuadrangle_1D2D::Compute nbe = "<Compute( aMesh, CircEdge, false, MeshDim_1D ); - if( !ok ) return false; - std::map< double, const SMDS_MeshNode* > theNodes; - GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes); + if(nbe==1) + { + if (!algo1d->ComputeCircularEdge( aMesh, CircEdge )) + return error( algo1d->GetComputeError() ); + map< double, const SMDS_MeshNode* > theNodes; + if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes)) + return error("Circular edge is incorrectly meshed"); CNodes.clear(); - std::map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin(); + map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin(); const SMDS_MeshNode* NF = (*itn).second; CNodes.push_back( (*itn).second ); double fang = (*itn).first; - itn++; - for(; itn != theNodes.end(); itn++ ) { - CNodes.push_back( (*itn).second ); - double ang = (*itn).first - fang; - if( ang>PI ) ang = ang - 2*PI; - if( ang<-PI ) ang = ang + 2*PI; - Angles.Append( ang ); + if ( itn != theNodes.end() ) { + itn++; + for(; itn != theNodes.end(); itn++ ) { + CNodes.push_back( (*itn).second ); + double ang = (*itn).first - fang; + if( ang>PI ) ang = ang - 2*PI; + if( ang<-PI ) ang = ang + 2*PI; + Angles.Append( ang ); + } } P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() ); P0 = aCirc->Location(); - myLayerPositions.clear(); - computeLayerPositions(P0,P1); + if ( !computeLayerPositions(P0,P1)) + return false; exp.Init( CircEdge, TopAbs_VERTEX ); TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() ); @@ -267,253 +517,174 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh, double V = PC.Y() + aVec2d.Y()*myLayerPositions[i]; meshDS->SetNodeOnFace( node, faceID, U, V ); Pnts2d1.Append(gp_Pnt2d(U,V)); - Pnts2d2.Append(gp_Pnt2d(U,V)); } Nodes1[Nodes1.size()-1] = NF; Nodes2[Nodes1.size()-1] = NF; } - else if(nbe==2) { + else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL ) + { // one curve must be a half of circle and other curve must be // a segment of line - Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1); - while( !tc.IsNull() ) { - C1 = tc->BasisCurve(); - tc = Handle(Geom_TrimmedCurve)::DownCast(C1); - } - tc = Handle(Geom_TrimmedCurve)::DownCast(C2); - while( !tc.IsNull() ) { - C2 = tc->BasisCurve(); - tc = Handle(Geom_TrimmedCurve)::DownCast(C2); - } - - Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1); - Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(C2); - CircEdge = E1; - LinEdge1 = E2; - double fp = f1; - double lp = l1; - if( aCirc.IsNull() ) { - aCirc = Handle(Geom_Circle)::DownCast(C2); - CircEdge = E2; - LinEdge1 = E1; - fp = f2; - lp = l2; - aLine = Handle(Geom_Line)::DownCast(C3); - } - if( aCirc.IsNull() ) { - // not circle - return error(COMPERR_BAD_SHAPE); - } + double fp, lp; + Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp )); if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) { // not half of circle return error(COMPERR_BAD_SHAPE); } + Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 )); if( aLine.IsNull() ) { // other curve not line return error(COMPERR_BAD_SHAPE); } - SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1); - if( sm1 ) { - SMESHDS_SubMesh* sdssm1 = sm1->GetSubMeshDS(); - if( sdssm1 ) { - if( sm1->GetSubMeshDS()->NbNodes()>0 ) { - SMESH_subMesh* sm = aMesh.GetSubMesh(F); - SMESH_ComputeErrorPtr& smError = sm->GetComputeError(); - smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED, - "Invalid set of hypothesises",this)); - return false; - } - } - } - bool ok = _gen->Compute( aMesh, CircEdge, false, MeshDim_1D ); - if( !ok ) return false; - std::map< double, const SMDS_MeshNode* > theNodes; - GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes); + if ( !algo1d->ComputeCircularEdge( aMesh, CircEdge )) + return error( algo1d->GetComputeError() ); + map< double, const SMDS_MeshNode* > theNodes; + if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes) ) + return error("Circular edge is incorrectly meshed"); + map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin(); CNodes.clear(); - std::map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin(); - const SMDS_MeshNode* NF = (*itn).second; - CNodes.push_back( (*itn).second ); + CNodes.push_back( itn->second ); double fang = (*itn).first; itn++; - const SMDS_MeshNode* NL; - int nbn = 1; for(; itn != theNodes.end(); itn++ ) { - nbn++; - if( nbn == theNodes.size() ) - NL = (*itn).second; CNodes.push_back( (*itn).second ); double ang = (*itn).first - fang; if( ang>PI ) ang = ang - 2*PI; if( ang<-PI ) ang = ang + 2*PI; Angles.Append( ang ); } + const SMDS_MeshNode* NF = theNodes.begin()->second; + const SMDS_MeshNode* NL = theNodes.rbegin()->second; P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() ); gp_Pnt P2( NL->X(), NL->Y(), NL->Z() ); P0 = aCirc->Location(); - myLayerPositions.clear(); - computeLayerPositions(P0,P1); + bool linEdgeComputed; + if ( !computeLayerPositions(P0,P1,LinEdge1,&linEdgeComputed)) + return false; - gp_Vec aVec(P0,P1); - int edgeID = meshDS->ShapeToIndex(LinEdge1); - // check orientation - Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp); - gp_Pnt Ptmp; - Crv->D0(fp,Ptmp); - bool ori = true; - if( P1.Distance(Ptmp) > Precision::Confusion() ) - ori = false; - // get UV points for edge - gp_Pnt2d PF,PL; - BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL ); - PC = gp_Pnt2d( (PF.X()+PL.X())/2, (PF.Y()+PL.Y())/2 ); - gp_Vec2d V2d; - if(ori) V2d = gp_Vec2d(PC,PF); - else V2d = gp_Vec2d(PC,PL); - // add nodes on edge - double cp = (fp+lp)/2; - double dp2 = (lp-fp)/2; - NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z()); - meshDS->SetNodeOnEdge(NC, edgeID, cp); - Nodes1.resize( myLayerPositions.size()+1 ); - Nodes2.resize( myLayerPositions.size()+1 ); - int i = 0; - for(; iAddNode(P.X(), P.Y(), P.Z()); - Nodes1[i] = node; - double param; - if(ori) - param = fp + dp2*(1-myLayerPositions[i]); - else - param = cp + dp2*myLayerPositions[i]; - meshDS->SetNodeOnEdge(node, edgeID, param); - P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i], - P0.Y() - aVec.Y()*myLayerPositions[i], - P0.Z() - aVec.Z()*myLayerPositions[i] ); - node = meshDS->AddNode(P.X(), P.Y(), P.Z()); - Nodes2[i] = node; - if(!ori) - param = fp + dp2*(1-myLayerPositions[i]); - else - param = cp + dp2*myLayerPositions[i]; - meshDS->SetNodeOnEdge(node, edgeID, param); - // parameters on face - gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i], - PC.Y() + V2d.Y()*myLayerPositions[i] ); - Pnts2d1.Append(P2d); - P2d = gp_Pnt2d( PC.X() - V2d.X()*myLayerPositions[i], - PC.Y() - V2d.Y()*myLayerPositions[i] ); - Pnts2d2.Append(P2d); + if ( linEdgeComputed ) + { + if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes)) + return error("Invalid mesh on a straight edge"); + + vector< const SMDS_MeshNode* > *pNodes1 = &Nodes1, *pNodes2 = &Nodes2; + bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF ); + if ( !nodesFromP0ToP1 ) std::swap( pNodes1, pNodes2 ); + + map< double, const SMDS_MeshNode* >::reverse_iterator ritn = theNodes.rbegin(); + itn = theNodes.begin(); + for ( int i = Nodes1.size()-1; i > -1; ++itn, ++ritn, --i ) + { + (*pNodes1)[i] = ritn->second; + (*pNodes2)[i] = itn->second; + Points.Append( gpXYZ( Nodes1[i])); + Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i])); + } + NC = const_cast( itn->second ); + Points.Remove( Nodes1.size() ); } - Nodes1[ myLayerPositions.size() ] = NF; - Nodes2[ myLayerPositions.size() ] = NL; - // create 1D elements on edge - std::vector< const SMDS_MeshNode* > tmpNodes; - tmpNodes.resize(2*Nodes1.size()+1); - for(i=0; iAddEdge( tmpNodes[i-1], tmpNodes[i] ); - if(ME) meshDS->SetMeshElementOnShape(ME, edgeID); + else + { + gp_Vec aVec(P0,P1); + int edgeID = meshDS->ShapeToIndex(LinEdge1); + // check orientation + Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp); + gp_Pnt Ptmp; + Crv->D0(fp,Ptmp); + bool ori = true; + if( P1.Distance(Ptmp) > Precision::Confusion() ) + ori = false; + // get UV points for edge + gp_Pnt2d PF,PL; + BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL ); + PC = gp_Pnt2d( (PF.X()+PL.X())/2, (PF.Y()+PL.Y())/2 ); + gp_Vec2d V2d; + if(ori) V2d = gp_Vec2d(PC,PF); + else V2d = gp_Vec2d(PC,PL); + // add nodes on edge + double cp = (fp+lp)/2; + double dp2 = (lp-fp)/2; + NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z()); + meshDS->SetNodeOnEdge(NC, edgeID, cp); + Nodes1.resize( myLayerPositions.size()+1 ); + Nodes2.resize( myLayerPositions.size()+1 ); + int i = 0; + for(; iAddNode(P.X(), P.Y(), P.Z()); + Nodes1[i] = node; + double param; + if(ori) + param = fp + dp2*(1-myLayerPositions[i]); + else + param = cp + dp2*myLayerPositions[i]; + meshDS->SetNodeOnEdge(node, edgeID, param); + P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i], + P0.Y() - aVec.Y()*myLayerPositions[i], + P0.Z() - aVec.Z()*myLayerPositions[i] ); + node = meshDS->AddNode(P.X(), P.Y(), P.Z()); + Nodes2[i] = node; + if(!ori) + param = fp + dp2*(1-myLayerPositions[i]); + else + param = cp + dp2*myLayerPositions[i]; + meshDS->SetNodeOnEdge(node, edgeID, param); + // parameters on face + gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i], + PC.Y() + V2d.Y()*myLayerPositions[i] ); + Pnts2d1.Append(P2d); + } + Nodes1[ myLayerPositions.size() ] = NF; + Nodes2[ myLayerPositions.size() ] = NL; + // create 1D elements on edge + vector< const SMDS_MeshNode* > tmpNodes; + tmpNodes.resize(2*Nodes1.size()+1); + for(i=0; iAddEdge( tmpNodes[i-1], tmpNodes[i] ); + if(ME) meshDS->SetMeshElementOnShape(ME, edgeID); + } + markEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F )); } } - else { // nbe==3 + else // nbe==3 or ( nbe==2 && linEdge is INTERNAL ) + { + if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL ) + LinEdge2 = LinEdge1; + // one curve must be a part of circle and other curves must be // segments of line - Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1); - while( !tc.IsNull() ) { - C1 = tc->BasisCurve(); - tc = Handle(Geom_TrimmedCurve)::DownCast(C1); - } - tc = Handle(Geom_TrimmedCurve)::DownCast(C2); - while( !tc.IsNull() ) { - C2 = tc->BasisCurve(); - tc = Handle(Geom_TrimmedCurve)::DownCast(C2); - } - tc = Handle(Geom_TrimmedCurve)::DownCast(C3); - while( !tc.IsNull() ) { - C3 = tc->BasisCurve(); - tc = Handle(Geom_TrimmedCurve)::DownCast(C3); - } - - Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1); - Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast(C2); - Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast(C3); - CircEdge = E1; - LinEdge1 = E2; - LinEdge2 = E3; - double fp = f1; - double lp = l1; - if( aCirc.IsNull() ) { - aCirc = Handle(Geom_Circle)::DownCast(C2); - CircEdge = E2; - LinEdge1 = E3; - LinEdge2 = E1; - fp = f2; - lp = l2; - aLine1 = Handle(Geom_Line)::DownCast(C3); - aLine2 = Handle(Geom_Line)::DownCast(C1); - if( aCirc.IsNull() ) { - aCirc = Handle(Geom_Circle)::DownCast(C3); - CircEdge = E3; - LinEdge1 = E1; - LinEdge2 = E2; - fp = f3; - lp = l3; - aLine1 = Handle(Geom_Line)::DownCast(C1); - aLine2 = Handle(Geom_Line)::DownCast(C2); - } - } - if( aCirc.IsNull() ) { - // not circle - return error(COMPERR_BAD_SHAPE); - } - if( aLine1.IsNull() || aLine2.IsNull() ) { - // other curve not line + double fp, lp; + Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge )); + Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 )); + Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 )); + if( aCirc.IsNull() || aLine1.IsNull() || aLine2.IsNull() ) return error(COMPERR_BAD_SHAPE); - } - SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1); - SMESH_subMesh* sm2 = aMesh.GetSubMesh(LinEdge2); - if( sm1 && sm2 ) { - SMESHDS_SubMesh* sdssm1 = sm1->GetSubMeshDS(); - SMESHDS_SubMesh* sdssm2 = sm2->GetSubMeshDS(); - if( sdssm1 && sdssm2 ) { - if( sm1->GetSubMeshDS()->NbNodes()>0 || sm2->GetSubMeshDS()->NbNodes()>0 ) { - SMESH_subMesh* sm = aMesh.GetSubMesh(F); - SMESH_ComputeErrorPtr& smError = sm->GetComputeError(); - smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED, - "Invalid set of hypothesises",this)); - return false; - } - } - } - bool ok = _gen->Compute( aMesh, CircEdge, false, MeshDim_1D ); - if( !ok ) return false; - std::map< double, const SMDS_MeshNode* > theNodes; - GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes); + if ( !algo1d->ComputeCircularEdge( aMesh, CircEdge )) + return error( algo1d->GetComputeError() ); + map< double, const SMDS_MeshNode* > theNodes; + if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes)) + return error("Circular edge is incorrectly meshed"); + const SMDS_MeshNode* NF = theNodes.begin()->second; + const SMDS_MeshNode* NL = theNodes.rbegin()->second; CNodes.clear(); - std::map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin(); - const SMDS_MeshNode* NF = (*itn).second; - CNodes.push_back( (*itn).second ); + CNodes.push_back( NF ); + map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin(); double fang = (*itn).first; itn++; - const SMDS_MeshNode* NL; - int nbn = 1; for(; itn != theNodes.end(); itn++ ) { - nbn++; - if( nbn == theNodes.size() ) - NL = (*itn).second; CNodes.push_back( (*itn).second ); double ang = (*itn).first - fang; if( ang>PI ) ang = ang - 2*PI; @@ -524,8 +695,18 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh, gp_Pnt P2( NL->X(), NL->Y(), NL->Z() ); P0 = aCirc->Location(); - myLayerPositions.clear(); - computeLayerPositions(P0,P1); + bool linEdge1Computed, linEdge2Computed; + if ( !computeLayerPositions(P0,P1,LinEdge1,&linEdge1Computed)) + return false; + + Nodes1.resize( myLayerPositions.size()+1 ); + Nodes2.resize( myLayerPositions.size()+1 ); + + // check that both linear edges have same hypotheses + if ( !computeLayerPositions(P0,P1,LinEdge2, &linEdge2Computed)) + return false; + if ( Nodes1.size() != myLayerPositions.size()+1 ) + return error("Different hypotheses apply to radial edges"); exp.Init( LinEdge1, TopAbs_VERTEX ); TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() ); @@ -534,118 +715,169 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh, gp_Pnt PE1 = BRep_Tool::Pnt(V1); gp_Pnt PE2 = BRep_Tool::Pnt(V2); if( ( P1.Distance(PE1) > Precision::Confusion() ) && - ( P1.Distance(PE2) > Precision::Confusion() ) ) { - TopoDS_Edge E = LinEdge1; - LinEdge1 = LinEdge2; - LinEdge2 = E; + ( P1.Distance(PE2) > Precision::Confusion() ) ) + { + std::swap( LinEdge1, LinEdge2 ); + std::swap( linEdge1Computed, linEdge2Computed ); } - TopoDS_Vertex VC; + TopoDS_Vertex VC = V2; if( ( P1.Distance(PE1) > Precision::Confusion() ) && - ( P2.Distance(PE1) > Precision::Confusion() ) ) { + ( P2.Distance(PE1) > Precision::Confusion() ) ) VC = V1; - } - else VC = V2; int vertID = meshDS->ShapeToIndex(VC); + // LinEdge1 - int edgeID = meshDS->ShapeToIndex(LinEdge1); - gp_Vec aVec(P0,P1); - // check orientation - Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp); - gp_Pnt Ptmp; - Crv->D0(fp,Ptmp); - bool ori = false; - if( P1.Distance(Ptmp) > Precision::Confusion() ) - ori = true; - // get UV points for edge - gp_Pnt2d PF,PL; - BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL ); - gp_Vec2d V2d; - if(ori) { - V2d = gp_Vec2d(PF,PL); - PC = PF; - } - else { - V2d = gp_Vec2d(PL,PF); - PC = PL; - } - NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z()); - meshDS->SetNodeOnVertex(NC, vertID); - double dp = lp-fp; - Nodes1.resize( myLayerPositions.size()+1 ); - int i = 0; - for(; iAddNode(P.X(), P.Y(), P.Z()); - Nodes1[i] = node; - double param; - if(!ori) - param = fp + dp*(1-myLayerPositions[i]); - else - param = fp + dp*myLayerPositions[i]; - meshDS->SetNodeOnEdge(node, edgeID, param); - // parameters on face - gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i], - PC.Y() + V2d.Y()*myLayerPositions[i] ); - Pnts2d1.Append(P2d); + if ( linEdge1Computed ) + { + if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes)) + return error("Invalid mesh on a straight edge"); + + bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF ); + NC = const_cast + ( nodesFromP0ToP1 ? theNodes.begin()->second : theNodes.rbegin()->second ); + int i = 0, ir = Nodes1.size()-1; + int * pi = nodesFromP0ToP1 ? &i : &ir; + itn = theNodes.begin(); + if ( nodesFromP0ToP1 ) ++itn; + for ( ; i < Nodes1.size(); ++i, --ir, ++itn ) + { + Nodes1[*pi] = itn->second; + } + for ( i = 0; i < Nodes1.size()-1; ++i ) + { + Points.Append( gpXYZ( Nodes1[i])); + Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i])); + } } - Nodes1[ myLayerPositions.size() ] = NF; - // create 1D elements on edge - SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] ); - if(ME) meshDS->SetMeshElementOnShape(ME, edgeID); - for(i=1; iAddEdge( Nodes1[i-1], Nodes1[i] ); + else + { + int edgeID = meshDS->ShapeToIndex(LinEdge1); + gp_Vec aVec(P0,P1); + // check orientation + Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp); + gp_Pnt Ptmp = Crv->Value(fp); + bool ori = false; + if( P1.Distance(Ptmp) > Precision::Confusion() ) + ori = true; + // get UV points for edge + gp_Pnt2d PF,PL; + BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL ); + gp_Vec2d V2d; + if(ori) { + V2d = gp_Vec2d(PF,PL); + PC = PF; + } + else { + V2d = gp_Vec2d(PL,PF); + PC = PL; + } + NC = const_cast( VertexNode( VC, meshDS )); + if ( !NC ) + { + NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z()); + meshDS->SetNodeOnVertex(NC, vertID); + } + double dp = lp-fp; + int i = 0; + for(; iAddNode(P.X(), P.Y(), P.Z()); + Nodes1[i] = node; + double param; + if(!ori) + param = fp + dp*(1-myLayerPositions[i]); + else + param = fp + dp*myLayerPositions[i]; + meshDS->SetNodeOnEdge(node, edgeID, param); + // parameters on face + gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i], + PC.Y() + V2d.Y()*myLayerPositions[i] ); + Pnts2d1.Append(P2d); + } + Nodes1[ myLayerPositions.size() ] = NF; + // create 1D elements on edge + SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] ); if(ME) meshDS->SetMeshElementOnShape(ME, edgeID); + for(i=1; iAddEdge( Nodes1[i-1], Nodes1[i] ); + if(ME) meshDS->SetMeshElementOnShape(ME, edgeID); + } + if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL ) + Nodes2 = Nodes1; } + markEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F )); + // LinEdge2 - edgeID = meshDS->ShapeToIndex(LinEdge1); - aVec = gp_Vec(P0,P2); - // check orientation - Crv = BRep_Tool::Curve(LinEdge2,fp,lp); - Crv->D0(fp,Ptmp); - ori = false; - if( P2.Distance(Ptmp) > Precision::Confusion() ) - ori = true; - // get UV points for edge - BRep_Tool::UVPoints( LinEdge2, TopoDS::Face(aShape), PF, PL ); - if(ori) { - V2d = gp_Vec2d(PF,PL); - PC = PF; - } - else { - V2d = gp_Vec2d(PL,PF); - PC = PL; - } - dp = lp-fp; - Nodes2.resize( myLayerPositions.size()+1 ); - for(i=0; iAddNode(P.X(), P.Y(), P.Z()); - Nodes2[i] = node; - double param; - if(!ori) - param = fp + dp*(1-myLayerPositions[i]); - else - param = fp + dp*myLayerPositions[i]; - meshDS->SetNodeOnEdge(node, edgeID, param); - // parameters on face - gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i], - PC.Y() + V2d.Y()*myLayerPositions[i] ); - Pnts2d2.Append(P2d); + if ( linEdge2Computed ) + { + if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge2,true,theNodes)) + return error("Invalid mesh on a straight edge"); + + bool nodesFromP0ToP2 = ( theNodes.rbegin()->second == NL ); + int i = 0, ir = Nodes1.size()-1; + int * pi = nodesFromP0ToP2 ? &i : &ir; + itn = theNodes.begin(); + if ( nodesFromP0ToP2 ) ++itn; + for ( ; i < Nodes2.size(); ++i, --ir, ++itn ) + Nodes2[*pi] = itn->second; } - Nodes2[ myLayerPositions.size() ] = NL; - // create 1D elements on edge - ME = myHelper->AddEdge( NC, Nodes2[0] ); - if(ME) meshDS->SetMeshElementOnShape(ME, edgeID); - for(i=1; iAddEdge( Nodes2[i-1], Nodes2[i] ); + else + { + int edgeID = meshDS->ShapeToIndex(LinEdge2); + gp_Vec aVec = gp_Vec(P0,P2); + // check orientation + Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge2,fp,lp); + gp_Pnt Ptmp = Crv->Value(fp); + bool ori = false; + if( P2.Distance(Ptmp) > Precision::Confusion() ) + ori = true; + // get UV points for edge + gp_Pnt2d PF,PL; + BRep_Tool::UVPoints( LinEdge2, TopoDS::Face(aShape), PF, PL ); + gp_Vec2d V2d; + if(ori) { + V2d = gp_Vec2d(PF,PL); + PC = PF; + } + else { + V2d = gp_Vec2d(PL,PF); + PC = PL; + } + double dp = lp-fp; + for(int i=0; iAddNode(P.X(), P.Y(), P.Z()); + Nodes2[i] = node; + double param; + if(!ori) + param = fp + dp*(1-myLayerPositions[i]); + else + param = fp + dp*myLayerPositions[i]; + meshDS->SetNodeOnEdge(node, edgeID, param); + // parameters on face + gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i], + PC.Y() + V2d.Y()*myLayerPositions[i] ); + } + Nodes2[ myLayerPositions.size() ] = NL; + // create 1D elements on edge + SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes2[0] ); if(ME) meshDS->SetMeshElementOnShape(ME, edgeID); + for(int i=1; iAddEdge( Nodes2[i-1], Nodes2[i] ); + if(ME) meshDS->SetMeshElementOnShape(ME, edgeID); + } } + markEdgeAsComputedByMe( LinEdge2, aMesh.GetSubMesh( F )); } + markEdgeAsComputedByMe( CircEdge, aMesh.GetSubMesh( F )); + + // orientation + bool IsForward = ( CircEdge.Orientation()==TopAbs_FORWARD ); // create nodes and mesh elements on face // find axis of rotation @@ -658,7 +890,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh, //cout<<"Angles.Length() = "<GetLayerDistribution() ) - return error( "Invalid LayerDistribution hypothesis"); - myUsedHyps.clear(); - myUsedHyps.push_back( hyp->GetLayerDistribution() ); + myLayerPositions.clear(); - TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut ); - SMESH_Hypothesis::Hypothesis_Status aStatus; - if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus )) - return error( "StdMeshers_Regular_1D::CheckHypothesis() failed " - "with LayerDistribution hypothesis"); + SMESH_Mesh * mesh = myHelper->GetMesh(); - BRepAdaptor_Curve C3D(edge); - double f = C3D.FirstParameter(), l = C3D.LastParameter(); - list< double > params; - if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false )) - return error("StdMeshers_Regular_1D failed to compute layers distribution"); + const SMESH_Hypothesis* hyp1D = myDistributionHypo ? myDistributionHypo->GetLayerDistribution() : 0; + int nbLayers = myNbLayerHypo ? myNbLayerHypo->GetNumberOfLayers() : 0; - positions.clear(); - positions.reserve( params.size() ); - for (list::iterator itU = params.begin(); itU != params.end(); itU++) - positions.push_back( *itU / len ); - return true; + if ( !hyp1D && !nbLayers ) + { + // No own algo hypotheses assigned, so first try to find any 1D hypothesis. + // We need some edge + TopoDS_Shape edge = linEdge; + if ( edge.IsNull() && !myHelper->GetSubShape().IsNull()) + for ( TopExp_Explorer e(myHelper->GetSubShape(), TopAbs_EDGE); e.More(); e.Next()) + edge = e.Current(); + if ( !edge.IsNull() ) + { + // find a hyp usable by TNodeDistributor + SMESH_HypoFilter hypKind; + TNodeDistributor::GetDistributor(*mesh)->InitCompatibleHypoFilter(hypKind,/*ignoreAux=*/1); + hyp1D = mesh->GetHypothesis( edge, hypKind, /*fromAncestors=*/true); + } } -protected: - // ----------------------------------------------------------------------------- - TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen) - : StdMeshers_Regular_1D( hypId, studyId, gen) + if ( hyp1D ) // try to compute with hyp1D { + if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions,p1,p2,*mesh,hyp1D )) + if ( myDistributionHypo ) { // bad hyp assigned + return error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() ); + } + else { + // bad hyp found, its Ok, lets try with default nb of segnents + } } - // ----------------------------------------------------------------------------- - virtual const list & - GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool) + + if ( myLayerPositions.empty() ) // try to use nb of layers { - return myUsedHyps; + if ( !nbLayers ) + nbLayers = _gen->GetDefaultNbSegments(); + + if ( nbLayers ) + { + myLayerPositions.resize( nbLayers - 1 ); + for ( int z = 1; z < nbLayers; ++z ) + myLayerPositions[ z - 1 ] = double( z )/ double( nbLayers ); + } } - // ----------------------------------------------------------------------------- -}; -//================================================================================ -/*! - * \brief Compute positions of nodes between the internal and the external surfaces - * \retval bool - is a success - */ -//================================================================================ + // Second, check presence of a mesh built by other algo on linEdge + // and mesh conformity to my hypothesis -bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt& pIn, - const gp_Pnt& pOut) -{ - if ( myNbLayerHypo ) + bool meshComputed = (!linEdge.IsNull() && !mesh->GetSubMesh(linEdge)->IsEmpty() ); + if ( linEdgeComputed ) *linEdgeComputed = meshComputed; + + if ( meshComputed ) { - int nbSegments = myNbLayerHypo->GetNumberOfLayers(); - myLayerPositions.resize( nbSegments - 1 ); - for ( int z = 1; z < nbSegments; ++z ) - myLayerPositions[ z - 1 ] = double( z )/ double( nbSegments ); - return true; - } - if ( myDistributionHypo ) { - SMESH_Mesh * mesh = myHelper->GetMesh(); - if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions, pIn, pOut, - *mesh, myDistributionHypo )) + vector< double > nodeParams; + GetNodeParamOnEdge( mesh->GetMeshDS(), linEdge, nodeParams ); + + if ( myLayerPositions.empty() ) { - error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() ); - return false; + myLayerPositions.resize( nodeParams.size() - 2 ); + } + else if ( myDistributionHypo || myNbLayerHypo ) + { + // linEdge is computed by other algo. Check if there is a meshed face + // using nodes on linEdge + bool nodesAreUsed = false; + TopTools_ListIteratorOfListOfShape ancestIt = mesh->GetAncestors( linEdge ); + for ( ; ancestIt.More() && !nodesAreUsed; ancestIt.Next() ) + if ( ancestIt.Value().ShapeType() == TopAbs_FACE ) + nodesAreUsed = (!mesh->GetSubMesh( ancestIt.Value() )->IsEmpty()); + if ( !nodesAreUsed ) { + // rebuild them + mesh->GetSubMesh( linEdge )->ComputeStateEngine( SMESH_subMesh::CLEAN ); + if ( linEdgeComputed ) *linEdgeComputed = false; + } + else if ( myLayerPositions.size() != nodeParams.size()-2 ) { + return error("Radial edge is meshed by other algorithm"); + } } } - RETURN_BAD_RESULT("Bad hypothesis"); + + return !myLayerPositions.empty(); } @@ -827,113 +1070,82 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh, if( aShape.ShapeType() != TopAbs_FACE ) { return false; } - SMESH_subMesh * smf = aMesh.GetSubMesh(aShape); - MapShapeNbElemsItr anIt = aResMap.find(smf); - if( anIt != aResMap.end() ) { + SMESH_subMesh * sm = aMesh.GetSubMesh(aShape); + if( aResMap.count(sm) ) return false; - } - myLayerPositions.clear(); - gp_Pnt P0(0,0,0); - gp_Pnt P1(100,0,0); - computeLayerPositions(P0,P1); - - TopoDS_Edge E1,E2,E3; - Handle(Geom_Curve) C1,C2,C3; - double f1,l1,f2,l2,f3,l3; - int nbe = 0; - TopExp_Explorer exp; - for ( exp.Init( aShape, TopAbs_EDGE ); exp.More(); exp.Next() ) { - nbe++; - TopoDS_Edge E = TopoDS::Edge( exp.Current() ); - if(nbe==1) { - E1 = E; - C1 = BRep_Tool::Curve(E,f1,l1); - } - else if(nbe==2) { - E2 = E; - C2 = BRep_Tool::Curve(E,f2,l2); - } - else if(nbe==3) { - E3 = E; - C3 = BRep_Tool::Curve(E,f3,l3); - } - } + vector& aResVec = + aResMap.insert( make_pair(sm, vector(SMDSEntity_Last,0))).first->second; + + myHelper = new SMESH_MesherHelper( aMesh ); + myHelper->SetSubShape( aShape ); + auto_ptr helperDeleter( myHelper ); + + TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh); TopoDS_Edge CircEdge, LinEdge1, LinEdge2; + int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 ); + if( nbe>3 || nbe < 1 || CircEdge.IsNull() ) + return false; + + Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge )); + if( aCirc.IsNull() ) + return error(COMPERR_BAD_SHAPE); + + gp_Pnt P0 = aCirc->Location(); + gp_Pnt P1 = aCirc->Value(0.); + computeLayerPositions( P0, P1, LinEdge1 ); + int nb0d=0, nb2d_tria=0, nb2d_quad=0; - bool isQuadratic = false; - if(nbe==1) { + bool isQuadratic = false, ok = true; + if(nbe==1) + { // C1 must be a circle - Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1); - if( !aCirc.IsNull() ) { - bool ok = _gen->Evaluate( aMesh, CircEdge, aResMap ); - if(ok) { - SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge); - MapShapeNbElemsItr anIt = aResMap.find(sm); - std::vector aVec = (*anIt).second; - isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge]; - if(isQuadratic) { - // main nodes - nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size(); - // radial medium nodes - nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1); - // other medium nodes - nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size(); - } - else { - nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size(); - } - nb2d_tria = aVec[SMDSEntity_Node] + 1; - nb2d_quad = nb0d; + ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap ); + if(ok) { + const vector& aVec = aResMap[aMesh.GetSubMesh(CircEdge)]; + isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge]; + if(isQuadratic) { + // main nodes + nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size(); + // radial medium nodes + nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1); + // other medium nodes + nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size(); } + else { + nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size(); + } + nb2d_tria = aVec[SMDSEntity_Node] + 1; + nb2d_quad = nb0d; } } - else if(nbe==2) { + else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL) + { // one curve must be a half of circle and other curve must be // a segment of line - Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1); - while( !tc.IsNull() ) { - C1 = tc->BasisCurve(); - tc = Handle(Geom_TrimmedCurve)::DownCast(C1); - } - tc = Handle(Geom_TrimmedCurve)::DownCast(C2); - while( !tc.IsNull() ) { - C2 = tc->BasisCurve(); - tc = Handle(Geom_TrimmedCurve)::DownCast(C2); - } - Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1); - Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(C2); - CircEdge = E1; - LinEdge1 = E2; - double fp = f1; - double lp = l1; - if( aCirc.IsNull() ) { - aCirc = Handle(Geom_Circle)::DownCast(C2); - CircEdge = E2; - LinEdge1 = E1; - fp = f2; - lp = l2; - aLine = Handle(Geom_Line)::DownCast(C3); - } - bool ok = !aCirc.IsNull() && !aLine.IsNull(); + double fp, lp; + Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp )); if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) { // not half of circle - ok = false; + return error(COMPERR_BAD_SHAPE); + } + Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 )); + if( aLine.IsNull() ) { + // other curve not line + return error(COMPERR_BAD_SHAPE); } - SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1); - MapShapeNbElemsItr anIt = aResMap.find(sm1); - if( anIt!=aResMap.end() ) { - ok = false; + ok = !aResMap.count( aMesh.GetSubMesh(LinEdge1) ); + if ( !ok ) { + const vector& aVec = aResMap[ aMesh.GetSubMesh(LinEdge1) ]; + ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() ); } if(ok) { - ok = _gen->Evaluate( aMesh, CircEdge, aResMap ); + ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap ); } if(ok) { - SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge); - MapShapeNbElemsItr anIt = aResMap.find(sm); - std::vector aVec = (*anIt).second; - isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge]; + const vector& aVec = aResMap[ aMesh.GetSubMesh(CircEdge) ]; + isQuadratic = aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]; if(isQuadratic) { // main nodes nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size(); @@ -948,8 +1160,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh, nb2d_tria = aVec[SMDSEntity_Node] + 1; nb2d_quad = nb2d_tria * myLayerPositions.size(); // add evaluation for edges - std::vector aResVec(SMDSEntity_Last); - for(int i=SMDSEntity_Node; i aResVec(SMDSEntity_Last,0); if(isQuadratic) { aResVec[SMDSEntity_Node] = 4*myLayerPositions.size() + 3; aResVec[SMDSEntity_Quad_Edge] = 2*myLayerPositions.size() + 2; @@ -958,74 +1169,49 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh, aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1; aResVec[SMDSEntity_Edge] = 2*myLayerPositions.size() + 2; } - sm = aMesh.GetSubMesh(LinEdge1); - aResMap.insert(std::make_pair(sm,aResVec)); + aResMap[ aMesh.GetSubMesh(LinEdge1) ] = aResVec; } } - else { // nbe==3 + else // nbe==3 or ( nbe==2 && linEdge is INTERNAL ) + { + if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL ) + LinEdge2 = LinEdge1; + // one curve must be a part of circle and other curves must be // segments of line - Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1); - while( !tc.IsNull() ) { - C1 = tc->BasisCurve(); - tc = Handle(Geom_TrimmedCurve)::DownCast(C1); - } - tc = Handle(Geom_TrimmedCurve)::DownCast(C2); - while( !tc.IsNull() ) { - C2 = tc->BasisCurve(); - tc = Handle(Geom_TrimmedCurve)::DownCast(C2); - } - tc = Handle(Geom_TrimmedCurve)::DownCast(C3); - while( !tc.IsNull() ) { - C3 = tc->BasisCurve(); - tc = Handle(Geom_TrimmedCurve)::DownCast(C3); + Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 )); + Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 )); + if( aLine1.IsNull() || aLine2.IsNull() ) { + // other curve not line + return error(COMPERR_BAD_SHAPE); } - Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1); - Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast(C2); - Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast(C3); - CircEdge = E1; - LinEdge1 = E2; - LinEdge2 = E3; - double fp = f1; - double lp = l1; - if( aCirc.IsNull() ) { - aCirc = Handle(Geom_Circle)::DownCast(C2); - CircEdge = E2; - LinEdge1 = E3; - LinEdge2 = E1; - fp = f2; - lp = l2; - aLine1 = Handle(Geom_Line)::DownCast(C3); - aLine2 = Handle(Geom_Line)::DownCast(C1); - if( aCirc.IsNull() ) { - aCirc = Handle(Geom_Circle)::DownCast(C3); - CircEdge = E3; - LinEdge1 = E1; - LinEdge2 = E2; - fp = f3; - lp = l3; - aLine1 = Handle(Geom_Line)::DownCast(C1); - aLine2 = Handle(Geom_Line)::DownCast(C2); + int nbLayers = myLayerPositions.size(); + computeLayerPositions( P0, P1, LinEdge2 ); + if ( nbLayers != myLayerPositions.size() ) + return error("Different hypotheses apply to radial edges"); + + bool ok = !aResMap.count( aMesh.GetSubMesh(LinEdge1)); + if ( !ok ) { + if ( myDistributionHypo || myNbLayerHypo ) + ok = true; // override other 1d hyps + else { + const vector& aVec = aResMap[ aMesh.GetSubMesh(LinEdge1) ]; + ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() ); } } - bool ok = !aCirc.IsNull() && !aLine1.IsNull() && !aLine1.IsNull(); - SMESH_subMesh* sm = aMesh.GetSubMesh(LinEdge1); - MapShapeNbElemsItr anIt = aResMap.find(sm); - if( anIt!=aResMap.end() ) { - ok = false; - } - sm = aMesh.GetSubMesh(LinEdge2); - anIt = aResMap.find(sm); - if( anIt!=aResMap.end() ) { - ok = false; + if( ok && aResMap.count( aMesh.GetSubMesh(LinEdge2) )) { + if ( myDistributionHypo || myNbLayerHypo ) + ok = true; // override other 1d hyps + else { + const vector& aVec = aResMap[ aMesh.GetSubMesh(LinEdge2) ]; + ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() ); + } } if(ok) { - ok = _gen->Evaluate( aMesh, CircEdge, aResMap ); + ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap ); } if(ok) { - SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge); - MapShapeNbElemsItr anIt = aResMap.find(sm); - std::vector aVec = (*anIt).second; + const vector& aVec = aResMap[ aMesh.GetSubMesh(CircEdge) ]; isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge]; if(isQuadratic) { // main nodes @@ -1041,8 +1227,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh, nb2d_tria = aVec[SMDSEntity_Node] + 1; nb2d_quad = nb2d_tria * myLayerPositions.size(); // add evaluation for edges - std::vector aResVec(SMDSEntity_Last); - for(int i=SMDSEntity_Node; i aResVec(SMDSEntity_Last, 0); if(isQuadratic) { aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1; aResVec[SMDSEntity_Quad_Edge] = myLayerPositions.size() + 1; @@ -1052,17 +1237,12 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh, aResVec[SMDSEntity_Edge] = myLayerPositions.size() + 1; } sm = aMesh.GetSubMesh(LinEdge1); - aResMap.insert(std::make_pair(sm,aResVec)); + aResMap[sm] = aResVec; sm = aMesh.GetSubMesh(LinEdge2); - aResMap.insert(std::make_pair(sm,aResVec)); + aResMap[sm] = aResVec; } } - std::vector aResVec(SMDSEntity_Last); - for(int i=SMDSEntity_Node; i