0020431: EDF 1020 SMESH : Radial Mesh of a cylinder

fix for the case where this algo is applied to disk segments sharing
  common straight edges

src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx

@@ -45,23 +45,15 @@
 
 #include <BRepAdaptor_Curve.hxx>
 #include <BRepBuilderAPI_MakeEdge.hxx>
-//#include <BRepTools.hxx>
 #include <BRep_Tool.hxx>
-#include <TopExp_Explorer.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
-#include <TopoDS.hxx>
-//#include <TopoDS_Shell.hxx>
-//#include <TopoDS_Solid.hxx>
-//#include <TopTools_MapOfShape.hxx>
-//#include <gp.hxx>
-//#include <gp_Pnt.hxx>
-
-
-#include <Geom_TrimmedCurve.hxx>
 #include <Geom_Circle.hxx>
 #include <Geom_Line.hxx>
+#include <Geom_TrimmedCurve.hxx>
 #include <TColgp_SequenceOfPnt.hxx>
 #include <TColgp_SequenceOfPnt2d.hxx>
-#include <GeomAPI_ProjectPointOnSurf.hxx>
+#include <TopExp_Explorer.hxx>
+#include <TopoDS.hxx>
 
 
 using namespace std;
@@ -148,6 +140,160 @@ bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis
   return true;
 }
 
+namespace
+{
+  // ------------------------------------------------------------------------------
+  /*!
+   * \brief Listener used to mark edges meshed by StdMeshers_RadialQuadrangle_1D2D
+   */
+  class TLinEdgeMarker : public SMESH_subMeshEventListener
+  {
+    TLinEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false) {}
+  public:
+    static SMESH_subMeshEventListener* getListener()
+    {
+      static TLinEdgeMarker theEdgeMarker;
+      return &theEdgeMarker;
+    }
+  };
+  
+  // ------------------------------------------------------------------------------
+  /*!
+   * \brief Mark an edge as computed by StdMeshers_RadialQuadrangle_1D2D
+   */
+  void markLinEdgeAsComputedByMe(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
+  {
+    if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
+    {
+      if ( !edgeSM->GetEventListenerData( TLinEdgeMarker::getListener() ))
+        faceSubMesh->SetEventListener( TLinEdgeMarker::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 isEdgeCompitaballyMeshed(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
+  {
+    if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
+    {
+      if ( SMESH_subMeshEventListenerData* otherFaceData =
+           edgeSM->GetEventListenerData( TLinEdgeMarker::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();
+        const list <const SMESHDS_Hypothesis *> & hyps1 =
+          radialQuadAlgo->GetUsedHypothesis( mesh, faceSubMesh->GetSubShape());
+        const list <const SMESHDS_Hypothesis *> & hyps2 =
+          radialQuadAlgo->GetUsedHypothesis( mesh, otherFaceSubMesh->GetSubShape());
+        if( hyps1.empty() && hyps2.empty() )
+          return true; // defaul hyps
+        if ( hyps1.size() != hyps2.size() ||
+             strcmp( hyps1.front()->GetName(), hyps2.front()->GetName() ))
+          return false;
+        ostringstream hypDump1, hypDump2;
+        list <const SMESHDS_Hypothesis*>::const_iterator hyp1 = hyps1.begin();
+        for ( ; hyp1 != hyps1.end(); ++hyp1 )
+          const_cast<SMESHDS_Hypothesis*>(*hyp1)->SaveTo( hypDump1 );
+        list <const SMESHDS_Hypothesis*>::const_iterator hyp2 = hyps2.begin();
+        for ( ; hyp2 != hyps2.end(); ++hyp2 )
+          const_cast<SMESHDS_Hypothesis*>(*hyp2)->SaveTo( hypDump2 );
+        return hypDump1.str() == hypDump2.str();
+      }
+    }
+    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 Allow algo to do something after persistent restoration
+ * \param subMesh - restored submesh
+ *
+ * call markLinEdgeAsComputedByMe()
+ */
+//=======================================================================
+
+void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh)
+{
+  if ( !faceSubMesh->IsEmpty() )
+  {
+    TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
+    analyseFace( faceSubMesh->GetSubShape(), CircEdge, LinEdge1, LinEdge2 );
+    if ( !LinEdge1.IsNull() ) markLinEdgeAsComputedByMe( LinEdge1, faceSubMesh );
+    if ( !LinEdge2.IsNull() ) markLinEdgeAsComputedByMe( LinEdge2, faceSubMesh );
+  }
+}
 
 //=======================================================================
 //function : Compute
@@ -162,31 +308,14 @@ 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<SMESH_MesherHelper> helperDeleter( myHelper );
 
   myLayerPositions.clear();
 
-  TopoDS_Edge E1,E2,E3;
+  TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
-  Handle(Geom_Curve) C1,C2,C3;
+  int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
-  double f1,l1,f2,l2,f3,l3;
+  if( nbe>3 || nbe < 1 || CircEdge.IsNull() )
-  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);
   
   gp_Pnt P0,P1;
@@ -196,36 +325,28 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
   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);
 
-  // orientation
+  if(nbe==1)
-  bool IsForward = F.Orientation()==TopAbs_FORWARD;
+  {
+    Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
 
-  //cout<<"RadialQuadrangle_1D2D::Compute   nbe = "<<nbe<<endl;
+    bool ok = _gen->Compute( aMesh, CircEdge );
-  TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
-  if(nbe==1) {
-    // C1 must be a circle
-    Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
-    if( aCirc.IsNull() )
-      return error(COMPERR_BAD_SHAPE);
-    CircEdge = E1;
-
-    bool ok = _gen->Compute( aMesh, CircEdge, false, MeshDim_1D );
     if( !ok ) return false;
-    std::map< double, const SMDS_MeshNode* > theNodes;
+    map< double, const SMDS_MeshNode* > theNodes;
     ok = GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
     if( !ok ) return false;
 
     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;
@@ -273,88 +394,53 @@ 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);
+    double fp, lp;
-    while( !tc.IsNull() ) {
+    Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
-      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);
-    }
     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);
+    bool linEdgeComputed = false;
-    if( sm1 ) {
+    if( SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1) ) {
-      SMESHDS_SubMesh* sdssm1 = sm1->GetSubMeshDS();
+      if( !sm1->IsEmpty() )
-      if( sdssm1 ) {
+        if( isEdgeCompitaballyMeshed( LinEdge1, aMesh.GetSubMesh(F) ))
-        if( sm1->GetSubMeshDS()->NbNodes()>0 ) {
+          linEdgeComputed = true;
-          SMESH_subMesh* sm = aMesh.GetSubMesh(F);
+        else
-          SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+          return error("Invalid set of hypotheses");
-          smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED,
-                                               "Invalid set of hypothesises",this));
-          return false;
-        }
-      }
     }
 
-    bool ok = _gen->Compute( aMesh, CircEdge, false, MeshDim_1D );
+    bool ok = _gen->Compute( aMesh, CircEdge );
     if( !ok ) return false;
-    std::map< double, const SMDS_MeshNode* > theNodes;
+    map< double, const SMDS_MeshNode* > theNodes;
     GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
 
     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++;
-    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;
+    CNodes.push_back( NF );
     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
     P0 = aCirc->Location();
@@ -362,6 +448,29 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
     myLayerPositions.clear();
     computeLayerPositions(P0,P1);
 
+    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<SMDS_MeshNode*>( itn->second );
+      Points.Remove( Nodes1.size() );
+    }
+    else
+    {
       gp_Vec aVec(P0,P1);
       int edgeID = meshDS->ShapeToIndex(LinEdge1);
       // check orientation
@@ -413,14 +522,11 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
         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);
       }
       Nodes1[ myLayerPositions.size() ] = NF;
       Nodes2[ myLayerPositions.size() ] = NL;
       // create 1D elements on edge
-    std::vector< const SMDS_MeshNode* > tmpNodes;
+      vector< const SMDS_MeshNode* > tmpNodes;
       tmpNodes.resize(2*Nodes1.size()+1);
       for(i=0; i<Nodes2.size(); i++)
         tmpNodes[Nodes2.size()-i-1] = Nodes2[i];
@@ -431,95 +537,54 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
         SMDS_MeshEdge* ME = myHelper->AddEdge( tmpNodes[i-1], tmpNodes[i] );
         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
       }
+      markLinEdgeAsComputedByMe( 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);
+    double fp, lp;
-    while( !tc.IsNull() ) {
+    Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
-      C1 = tc->BasisCurve();
+    Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
-      tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
+    Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
-    }
-    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
       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 );
+    bool linEdge1Computed = false;
+    if ( SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1))
+      if( !sm1->IsEmpty() )
+        if( isEdgeCompitaballyMeshed( LinEdge1, aMesh.GetSubMesh(F) ))
+          linEdge1Computed = true;
+        else
+          return error("Invalid set of hypotheses");
+
+    bool linEdge2Computed = false;
+    if ( SMESH_subMesh* sm2 = aMesh.GetSubMesh(LinEdge2))
+      if( !sm2->IsEmpty() )
+        if( isEdgeCompitaballyMeshed( LinEdge2, aMesh.GetSubMesh(F)  ))
+          linEdge2Computed = true;
+        else
+          return error("Invalid set of hypotheses");
+
+    bool ok = _gen->Compute( aMesh, CircEdge );
     if( !ok ) return false;
-    std::map< double, const SMDS_MeshNode* > theNodes;
+    map< double, const SMDS_MeshNode* > theNodes;
     GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
 
+    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();
+    CNodes.push_back( NF );
-    const SMDS_MeshNode* NF = (*itn).second;
+    map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
-    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;
@@ -533,6 +598,9 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
     myLayerPositions.clear();
     computeLayerPositions(P0,P1);
 
+    Nodes1.resize( myLayerPositions.size()+1 );
+    Nodes2.resize( myLayerPositions.size()+1 );
+
     exp.Init( LinEdge1, TopAbs_VERTEX );
     TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
     exp.Next();
@@ -540,25 +608,47 @@ 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() ) ) {
+        ( P1.Distance(PE2) > Precision::Confusion() ) )
-      TopoDS_Edge E = LinEdge1;
+    {
-      LinEdge1 = LinEdge2;
+      std::swap( LinEdge1, LinEdge2 );
-      LinEdge2 = E;
+      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
+    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<SMDS_MeshNode*>
+        ( 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]));
+      }
+    }
+    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;
+      gp_Pnt Ptmp = Crv->Value(fp);
-    Crv->D0(fp,Ptmp);
       bool ori = false;
       if( P1.Distance(Ptmp) > Precision::Confusion() )
         ori = true;
@@ -574,10 +664,13 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
         V2d = gp_Vec2d(PL,PF);
         PC = PL;
       }
+      NC = const_cast<SMDS_MeshNode*>( VertexNode( VC, meshDS ));
+      if ( !NC )
+      {
         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(; i<myLayerPositions.size(); i++) {
         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
@@ -602,20 +695,42 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
       SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] );
       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
       for(i=1; i<Nodes1.size(); i++) {
-      SMDS_MeshEdge* ME = myHelper->AddEdge( Nodes1[i-1], Nodes1[i] );
+        ME = myHelper->AddEdge( Nodes1[i-1], Nodes1[i] );
         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
       }
+      if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
+        Nodes2 = Nodes1;
+    }
+    markLinEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
+
     // LinEdge2
-    edgeID = meshDS->ShapeToIndex(LinEdge2);
+    if ( linEdge2Computed )
-    aVec = gp_Vec(P0,P2);
+    {
+      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;
+    }
+    else
+    {
+      int edgeID = meshDS->ShapeToIndex(LinEdge2);
+      gp_Vec aVec = gp_Vec(P0,P2);
       // check orientation
-    Crv = BRep_Tool::Curve(LinEdge2,fp,lp);
+      Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge2,fp,lp);
-    Crv->D0(fp,Ptmp);
+      gp_Pnt Ptmp = Crv->Value(fp);
-    ori = false;
+      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;
@@ -624,9 +739,8 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
         V2d = gp_Vec2d(PL,PF);
         PC = PL;
       }
-    dp = lp-fp;
+      double dp = lp-fp;
-    Nodes2.resize( myLayerPositions.size()+1 );
+      for(int i=0; i<myLayerPositions.size(); i++) {
-    for(i=0; i<myLayerPositions.size(); i++) {
         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
                   P0.Y() + aVec.Y()*myLayerPositions[i],
                   P0.Z() + aVec.Z()*myLayerPositions[i] );
@@ -641,17 +755,21 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
         // parameters on face
         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
                       PC.Y() + V2d.Y()*myLayerPositions[i] );
-      Pnts2d2.Append(P2d);
       }
       Nodes2[ myLayerPositions.size() ] = NL;
       // create 1D elements on edge
-    ME = myHelper->AddEdge( NC, Nodes2[0] );
+      SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes2[0] );
       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
-    for(i=1; i<Nodes2.size(); i++) {
+      for(int i=1; i<Nodes2.size(); i++) {
-      SMDS_MeshEdge* ME = myHelper->AddEdge( Nodes2[i-1], Nodes2[i] );
+        ME = myHelper->AddEdge( Nodes2[i-1], Nodes2[i] );
         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
       }
     }
+    markLinEdgeAsComputedByMe( LinEdge2, aMesh.GetSubMesh( F ));
+  }
+
+  // orientation
+  bool IsForward = ( CircEdge.Orientation()==TopAbs_FORWARD );
 
   // create nodes and mesh elements on face
   // find axis of rotation
@@ -664,7 +782,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
   //cout<<"Angles.Length() = "<<Angles.Length()<<"   Points.Length() = "<<Points.Length()<<endl;
   //cout<<"Nodes1.size() = "<<Nodes1.size()<<"   Pnts2d1.Length() = "<<Pnts2d1.Length()<<endl;
   for(; i<Angles.Length(); i++) {
-    std::vector< const SMDS_MeshNode* > tmpNodes;
+    vector< const SMDS_MeshNode* > tmpNodes;
     tmpNodes.reserve(Nodes1.size());
     gp_Trsf aTrsf;
     gp_Ax1 theAxis(P0,gp_Dir(Axis));
@@ -729,10 +847,6 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
     MF = myHelper->AddFace( NC, Nodes2[0], Nodes1[0] );
   if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 
-
-  // to delete helper at exit from Compute()
-  std::auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
-
   return true;
 }
 
@@ -895,7 +1009,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
       if(ok) {
         SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
         MapShapeNbElemsItr anIt = aResMap.find(sm);
-        std::vector<int> aVec = (*anIt).second;
+        vector<int> aVec = (*anIt).second;
         isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
         if(isQuadratic) {
           // main nodes
@@ -956,7 +1070,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
     if(ok) {
       SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
       MapShapeNbElemsItr anIt = aResMap.find(sm);
-      std::vector<int> aVec = (*anIt).second;
+      vector<int> aVec = (*anIt).second;
       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
       if(isQuadratic) {
         // main nodes
@@ -972,7 +1086,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<int> aResVec(SMDSEntity_Last);
+      vector<int> aResVec(SMDSEntity_Last);
       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
       if(isQuadratic) {
         aResVec[SMDSEntity_Node] = 4*myLayerPositions.size() + 3;
@@ -983,7 +1097,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
         aResVec[SMDSEntity_Edge] = 2*myLayerPositions.size() + 2;
       }
       sm = aMesh.GetSubMesh(LinEdge1);
-      aResMap.insert(std::make_pair(sm,aResVec));
+      aResMap.insert(make_pair(sm,aResVec));
     }
   }
   else { // nbe==3
@@ -1049,7 +1163,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
     if(ok) {
       SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
       MapShapeNbElemsItr anIt = aResMap.find(sm);
-      std::vector<int> aVec = (*anIt).second;
+      vector<int> aVec = (*anIt).second;
       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
       if(isQuadratic) {
         // main nodes
@@ -1065,7 +1179,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<int> aResVec(SMDSEntity_Last);
+      vector<int> aResVec(SMDSEntity_Last);
       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
       if(isQuadratic) {
         aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
@@ -1076,13 +1190,13 @@ 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.insert(make_pair(sm,aResVec));
       sm = aMesh.GetSubMesh(LinEdge2);
-      aResMap.insert(std::make_pair(sm,aResVec));
+      aResMap.insert(make_pair(sm,aResVec));
     }
   }
 
-  std::vector<int> aResVec(SMDSEntity_Last);
+  vector<int> aResVec(SMDSEntity_Last);
   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 
@@ -1097,12 +1211,12 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
       aResVec[SMDSEntity_Triangle] = nb2d_tria;
       aResVec[SMDSEntity_Quadrangle] = nb2d_quad;
     }
-    aResMap.insert(std::make_pair(sm,aResVec));
+    aResMap.insert(make_pair(sm,aResVec));
     return true;
   }
 
   // invalid case
-  aResMap.insert(std::make_pair(sm,aResVec));
+  aResMap.insert(make_pair(sm,aResVec));
   SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
   smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
                                         "Submesh can not be evaluated",this));