src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 //  SMESH SMESH : implementaion of SMESH idl descriptions
  23 // File      : StdMeshers_RadialQuadrangle_1D2D.cxx
  24 // Module    : SMESH
  25 // Created   : Fri Oct 20 11:37:07 2006
  26 // Author    : Edward AGAPOV (eap)
  27 //
  28 #include "StdMeshers_RadialQuadrangle_1D2D.hxx"
  29
  30 //#include "StdMeshers_ProjectionUtils.hxx"
  31 #include "StdMeshers_NumberOfLayers.hxx"
  32 #include "StdMeshers_LayerDistribution.hxx"
  33 //#include "StdMeshers_Prism_3D.hxx"
  34 #include "StdMeshers_Regular_1D.hxx"
  35
  36 #include "SMDS_MeshNode.hxx"
  37 #include "SMESHDS_SubMesh.hxx"
  38 #include "SMESH_Gen.hxx"
  39 #include "SMESH_Mesh.hxx"
  40 #include "SMESH_MesherHelper.hxx"
  41 #include "SMESH_subMesh.hxx"
  42 #include "SMESH_subMeshEventListener.hxx"
  43
  44 #include "utilities.h"
  45
  46 #include <BRepAdaptor_Curve.hxx>
  47 #include <BRepBuilderAPI_MakeEdge.hxx>
  48 #include <BRep_Tool.hxx>
  49 #include <GeomAPI_ProjectPointOnSurf.hxx>
  50 #include <Geom_Circle.hxx>
  51 #include <Geom_Line.hxx>
  52 #include <Geom_TrimmedCurve.hxx>
  53 #include <TColgp_SequenceOfPnt.hxx>
  54 #include <TColgp_SequenceOfPnt2d.hxx>
  55 #include <TopExp_Explorer.hxx>
  56 #include <TopoDS.hxx>
  57
  58
  59 using namespace std;
  60
  61 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
  62 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
  63
  64 //typedef StdMeshers_ProjectionUtils TAssocTool;
  65
  66
  67 //=======================================================================
  68 //function : StdMeshers_RadialQuadrangle_1D2D
  69 //purpose  :
  70 //=======================================================================
  71
  72 StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int hypId,
  73                                                                    int studyId,
  74                                                                    SMESH_Gen* gen)
  75   :SMESH_2D_Algo(hypId, studyId, gen)
  76 {
  77   _name = "RadialQuadrangle_1D2D";
  78   _shapeType = (1 << TopAbs_FACE);        // 1 bit per shape type
  79
  80   _compatibleHypothesis.push_back("LayerDistribution2D");
  81   _compatibleHypothesis.push_back("NumberOfLayers2D");
  82   myNbLayerHypo = 0;
  83   myDistributionHypo = 0;
  84   _requireDescretBoundary = false;
  85 }
  86
  87
  88 //================================================================================
  89 /*!
  90  * \brief Destructor
  91  */
  92 //================================================================================
  93
  94 StdMeshers_RadialQuadrangle_1D2D::~StdMeshers_RadialQuadrangle_1D2D()
  95 {}
  96
  97
  98 //=======================================================================
  99 //function : CheckHypothesis
 100 //purpose  :
 101 //=======================================================================
 102
 103 bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis
 104                            (SMESH_Mesh&                          aMesh,
 105                             const TopoDS_Shape&                  aShape,
 106                             SMESH_Hypothesis::Hypothesis_Status& aStatus)
 107 {
 108   // check aShape
 109   myNbLayerHypo = 0;
 110   myDistributionHypo = 0;
 111
 112   list <const SMESHDS_Hypothesis * >::const_iterator itl;
 113
 114   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
 115   if ( hyps.size() == 0 ) {
 116     aStatus = SMESH_Hypothesis::HYP_MISSING;
 117     return false;  // can't work with no hypothesis
 118   }
 119
 120   if ( hyps.size() > 1 ) {
 121     aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
 122     return false;
 123   }
 124
 125   const SMESHDS_Hypothesis *theHyp = hyps.front();
 126
 127   string hypName = theHyp->GetName();
 128
 129   if (hypName == "NumberOfLayers2D") {
 130     myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
 131     aStatus = SMESH_Hypothesis::HYP_OK;
 132     return true;
 133   }
 134   if (hypName == "LayerDistribution2D") {
 135     myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
 136     aStatus = SMESH_Hypothesis::HYP_OK;
 137     return true;
 138   }
 139   aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 140   return true;
 141 }
 142
 143 namespace
 144 {
 145   // ------------------------------------------------------------------------------
 146   /*!
 147    * \brief Listener used to mark edges meshed by StdMeshers_RadialQuadrangle_1D2D
 148    */
 149   class TLinEdgeMarker : public SMESH_subMeshEventListener
 150   {
 151     TLinEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false) {}
 152   public:
 153     static SMESH_subMeshEventListener* getListener()
 154     {
 155       static TLinEdgeMarker theEdgeMarker;
 156       return &theEdgeMarker;
 157     }
 158   };
 159
 160   // ------------------------------------------------------------------------------
 161   /*!
 162    * \brief Mark an edge as computed by StdMeshers_RadialQuadrangle_1D2D
 163    */
 164   void markLinEdgeAsComputedByMe(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
 165   {
 166     if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
 167     {
 168       if ( !edgeSM->GetEventListenerData( TLinEdgeMarker::getListener() ))
 169         faceSubMesh->SetEventListener( TLinEdgeMarker::getListener(),
 170                                        SMESH_subMeshEventListenerData::MakeData(faceSubMesh),
 171                                        edgeSM);
 172     }
 173   }
 174   // ------------------------------------------------------------------------------
 175   /*!
 176    * \brief Return true if a radial edge was meshed with StdMeshers_RadialQuadrangle_1D2D with
 177    * the same radial distribution
 178    */
 179   bool isEdgeCompitaballyMeshed(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
 180   {
 181     if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
 182     {
 183       if ( SMESH_subMeshEventListenerData* otherFaceData =
 184            edgeSM->GetEventListenerData( TLinEdgeMarker::getListener() ))
 185       {
 186         // compare hypothesis aplied to two disk faces sharing radial edges
 187         SMESH_Mesh& mesh = *faceSubMesh->GetFather();
 188         SMESH_Algo* radialQuadAlgo = mesh.GetGen()->GetAlgo(mesh, faceSubMesh->GetSubShape() );
 189         SMESH_subMesh* otherFaceSubMesh = otherFaceData->mySubMeshes.front();
 190         const list <const SMESHDS_Hypothesis *> & hyps1 =
 191           radialQuadAlgo->GetUsedHypothesis( mesh, faceSubMesh->GetSubShape());
 192         const list <const SMESHDS_Hypothesis *> & hyps2 =
 193           radialQuadAlgo->GetUsedHypothesis( mesh, otherFaceSubMesh->GetSubShape());
 194         if( hyps1.empty() && hyps2.empty() )
 195           return true; // defaul hyps
 196         if ( hyps1.size() != hyps2.size() ||
 197              strcmp( hyps1.front()->GetName(), hyps2.front()->GetName() ))
 198           return false;
 199         ostringstream hypDump1, hypDump2;
 200         list <const SMESHDS_Hypothesis*>::const_iterator hyp1 = hyps1.begin();
 201         for ( ; hyp1 != hyps1.end(); ++hyp1 )
 202           const_cast<SMESHDS_Hypothesis*>(*hyp1)->SaveTo( hypDump1 );
 203         list <const SMESHDS_Hypothesis*>::const_iterator hyp2 = hyps2.begin();
 204         for ( ; hyp2 != hyps2.end(); ++hyp2 )
 205           const_cast<SMESHDS_Hypothesis*>(*hyp2)->SaveTo( hypDump2 );
 206         return hypDump1.str() == hypDump2.str();
 207       }
 208     }
 209     return false;
 210   }
 211
 212   //================================================================================
 213   /*!
 214    * \brief Return base curve of the edge and extremum parameters
 215    */
 216   //================================================================================
 217
 218   Handle(Geom_Curve) getCurve(const TopoDS_Edge& edge, double* f=0, double* l=0)
 219   {
 220     Handle(Geom_Curve) C;
 221     if ( !edge.IsNull() )
 222     {
 223       double first = 0., last = 0.;
 224       C = BRep_Tool::Curve(edge, first, last);
 225       if ( !C.IsNull() )
 226       {
 227         Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C);
 228         while( !tc.IsNull() ) {
 229           C = tc->BasisCurve();
 230           tc = Handle(Geom_TrimmedCurve)::DownCast(C);
 231         }
 232         if ( f ) *f = first;
 233         if ( l ) *l = last;
 234       }
 235     }
 236     return C;
 237   }
 238
 239   //================================================================================
 240   /*!
 241    * \brief Return edges of the face
 242    *  \retval int - nb of edges
 243    */
 244   //================================================================================
 245
 246   int analyseFace(const TopoDS_Shape& face,
 247                   TopoDS_Edge&        CircEdge,
 248                   TopoDS_Edge&        LinEdge1,
 249                   TopoDS_Edge&        LinEdge2)
 250   {
 251     CircEdge.Nullify(); LinEdge1.Nullify(); LinEdge2.Nullify();
 252     int nbe = 0;
 253
 254     for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next(), ++nbe )
 255     {
 256       const TopoDS_Edge& E = TopoDS::Edge( exp.Current() );
 257       double f,l;
 258       Handle(Geom_Curve) C = getCurve(E,&f,&l);
 259       if ( !C.IsNull() )
 260       {
 261         if ( C->IsKind( STANDARD_TYPE(Geom_Circle)))
 262         {
 263           if ( CircEdge.IsNull() )
 264             CircEdge = E;
 265           else
 266             return 0;
 267         }
 268         else if ( LinEdge1.IsNull() )
 269           LinEdge1 = E;
 270         else
 271           LinEdge2 = E;
 272       }
 273     }
 274     return nbe;
 275   }
 276 }
 277
 278 //=======================================================================
 279 /*!
 280  * \brief Allow algo to do something after persistent restoration
 281  * \param subMesh - restored submesh
 282  *
 283  * call markLinEdgeAsComputedByMe()
 284  */
 285 //=======================================================================
 286
 287 void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh)
 288 {
 289   if ( !faceSubMesh->IsEmpty() )
 290   {
 291     TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
 292     analyseFace( faceSubMesh->GetSubShape(), CircEdge, LinEdge1, LinEdge2 );
 293     if ( !LinEdge1.IsNull() ) markLinEdgeAsComputedByMe( LinEdge1, faceSubMesh );
 294     if ( !LinEdge2.IsNull() ) markLinEdgeAsComputedByMe( LinEdge2, faceSubMesh );
 295   }
 296 }
 297
 298 //=======================================================================
 299 //function : Compute
 300 //purpose  :
 301 //=======================================================================
 302
 303 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
 304                                                const TopoDS_Shape& aShape)
 305 {
 306   TopExp_Explorer exp;
 307   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 308
 309   myHelper = new SMESH_MesherHelper( aMesh );
 310   myHelper->IsQuadraticSubMesh( aShape );
 311   // to delete helper at exit from Compute()
 312   auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
 313
 314   myLayerPositions.clear();
 315
 316   TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
 317   int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
 318   if( nbe>3 || nbe < 1 || CircEdge.IsNull() )
 319     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)");
 320
 321   gp_Pnt P0,P1;
 322   // points for rotation
 323   TColgp_SequenceOfPnt Points;
 324   // angles for rotation
 325   TColStd_SequenceOfReal Angles;
 326   // Nodes1 and Nodes2 - nodes along radiuses
 327   // CNodes - nodes on circle edge
 328   vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes;
 329   SMDS_MeshNode * NC;
 330   // parameters edge nodes on face
 331   TColgp_SequenceOfPnt2d Pnts2d1;
 332   gp_Pnt2d PC;
 333
 334   int faceID = meshDS->ShapeToIndex(aShape);
 335   TopoDS_Face F = TopoDS::Face(aShape);
 336   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 337
 338   if(nbe==1)
 339   {
 340     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
 341
 342     bool ok = _gen->Compute( aMesh, CircEdge );
 343     if( !ok ) return false;
 344     map< double, const SMDS_MeshNode* > theNodes;
 345     ok = GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
 346     if( !ok ) return false;
 347
 348     CNodes.clear();
 349     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
 350     const SMDS_MeshNode* NF = (*itn).second;
 351     CNodes.push_back( (*itn).second );
 352     double fang = (*itn).first;
 353     if ( itn != theNodes.end() ) {
 354       itn++;
 355       for(; itn != theNodes.end(); itn++ ) {
 356         CNodes.push_back( (*itn).second );
 357         double ang = (*itn).first - fang;
 358         if( ang>PI ) ang = ang - 2*PI;
 359         if( ang<-PI ) ang = ang + 2*PI;
 360         Angles.Append( ang );
 361       }
 362     }
 363     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
 364     P0 = aCirc->Location();
 365
 366     myLayerPositions.clear();
 367     computeLayerPositions(P0,P1);
 368
 369     exp.Init( CircEdge, TopAbs_VERTEX );
 370     TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
 371     gp_Pnt2d p2dV = BRep_Tool::Parameters( V1, TopoDS::Face(aShape) );
 372
 373     NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
 374     GeomAPI_ProjectPointOnSurf PPS(P0,S);
 375     double U0,V0;
 376     PPS.Parameters(1,U0,V0);
 377     meshDS->SetNodeOnFace(NC, faceID, U0, V0);
 378     PC = gp_Pnt2d(U0,V0);
 379
 380     gp_Vec aVec(P0,P1);
 381     gp_Vec2d aVec2d(PC,p2dV);
 382     Nodes1.resize( myLayerPositions.size()+1 );
 383     Nodes2.resize( myLayerPositions.size()+1 );
 384     int i = 0;
 385     for(; i<myLayerPositions.size(); i++) {
 386       gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 387                 P0.Y() + aVec.Y()*myLayerPositions[i],
 388                 P0.Z() + aVec.Z()*myLayerPositions[i] );
 389       Points.Append(P);
 390       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 391       Nodes1[i] = node;
 392       Nodes2[i] = node;
 393       double U = PC.X() + aVec2d.X()*myLayerPositions[i];
 394       double V = PC.Y() + aVec2d.Y()*myLayerPositions[i];
 395       meshDS->SetNodeOnFace( node, faceID, U, V );
 396       Pnts2d1.Append(gp_Pnt2d(U,V));
 397     }
 398     Nodes1[Nodes1.size()-1] = NF;
 399     Nodes2[Nodes1.size()-1] = NF;
 400   }
 401   else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL )
 402   {
 403     // one curve must be a half of circle and other curve must be
 404     // a segment of line
 405     double fp, lp;
 406     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
 407     if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
 408       // not half of circle
 409       return error(COMPERR_BAD_SHAPE);
 410     }
 411     Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
 412     if( aLine.IsNull() ) {
 413       // other curve not line
 414       return error(COMPERR_BAD_SHAPE);
 415     }
 416     bool linEdgeComputed = false;
 417     if( SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1) ) {
 418       if( !sm1->IsEmpty() )
 419         if( isEdgeCompitaballyMeshed( LinEdge1, aMesh.GetSubMesh(F) ))
 420           linEdgeComputed = true;
 421         else
 422           return error("Invalid set of hypotheses");
 423     }
 424
 425     bool ok = _gen->Compute( aMesh, CircEdge );
 426     if( !ok ) return false;
 427     map< double, const SMDS_MeshNode* > theNodes;
 428     GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
 429
 430     CNodes.clear();
 431     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
 432     double fang = (*itn).first;
 433     itn++;
 434     for(; itn != theNodes.end(); itn++ ) {
 435       CNodes.push_back( (*itn).second );
 436       double ang = (*itn).first - fang;
 437       if( ang>PI ) ang = ang - 2*PI;
 438       if( ang<-PI ) ang = ang + 2*PI;
 439       Angles.Append( ang );
 440     }
 441     const SMDS_MeshNode* NF = theNodes.begin()->second;
 442     const SMDS_MeshNode* NL = theNodes.rbegin()->second;
 443     CNodes.push_back( NF );
 444     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
 445     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
 446     P0 = aCirc->Location();
 447
 448     myLayerPositions.clear();
 449     computeLayerPositions(P0,P1);
 450
 451     if ( linEdgeComputed )
 452     {
 453       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
 454         return error("Invalid mesh on a straight edge");
 455
 456       vector< const SMDS_MeshNode* > *pNodes1 = &Nodes1, *pNodes2 = &Nodes2;
 457       bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
 458       if ( !nodesFromP0ToP1 ) std::swap( pNodes1, pNodes2 );
 459
 460       map< double, const SMDS_MeshNode* >::reverse_iterator ritn = theNodes.rbegin();
 461       itn = theNodes.begin();
 462       for ( int i = Nodes1.size()-1; i > -1; ++itn, ++ritn, --i )
 463       {
 464         (*pNodes1)[i] = ritn->second;
 465         (*pNodes2)[i] =  itn->second;
 466         Points.Append( gpXYZ( Nodes1[i]));
 467         Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i]));
 468       }
 469       NC = const_cast<SMDS_MeshNode*>( itn->second );
 470       Points.Remove( Nodes1.size() );
 471     }
 472     else
 473     {
 474       gp_Vec aVec(P0,P1);
 475       int edgeID = meshDS->ShapeToIndex(LinEdge1);
 476       // check orientation
 477       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
 478       gp_Pnt Ptmp;
 479       Crv->D0(fp,Ptmp);
 480       bool ori = true;
 481       if( P1.Distance(Ptmp) > Precision::Confusion() )
 482         ori = false;
 483       // get UV points for edge
 484       gp_Pnt2d PF,PL;
 485       BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
 486       PC = gp_Pnt2d( (PF.X()+PL.X())/2, (PF.Y()+PL.Y())/2 );
 487       gp_Vec2d V2d;
 488       if(ori) V2d = gp_Vec2d(PC,PF);
 489       else V2d = gp_Vec2d(PC,PL);
 490       // add nodes on edge
 491       double cp = (fp+lp)/2;
 492       double dp2 = (lp-fp)/2;
 493       NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
 494       meshDS->SetNodeOnEdge(NC, edgeID, cp);
 495       Nodes1.resize( myLayerPositions.size()+1 );
 496       Nodes2.resize( myLayerPositions.size()+1 );
 497       int i = 0;
 498       for(; i<myLayerPositions.size(); i++) {
 499         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 500                   P0.Y() + aVec.Y()*myLayerPositions[i],
 501                   P0.Z() + aVec.Z()*myLayerPositions[i] );
 502         Points.Append(P);
 503         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 504         Nodes1[i] = node;
 505         double param;
 506         if(ori)
 507           param = fp + dp2*(1-myLayerPositions[i]);
 508         else
 509           param = cp + dp2*myLayerPositions[i];
 510         meshDS->SetNodeOnEdge(node, edgeID, param);
 511         P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i],
 512                     P0.Y() - aVec.Y()*myLayerPositions[i],
 513                     P0.Z() - aVec.Z()*myLayerPositions[i] );
 514         node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 515         Nodes2[i] = node;
 516         if(!ori)
 517           param = fp + dp2*(1-myLayerPositions[i]);
 518         else
 519           param = cp + dp2*myLayerPositions[i];
 520         meshDS->SetNodeOnEdge(node, edgeID, param);
 521         // parameters on face
 522         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
 523                       PC.Y() + V2d.Y()*myLayerPositions[i] );
 524         Pnts2d1.Append(P2d);
 525       }
 526       Nodes1[ myLayerPositions.size() ] = NF;
 527       Nodes2[ myLayerPositions.size() ] = NL;
 528       // create 1D elements on edge
 529       vector< const SMDS_MeshNode* > tmpNodes;
 530       tmpNodes.resize(2*Nodes1.size()+1);
 531       for(i=0; i<Nodes2.size(); i++)
 532         tmpNodes[Nodes2.size()-i-1] = Nodes2[i];
 533       tmpNodes[Nodes2.size()] = NC;
 534       for(i=0; i<Nodes1.size(); i++)
 535         tmpNodes[Nodes2.size()+1+i] = Nodes1[i];
 536       for(i=1; i<tmpNodes.size(); i++) {
 537         SMDS_MeshEdge* ME = myHelper->AddEdge( tmpNodes[i-1], tmpNodes[i] );
 538         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 539       }
 540       markLinEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
 541     }
 542   }
 543   else // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
 544   {
 545     if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
 546       LinEdge2 = LinEdge1;
 547
 548     // one curve must be a part of circle and other curves must be
 549     // segments of line
 550     double fp, lp;
 551     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
 552     Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
 553     Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
 554     if( aLine1.IsNull() || aLine2.IsNull() ) {
 555       // other curve not line
 556       return error(COMPERR_BAD_SHAPE);
 557     }
 558
 559     bool linEdge1Computed = false;
 560     if ( SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1))
 561       if( !sm1->IsEmpty() )
 562         if( isEdgeCompitaballyMeshed( LinEdge1, aMesh.GetSubMesh(F) ))
 563           linEdge1Computed = true;
 564         else
 565           return error("Invalid set of hypotheses");
 566
 567     bool linEdge2Computed = false;
 568     if ( SMESH_subMesh* sm2 = aMesh.GetSubMesh(LinEdge2))
 569       if( !sm2->IsEmpty() )
 570         if( isEdgeCompitaballyMeshed( LinEdge2, aMesh.GetSubMesh(F)  ))
 571           linEdge2Computed = true;
 572         else
 573           return error("Invalid set of hypotheses");
 574
 575     bool ok = _gen->Compute( aMesh, CircEdge );
 576     if( !ok ) return false;
 577     map< double, const SMDS_MeshNode* > theNodes;
 578     GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
 579
 580     const SMDS_MeshNode* NF = theNodes.begin()->second;
 581     const SMDS_MeshNode* NL = theNodes.rbegin()->second;
 582     CNodes.clear();
 583     CNodes.push_back( NF );
 584     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
 585     double fang = (*itn).first;
 586     itn++;
 587     for(; itn != theNodes.end(); itn++ ) {
 588       CNodes.push_back( (*itn).second );
 589       double ang = (*itn).first - fang;
 590       if( ang>PI ) ang = ang - 2*PI;
 591       if( ang<-PI ) ang = ang + 2*PI;
 592       Angles.Append( ang );
 593     }
 594     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
 595     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
 596     P0 = aCirc->Location();
 597
 598     myLayerPositions.clear();
 599     computeLayerPositions(P0,P1);
 600
 601     Nodes1.resize( myLayerPositions.size()+1 );
 602     Nodes2.resize( myLayerPositions.size()+1 );
 603
 604     exp.Init( LinEdge1, TopAbs_VERTEX );
 605     TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
 606     exp.Next();
 607     TopoDS_Vertex V2 = TopoDS::Vertex( exp.Current() );
 608     gp_Pnt PE1 = BRep_Tool::Pnt(V1);
 609     gp_Pnt PE2 = BRep_Tool::Pnt(V2);
 610     if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
 611         ( P1.Distance(PE2) > Precision::Confusion() ) )
 612     {
 613       std::swap( LinEdge1, LinEdge2 );
 614       std::swap( linEdge1Computed, linEdge2Computed );
 615     }
 616     TopoDS_Vertex VC = V2;
 617     if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
 618         ( P2.Distance(PE1) > Precision::Confusion() ) )
 619       VC = V1;
 620     int vertID = meshDS->ShapeToIndex(VC);
 621
 622     // LinEdge1
 623     if ( linEdge1Computed )
 624     {
 625       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
 626         return error("Invalid mesh on a straight edge");
 627
 628       bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
 629       NC = const_cast<SMDS_MeshNode*>
 630         ( nodesFromP0ToP1 ? theNodes.begin()->second : theNodes.rbegin()->second );
 631       int i = 0, ir = Nodes1.size()-1;
 632       int * pi = nodesFromP0ToP1 ? &i : &ir;
 633       itn = theNodes.begin();
 634       if ( nodesFromP0ToP1 ) ++itn;
 635       for ( ; i < Nodes1.size(); ++i, --ir, ++itn )
 636       {
 637         Nodes1[*pi] = itn->second;
 638       }
 639       for ( i = 0; i < Nodes1.size()-1; ++i )
 640       {
 641         Points.Append( gpXYZ( Nodes1[i]));
 642         Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i]));
 643       }
 644     }
 645     else
 646     {
 647       int edgeID = meshDS->ShapeToIndex(LinEdge1);
 648       gp_Vec aVec(P0,P1);
 649       // check orientation
 650       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
 651       gp_Pnt Ptmp = Crv->Value(fp);
 652       bool ori = false;
 653       if( P1.Distance(Ptmp) > Precision::Confusion() )
 654         ori = true;
 655       // get UV points for edge
 656       gp_Pnt2d PF,PL;
 657       BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
 658       gp_Vec2d V2d;
 659       if(ori) {
 660         V2d = gp_Vec2d(PF,PL);
 661         PC = PF;
 662       }
 663       else {
 664         V2d = gp_Vec2d(PL,PF);
 665         PC = PL;
 666       }
 667       NC = const_cast<SMDS_MeshNode*>( VertexNode( VC, meshDS ));
 668       if ( !NC )
 669       {
 670         NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
 671         meshDS->SetNodeOnVertex(NC, vertID);
 672       }
 673       double dp = lp-fp;
 674       int i = 0;
 675       for(; i<myLayerPositions.size(); i++) {
 676         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 677                   P0.Y() + aVec.Y()*myLayerPositions[i],
 678                   P0.Z() + aVec.Z()*myLayerPositions[i] );
 679         Points.Append(P);
 680         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 681         Nodes1[i] = node;
 682         double param;
 683         if(!ori)
 684           param = fp + dp*(1-myLayerPositions[i]);
 685         else
 686           param = fp + dp*myLayerPositions[i];
 687         meshDS->SetNodeOnEdge(node, edgeID, param);
 688         // parameters on face
 689         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
 690                       PC.Y() + V2d.Y()*myLayerPositions[i] );
 691         Pnts2d1.Append(P2d);
 692       }
 693       Nodes1[ myLayerPositions.size() ] = NF;
 694       // create 1D elements on edge
 695       SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] );
 696       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 697       for(i=1; i<Nodes1.size(); i++) {
 698         ME = myHelper->AddEdge( Nodes1[i-1], Nodes1[i] );
 699         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 700       }
 701       if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
 702         Nodes2 = Nodes1;
 703     }
 704     markLinEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
 705
 706     // LinEdge2
 707     if ( linEdge2Computed )
 708     {
 709       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge2,true,theNodes))
 710         return error("Invalid mesh on a straight edge");
 711
 712       bool nodesFromP0ToP2 = ( theNodes.rbegin()->second == NL );
 713       int i = 0, ir = Nodes1.size()-1;
 714       int * pi = nodesFromP0ToP2 ? &i : &ir;
 715       itn = theNodes.begin();
 716       if ( nodesFromP0ToP2 ) ++itn;
 717       for ( ; i < Nodes2.size(); ++i, --ir, ++itn )
 718         Nodes2[*pi] = itn->second;
 719     }
 720     else
 721     {
 722       int edgeID = meshDS->ShapeToIndex(LinEdge2);
 723       gp_Vec aVec = gp_Vec(P0,P2);
 724       // check orientation
 725       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge2,fp,lp);
 726       gp_Pnt Ptmp = Crv->Value(fp);
 727       bool ori = false;
 728       if( P2.Distance(Ptmp) > Precision::Confusion() )
 729         ori = true;
 730       // get UV points for edge
 731       gp_Pnt2d PF,PL;
 732       BRep_Tool::UVPoints( LinEdge2, TopoDS::Face(aShape), PF, PL );
 733       gp_Vec2d V2d;
 734       if(ori) {
 735         V2d = gp_Vec2d(PF,PL);
 736         PC = PF;
 737       }
 738       else {
 739         V2d = gp_Vec2d(PL,PF);
 740         PC = PL;
 741       }
 742       double dp = lp-fp;
 743       for(int i=0; i<myLayerPositions.size(); i++) {
 744         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 745                   P0.Y() + aVec.Y()*myLayerPositions[i],
 746                   P0.Z() + aVec.Z()*myLayerPositions[i] );
 747         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 748         Nodes2[i] = node;
 749         double param;
 750         if(!ori)
 751           param = fp + dp*(1-myLayerPositions[i]);
 752         else
 753           param = fp + dp*myLayerPositions[i];
 754         meshDS->SetNodeOnEdge(node, edgeID, param);
 755         // parameters on face
 756         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
 757                       PC.Y() + V2d.Y()*myLayerPositions[i] );
 758       }
 759       Nodes2[ myLayerPositions.size() ] = NL;
 760       // create 1D elements on edge
 761       SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes2[0] );
 762       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 763       for(int i=1; i<Nodes2.size(); i++) {
 764         ME = myHelper->AddEdge( Nodes2[i-1], Nodes2[i] );
 765         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 766       }
 767     }
 768     markLinEdgeAsComputedByMe( LinEdge2, aMesh.GetSubMesh( F ));
 769   }
 770
 771   // orientation
 772   bool IsForward = ( CircEdge.Orientation()==TopAbs_FORWARD );
 773
 774   // create nodes and mesh elements on face
 775   // find axis of rotation
 776   gp_Pnt P2 = gp_Pnt( CNodes[1]->X(), CNodes[1]->Y(), CNodes[1]->Z() );
 777   gp_Vec Vec1(P0,P1);
 778   gp_Vec Vec2(P0,P2);
 779   gp_Vec Axis = Vec1.Crossed(Vec2);
 780   // create elements
 781   int i = 1;
 782   //cout<<"Angles.Length() = "<<Angles.Length()<<"   Points.Length() = "<<Points.Length()<<endl;
 783   //cout<<"Nodes1.size() = "<<Nodes1.size()<<"   Pnts2d1.Length() = "<<Pnts2d1.Length()<<endl;
 784   for(; i<Angles.Length(); i++) {
 785     vector< const SMDS_MeshNode* > tmpNodes;
 786     tmpNodes.reserve(Nodes1.size());
 787     gp_Trsf aTrsf;
 788     gp_Ax1 theAxis(P0,gp_Dir(Axis));
 789     aTrsf.SetRotation( theAxis, Angles.Value(i) );
 790     gp_Trsf2d aTrsf2d;
 791     aTrsf2d.SetRotation( PC, Angles.Value(i) );
 792     // create nodes
 793     int j = 1;
 794     for(; j<=Points.Length(); j++) {
 795       double cx,cy,cz;
 796       Points.Value(j).Coord( cx, cy, cz );
 797       aTrsf.Transforms( cx, cy, cz );
 798       SMDS_MeshNode* node = myHelper->AddNode( cx, cy, cz );
 799       // find parameters on face
 800       Pnts2d1.Value(j).Coord( cx, cy );
 801       aTrsf2d.Transforms( cx, cy );
 802       // set node on face
 803       meshDS->SetNodeOnFace( node, faceID, cx, cy );
 804       tmpNodes[j-1] = node;
 805     }
 806     // create faces
 807     tmpNodes[Points.Length()] = CNodes[i];
 808     // quad
 809     for(j=0; j<Nodes1.size()-1; j++) {
 810       SMDS_MeshFace* MF;
 811       if(IsForward)
 812         MF = myHelper->AddFace( tmpNodes[j], Nodes1[j],
 813                                 Nodes1[j+1], tmpNodes[j+1] );
 814       else
 815         MF = myHelper->AddFace( tmpNodes[j], tmpNodes[j+1],
 816                                 Nodes1[j+1], Nodes1[j] );
 817       if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 818     }
 819     // tria
 820     SMDS_MeshFace* MF;
 821     if(IsForward)
 822       MF = myHelper->AddFace( NC, Nodes1[0], tmpNodes[0] );
 823     else
 824       MF = myHelper->AddFace( NC, tmpNodes[0], Nodes1[0] );
 825     if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 826     for(j=0; j<Nodes1.size(); j++) {
 827       Nodes1[j] = tmpNodes[j];
 828     }
 829   }
 830   // create last faces
 831   // quad
 832   for(i=0; i<Nodes1.size()-1; i++) {
 833     SMDS_MeshFace* MF;
 834     if(IsForward)
 835       MF = myHelper->AddFace( Nodes2[i], Nodes1[i],
 836                               Nodes1[i+1], Nodes2[i+1] );
 837     else
 838       MF = myHelper->AddFace( Nodes2[i],  Nodes2[i+1],
 839                               Nodes1[i+1], Nodes1[i] );
 840     if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 841   }
 842   // tria
 843   SMDS_MeshFace* MF;
 844   if(IsForward)
 845     MF = myHelper->AddFace( NC, Nodes1[0], Nodes2[0] );
 846   else
 847     MF = myHelper->AddFace( NC, Nodes2[0], Nodes1[0] );
 848   if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 849
 850   return true;
 851 }
 852
 853
 854 //================================================================================
 855 //================================================================================
 856 /*!
 857  * \brief Class computing layers distribution using data of
 858  *        StdMeshers_LayerDistribution hypothesis
 859  */
 860 //================================================================================
 861 //================================================================================
 862
 863 class TNodeDistributor: public StdMeshers_Regular_1D
 864 {
 865   list <const SMESHDS_Hypothesis *> myUsedHyps;
 866 public:
 867   // -----------------------------------------------------------------------------
 868   static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
 869   {
 870     const int myID = -1000;
 871     map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo;
 872     map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID );
 873     if ( id_algo == algoMap.end() )
 874       return new TNodeDistributor( myID, 0, aMesh.GetGen() );
 875     return static_cast< TNodeDistributor* >( id_algo->second );
 876   }
 877   // -----------------------------------------------------------------------------
 878   bool Compute( vector< double > &                  positions,
 879                 gp_Pnt                              pIn,
 880                 gp_Pnt                              pOut,
 881                 SMESH_Mesh&                         aMesh,
 882                 const StdMeshers_LayerDistribution* hyp)
 883   {
 884     double len = pIn.Distance( pOut );
 885     if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
 886
 887     if ( !hyp || !hyp->GetLayerDistribution() )
 888       return error( "Invalid LayerDistribution hypothesis");
 889     myUsedHyps.clear();
 890     myUsedHyps.push_back( hyp->GetLayerDistribution() );
 891
 892     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
 893     SMESH_Hypothesis::Hypothesis_Status aStatus;
 894     if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
 895       return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
 896                     "with LayerDistribution hypothesis");
 897
 898     BRepAdaptor_Curve C3D(edge);
 899     double f = C3D.FirstParameter(), l = C3D.LastParameter();
 900     list< double > params;
 901     if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
 902       return error("StdMeshers_Regular_1D failed to compute layers distribution");
 903
 904     positions.clear();
 905     positions.reserve( params.size() );
 906     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
 907       positions.push_back( *itU / len );
 908     return true;
 909   }
 910 protected:
 911   // -----------------------------------------------------------------------------
 912   TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
 913     : StdMeshers_Regular_1D( hypId, studyId, gen)
 914   {
 915   }
 916   // -----------------------------------------------------------------------------
 917   virtual const list <const SMESHDS_Hypothesis *> &
 918     GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
 919   {
 920     return myUsedHyps;
 921   }
 922   // -----------------------------------------------------------------------------
 923 };
 924
 925 //================================================================================
 926 /*!
 927  * \brief Compute positions of nodes between the internal and the external surfaces
 928   * \retval bool - is a success
 929  */
 930 //================================================================================
 931
 932 bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt& pIn,
 933                                                              const gp_Pnt& pOut)
 934 {
 935   if ( myNbLayerHypo )
 936   {
 937     int nbSegments = myNbLayerHypo->GetNumberOfLayers();
 938     myLayerPositions.resize( nbSegments - 1 );
 939     for ( int z = 1; z < nbSegments; ++z )
 940       myLayerPositions[ z - 1 ] = double( z )/ double( nbSegments );
 941     return true;
 942   }
 943   if ( myDistributionHypo ) {
 944     SMESH_Mesh * mesh = myHelper->GetMesh();
 945     if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions, pIn, pOut,
 946                                                             *mesh, myDistributionHypo ))
 947     {
 948       error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
 949       return false;
 950     }
 951   }
 952   RETURN_BAD_RESULT("Bad hypothesis");
 953 }
 954
 955
 956 //=======================================================================
 957 //function : Evaluate
 958 //purpose  :
 959 //=======================================================================
 960
 961 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
 962                                                 const TopoDS_Shape& aShape,
 963                                                 MapShapeNbElems& aResMap)
 964 {
 965   if( aShape.ShapeType() != TopAbs_FACE ) {
 966     return false;
 967   }
 968   SMESH_subMesh * smf = aMesh.GetSubMesh(aShape);
 969   MapShapeNbElemsItr anIt = aResMap.find(smf);
 970   if( anIt != aResMap.end() ) {
 971     return false;
 972   }
 973
 974   myLayerPositions.clear();
 975   gp_Pnt P0(0,0,0);
 976   gp_Pnt P1(100,0,0);
 977   computeLayerPositions(P0,P1);
 978
 979   TopoDS_Edge E1,E2,E3;
 980   Handle(Geom_Curve) C1,C2,C3;
 981   double f1,l1,f2,l2,f3,l3;
 982   int nbe = 0;
 983   TopExp_Explorer exp;
 984   for ( exp.Init( aShape, TopAbs_EDGE ); exp.More(); exp.Next() ) {
 985     nbe++;
 986     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
 987     if(nbe==1) {
 988       E1 = E;
 989       C1 = BRep_Tool::Curve(E,f1,l1);
 990     }
 991     else if(nbe==2) {
 992       E2 = E;
 993       C2 = BRep_Tool::Curve(E,f2,l2);
 994     }
 995     else if(nbe==3) {
 996       E3 = E;
 997       C3 = BRep_Tool::Curve(E,f3,l3);
 998     }
 999   }
1000
1001   TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
1002   int nb0d=0, nb2d_tria=0, nb2d_quad=0;
1003   bool isQuadratic = false;
1004   if(nbe==1) {
1005     // C1 must be a circle
1006     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
1007     if( !aCirc.IsNull() ) {
1008       bool ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
1009       if(ok) {
1010         SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
1011         MapShapeNbElemsItr anIt = aResMap.find(sm);
1012         vector<int> aVec = (*anIt).second;
1013         isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1014         if(isQuadratic) {
1015           // main nodes
1016           nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1017           // radial medium nodes
1018           nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1);
1019           // other medium nodes
1020           nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1021         }
1022         else {
1023           nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1024         }
1025         nb2d_tria = aVec[SMDSEntity_Node] + 1;
1026         nb2d_quad = nb0d;
1027       }
1028     }
1029   }
1030   else if(nbe==2) {
1031     // one curve must be a half of circle and other curve must be
1032     // a segment of line
1033     Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1034     while( !tc.IsNull() ) {
1035       C1 = tc->BasisCurve();
1036       tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1037     }
1038     tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1039     while( !tc.IsNull() ) {
1040       C2 = tc->BasisCurve();
1041       tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1042     }
1043     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
1044     Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(C2);
1045     CircEdge = E1;
1046     LinEdge1 = E2;
1047     double fp = f1;
1048     double lp = l1;
1049     if( aCirc.IsNull() ) {
1050       aCirc = Handle(Geom_Circle)::DownCast(C2);
1051       CircEdge = E2;
1052       LinEdge1 = E1;
1053       fp = f2;
1054       lp = l2;
1055       aLine = Handle(Geom_Line)::DownCast(C3);
1056     }
1057     bool ok = !aCirc.IsNull() && !aLine.IsNull();
1058     if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
1059       // not half of circle
1060       ok = false;
1061     }
1062     SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1);
1063     MapShapeNbElemsItr anIt = aResMap.find(sm1);
1064     if( anIt!=aResMap.end() ) {
1065       ok = false;
1066     }
1067     if(ok) {
1068       ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
1069     }
1070     if(ok) {
1071       SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
1072       MapShapeNbElemsItr anIt = aResMap.find(sm);
1073       vector<int> aVec = (*anIt).second;
1074       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1075       if(isQuadratic) {
1076         // main nodes
1077         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1078         // radial medium nodes
1079         nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1080         // other medium nodes
1081         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1082       }
1083       else {
1084         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1085       }
1086       nb2d_tria = aVec[SMDSEntity_Node] + 1;
1087       nb2d_quad = nb2d_tria * myLayerPositions.size();
1088       // add evaluation for edges
1089       vector<int> aResVec(SMDSEntity_Last);
1090       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1091       if(isQuadratic) {
1092         aResVec[SMDSEntity_Node] = 4*myLayerPositions.size() + 3;
1093         aResVec[SMDSEntity_Quad_Edge] = 2*myLayerPositions.size() + 2;
1094       }
1095       else {
1096         aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1097         aResVec[SMDSEntity_Edge] = 2*myLayerPositions.size() + 2;
1098       }
1099       sm = aMesh.GetSubMesh(LinEdge1);
1100       aResMap.insert(make_pair(sm,aResVec));
1101     }
1102   }
1103   else { // nbe==3
1104     // one curve must be a part of circle and other curves must be
1105     // segments of line
1106     Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1107     while( !tc.IsNull() ) {
1108       C1 = tc->BasisCurve();
1109       tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1110     }
1111     tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1112     while( !tc.IsNull() ) {
1113       C2 = tc->BasisCurve();
1114       tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1115     }
1116     tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
1117     while( !tc.IsNull() ) {
1118       C3 = tc->BasisCurve();
1119       tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
1120     }
1121     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
1122     Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast(C2);
1123     Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast(C3);
1124     CircEdge = E1;
1125     LinEdge1 = E2;
1126     LinEdge2 = E3;
1127     double fp = f1;
1128     double lp = l1;
1129     if( aCirc.IsNull() ) {
1130       aCirc = Handle(Geom_Circle)::DownCast(C2);
1131       CircEdge = E2;
1132       LinEdge1 = E3;
1133       LinEdge2 = E1;
1134       fp = f2;
1135       lp = l2;
1136       aLine1 = Handle(Geom_Line)::DownCast(C3);
1137       aLine2 = Handle(Geom_Line)::DownCast(C1);
1138       if( aCirc.IsNull() ) {
1139         aCirc = Handle(Geom_Circle)::DownCast(C3);
1140         CircEdge = E3;
1141         LinEdge1 = E1;
1142         LinEdge2 = E2;
1143         fp = f3;
1144         lp = l3;
1145         aLine1 = Handle(Geom_Line)::DownCast(C1);
1146         aLine2 = Handle(Geom_Line)::DownCast(C2);
1147       }
1148     }
1149     bool ok = !aCirc.IsNull() && !aLine1.IsNull() && !aLine1.IsNull();
1150     SMESH_subMesh* sm = aMesh.GetSubMesh(LinEdge1);
1151     MapShapeNbElemsItr anIt = aResMap.find(sm);
1152     if( anIt!=aResMap.end() ) {
1153       ok = false;
1154     }
1155     sm = aMesh.GetSubMesh(LinEdge2);
1156     anIt = aResMap.find(sm);
1157     if( anIt!=aResMap.end() ) {
1158       ok = false;
1159     }
1160     if(ok) {
1161       ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
1162     }
1163     if(ok) {
1164       SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
1165       MapShapeNbElemsItr anIt = aResMap.find(sm);
1166       vector<int> aVec = (*anIt).second;
1167       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1168       if(isQuadratic) {
1169         // main nodes
1170         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1171         // radial medium nodes
1172         nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1173         // other medium nodes
1174         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1175       }
1176       else {
1177         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1178       }
1179       nb2d_tria = aVec[SMDSEntity_Node] + 1;
1180       nb2d_quad = nb2d_tria * myLayerPositions.size();
1181       // add evaluation for edges
1182       vector<int> aResVec(SMDSEntity_Last);
1183       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1184       if(isQuadratic) {
1185         aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1186         aResVec[SMDSEntity_Quad_Edge] = myLayerPositions.size() + 1;
1187       }
1188       else {
1189         aResVec[SMDSEntity_Node] = myLayerPositions.size();
1190         aResVec[SMDSEntity_Edge] = myLayerPositions.size() + 1;
1191       }
1192       sm = aMesh.GetSubMesh(LinEdge1);
1193       aResMap.insert(make_pair(sm,aResVec));
1194       sm = aMesh.GetSubMesh(LinEdge2);
1195       aResMap.insert(make_pair(sm,aResVec));
1196     }
1197   }
1198
1199   vector<int> aResVec(SMDSEntity_Last);
1200   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1201   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1202
1203   //cout<<"nb0d = "<<nb0d<<"   nb2d_tria = "<<nb2d_tria<<"   nb2d_quad = "<<nb2d_quad<<endl;
1204   if(nb0d>0) {
1205     aResVec[0] = nb0d;
1206     if(isQuadratic) {
1207       aResVec[SMDSEntity_Quad_Triangle] = nb2d_tria;
1208       aResVec[SMDSEntity_Quad_Quadrangle] = nb2d_quad;
1209     }
1210     else {
1211       aResVec[SMDSEntity_Triangle] = nb2d_tria;
1212       aResVec[SMDSEntity_Quadrangle] = nb2d_quad;
1213     }
1214     aResMap.insert(make_pair(sm,aResVec));
1215     return true;
1216   }
1217
1218   // invalid case
1219   aResMap.insert(make_pair(sm,aResVec));
1220   SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1221   smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
1222                                         "Submesh can not be evaluated",this));
1223   return false;
1224
1225 }
1226