src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx

   1 // Copyright (C) 2007-2015  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, or (at your option) any later version.
  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 // File      : StdMeshers_QuadFromMedialAxis_1D2D.cxx
  23 // Created   : Wed Jun  3 17:33:45 2015
  24 // Author    : Edward AGAPOV (eap)
  25
  26 #include "StdMeshers_QuadFromMedialAxis_1D2D.hxx"
  27
  28 #include "SMESH_Block.hxx"
  29 #include "SMESH_Gen.hxx"
  30 #include "SMESH_MAT2d.hxx"
  31 #include "SMESH_Mesh.hxx"
  32 #include "SMESH_MeshEditor.hxx"
  33 #include "SMESH_MesherHelper.hxx"
  34 #include "SMESH_ProxyMesh.hxx"
  35 #include "SMESH_subMesh.hxx"
  36 #include "StdMeshers_FaceSide.hxx"
  37 #include "StdMeshers_Regular_1D.hxx"
  38 #include "StdMeshers_ViscousLayers2D.hxx"
  39
  40 #include <BRepBuilderAPI_MakeEdge.hxx>
  41 #include <BRepTools.hxx>
  42 #include <BRep_Tool.hxx>
  43 #include <GeomAPI_Interpolate.hxx>
  44 #include <Geom_Surface.hxx>
  45 #include <Precision.hxx>
  46 #include <TColgp_HArray1OfPnt.hxx>
  47 #include <TopExp.hxx>
  48 #include <TopLoc_Location.hxx>
  49 #include <TopTools_MapOfShape.hxx>
  50 #include <TopoDS.hxx>
  51 #include <TopoDS_Edge.hxx>
  52 #include <TopoDS_Face.hxx>
  53 #include <TopoDS_Vertex.hxx>
  54 #include <gp_Pnt.hxx>
  55
  56 #include <list>
  57 #include <vector>
  58
  59 //================================================================================
  60 /*!
  61  * \brief 1D algo
  62  */
  63 class StdMeshers_QuadFromMedialAxis_1D2D::Algo1D : public StdMeshers_Regular_1D
  64 {
  65 public:
  66   Algo1D(int studyId, SMESH_Gen* gen):
  67     StdMeshers_Regular_1D( gen->GetANewId(), studyId, gen )
  68   {
  69   }
  70   void SetSegmentLength( double len )
  71   {
  72     _value[ BEG_LENGTH_IND ] = len;
  73     _value[ PRECISION_IND  ] = 1e-7;
  74     _hypType = LOCAL_LENGTH;
  75   }
  76 };
  77
  78 //================================================================================
  79 /*!
  80  * \brief Constructor sets algo features
  81  */
  82 //================================================================================
  83
  84 StdMeshers_QuadFromMedialAxis_1D2D::StdMeshers_QuadFromMedialAxis_1D2D(int        hypId,
  85                                                                        int        studyId,
  86                                                                        SMESH_Gen* gen)
  87   : StdMeshers_Quadrangle_2D(hypId, studyId, gen),
  88     _regular1D( 0 )
  89 {
  90   _name = "QuadFromMedialAxis_1D2D";
  91   _shapeType = (1 << TopAbs_FACE);
  92   _onlyUnaryInput          = true;  // FACE by FACE so far
  93   _requireDiscreteBoundary = false; // make 1D by myself
  94   _supportSubmeshes        = true; // make 1D by myself
  95   _neededLowerHyps[ 1 ]    = true;  // suppress warning on hiding a global 1D algo
  96   _neededLowerHyps[ 2 ]    = true;  // suppress warning on hiding a global 2D algo
  97   _compatibleHypothesis.clear();
  98   _compatibleHypothesis.push_back("ViscousLayers2D");
  99 }
 100
 101 //================================================================================
 102 /*!
 103  * \brief Destructor
 104  */
 105 //================================================================================
 106
 107 StdMeshers_QuadFromMedialAxis_1D2D::~StdMeshers_QuadFromMedialAxis_1D2D()
 108 {
 109   delete _regular1D;
 110   _regular1D = 0;
 111 }
 112
 113 //================================================================================
 114 /*!
 115  * \brief Check if needed hypotheses are present
 116  */
 117 //================================================================================
 118
 119 bool StdMeshers_QuadFromMedialAxis_1D2D::CheckHypothesis(SMESH_Mesh&         aMesh,
 120                                                          const TopoDS_Shape& aShape,
 121                                                          Hypothesis_Status&  aStatus)
 122 {
 123   aStatus = HYP_OK;
 124   return true; // does not require hypothesis
 125 }
 126
 127 namespace
 128 {
 129   typedef map< const SMDS_MeshNode*, list< const SMDS_MeshNode* > > TMergeMap;
 130
 131   //================================================================================
 132   /*!
 133    * \brief Sinuous face
 134    */
 135   struct SinuousFace
 136   {
 137     FaceQuadStruct::Ptr   _quad;
 138     vector< TopoDS_Edge > _edges;
 139     vector< TopoDS_Edge > _sinuSide[2], _shortSide[2];
 140     vector< TopoDS_Edge > _sinuEdges;
 141     int                   _nbWires;
 142     list< int >           _nbEdgesInWire;
 143     TMergeMap             _nodesToMerge;
 144
 145     SinuousFace( const TopoDS_Face& f ): _quad( new FaceQuadStruct )
 146     {
 147       list< TopoDS_Edge > edges;
 148       _nbWires = SMESH_Block::GetOrderedEdges (f, edges, _nbEdgesInWire);
 149       _edges.assign( edges.begin(), edges.end() );
 150
 151       _quad->side.resize( 4 );
 152       _quad->face = f;
 153     }
 154     const TopoDS_Face& Face() const { return _quad->face; }
 155   };
 156
 157   //================================================================================
 158   /*!
 159    * \brief Temporary mesh
 160    */
 161   struct TmpMesh : public SMESH_Mesh
 162   {
 163     TmpMesh()
 164     {
 165       _myMeshDS = new SMESHDS_Mesh(/*id=*/0, /*isEmbeddedMode=*/true);
 166     }
 167   };
 168
 169   //================================================================================
 170   /*!
 171    * \brief Return a member of a std::pair
 172    */
 173   //================================================================================
 174
 175   template< typename T >
 176   T& get( std::pair< T, T >& thePair, bool is2nd )
 177   {
 178     return is2nd ? thePair.second : thePair.first;
 179   }
 180
 181   //================================================================================
 182   /*!
 183    * \brief Select two EDGEs from a map, either mapped to least values or to max values
 184    */
 185   //================================================================================
 186
 187   // template< class TVal2EdgesMap >
 188   // void getTwo( bool                 least,
 189   //              TVal2EdgesMap&       map,
 190   //              vector<TopoDS_Edge>& twoEdges,
 191   //              vector<TopoDS_Edge>& otherEdges)
 192   // {
 193   //   twoEdges.clear();
 194   //   otherEdges.clear();
 195   //   if ( least )
 196   //   {
 197   //     TVal2EdgesMap::iterator i = map.begin();
 198   //     twoEdges.push_back( i->second );
 199   //     twoEdges.push_back( ++i->second );
 200   //     for ( ; i != map.end(); ++i )
 201   //       otherEdges.push_back( i->second );
 202   //   }
 203   //   else
 204   //   {
 205   //     TVal2EdgesMap::reverse_iterator i = map.rbegin();
 206   //     twoEdges.push_back( i->second );
 207   //     twoEdges.push_back( ++i->second );
 208   //     for ( ; i != map.rend(); ++i )
 209   //       otherEdges.push_back( i->second );
 210   //   }
 211   //   TopoDS_Vertex v;
 212   //   if ( TopExp::CommonVertex( twoEdges[0], twoEdges[1], v ))
 213   //   {
 214   //     twoEdges.clear(); // two EDGEs must not be connected
 215   //     otherEdges.clear();
 216   //   }
 217   // }
 218
 219   //================================================================================
 220   /*!
 221    * \brief Finds out a minimal segment length given EDGEs will be divided into.
 222    *        This length is further used to discretize the Medial Axis
 223    */
 224   //================================================================================
 225
 226   double getMinSegLen(SMESH_MesherHelper&        theHelper,
 227                       const vector<TopoDS_Edge>& theEdges)
 228   {
 229     TmpMesh tmpMesh;
 230     SMESH_Mesh* mesh = theHelper.GetMesh();
 231
 232     vector< SMESH_Algo* > algos( theEdges.size() );
 233     for ( size_t i = 0; i < theEdges.size(); ++i )
 234     {
 235       SMESH_subMesh* sm = mesh->GetSubMesh( theEdges[i] );
 236       algos[i] = sm->GetAlgo();
 237     }
 238
 239     const int nbSegDflt = mesh->GetGen()->GetDefaultNbSegments();
 240     double minSegLen    = Precision::Infinite();
 241
 242     for ( size_t i = 0; i < theEdges.size(); ++i )
 243     {
 244       SMESH_subMesh* sm = mesh->GetSubMesh( theEdges[i] );
 245       if ( SMESH_Algo::IsStraight( theEdges[i], /*degenResult=*/true ))
 246         continue;
 247       // get algo
 248       size_t iOpp = ( theEdges.size() == 4 ? (i+2)%4 : i );
 249       SMESH_Algo*  algo = sm->GetAlgo();
 250       if ( !algo ) algo = algos[ iOpp ];
 251       // get hypo
 252       SMESH_Hypothesis::Hypothesis_Status status = SMESH_Hypothesis::HYP_MISSING;
 253       if ( algo )
 254       {
 255         if ( !algo->CheckHypothesis( *mesh, theEdges[i], status ))
 256           algo->CheckHypothesis( *mesh, theEdges[iOpp], status );
 257       }
 258       // compute
 259       if ( status != SMESH_Hypothesis::HYP_OK )
 260       {
 261         minSegLen = Min( minSegLen, SMESH_Algo::EdgeLength( theEdges[i] ) / nbSegDflt );
 262       }
 263       else
 264       {
 265         tmpMesh.Clear();
 266         tmpMesh.ShapeToMesh( TopoDS_Shape());
 267         tmpMesh.ShapeToMesh( theEdges[i] );
 268         try {
 269           mesh->GetGen()->Compute( tmpMesh, theEdges[i], true, true ); // make nodes on VERTEXes
 270           if ( !algo->Compute( tmpMesh, theEdges[i] ))
 271             continue;
 272         }
 273         catch (...) {
 274           continue;
 275         }
 276         SMDS_EdgeIteratorPtr segIt = tmpMesh.GetMeshDS()->edgesIterator();
 277         while ( segIt->more() )
 278         {
 279           const SMDS_MeshElement* seg = segIt->next();
 280           double len = SMESH_TNodeXYZ( seg->GetNode(0) ).Distance( seg->GetNode(1) );
 281           minSegLen = Min( minSegLen, len );
 282         }
 283       }
 284     }
 285     if ( Precision::IsInfinite( minSegLen ))
 286       minSegLen = mesh->GetShapeDiagonalSize() / nbSegDflt;
 287
 288     return minSegLen;
 289   }
 290
 291   //================================================================================
 292   /*!
 293    * \brief Returns EDGEs located between two VERTEXes at which given MA branches end
 294    *  \param [in] br1 - one MA branch
 295    *  \param [in] br2 - one more MA branch
 296    *  \param [in] allEdges - all EDGEs of a FACE
 297    *  \param [out] shortEdges - the found EDGEs
 298    *  \return bool - is OK or not
 299    */
 300   //================================================================================
 301
 302   bool getConnectedEdges( const SMESH_MAT2d::Branch* br1,
 303                           const SMESH_MAT2d::Branch* br2,
 304                           const vector<TopoDS_Edge>& allEdges,
 305                           vector<TopoDS_Edge>&       shortEdges)
 306   {
 307     vector< size_t > edgeIDs[4];
 308     br1->getGeomEdges( edgeIDs[0], edgeIDs[1] );
 309     br2->getGeomEdges( edgeIDs[2], edgeIDs[3] );
 310
 311     // EDGEs returned by a Branch form a connected chain with a VERTEX where
 312     // the Branch ends at the chain middle. One of end EDGEs of the chain is common
 313     // with either end EDGE of the chain of the other Branch, or the chains are connected
 314     // at a common VERTEX;
 315
 316     // Get indices of end EDGEs of the branches
 317     bool vAtStart1 = ( br1->getEnd(0)->_type == SMESH_MAT2d::BE_ON_VERTEX );
 318     bool vAtStart2 = ( br2->getEnd(0)->_type == SMESH_MAT2d::BE_ON_VERTEX );
 319     size_t iEnd[4] = {
 320       vAtStart1 ? edgeIDs[0].back() : edgeIDs[0][0],
 321       vAtStart1 ? edgeIDs[1].back() : edgeIDs[1][0],
 322       vAtStart2 ? edgeIDs[2].back() : edgeIDs[2][0],
 323       vAtStart2 ? edgeIDs[3].back() : edgeIDs[3][0]
 324     };
 325
 326     set< size_t > connectedIDs;
 327     TopoDS_Vertex vCommon;
 328     // look for the same EDGEs
 329     for ( int i = 0; i < 2; ++i )
 330       for ( int j = 2; j < 4; ++j )
 331         if ( iEnd[i] == iEnd[j] )
 332         {
 333           connectedIDs.insert( edgeIDs[i].begin(), edgeIDs[i].end() );
 334           connectedIDs.insert( edgeIDs[j].begin(), edgeIDs[j].end() );
 335           i = j = 4;
 336         }
 337     if ( connectedIDs.empty() )
 338       // look for connected EDGEs
 339       for ( int i = 0; i < 2; ++i )
 340         for ( int j = 2; j < 4; ++j )
 341           if ( TopExp::CommonVertex( allEdges[ iEnd[i]], allEdges[ iEnd[j]], vCommon ))
 342           {
 343             connectedIDs.insert( edgeIDs[i].begin(), edgeIDs[i].end() );
 344             connectedIDs.insert( edgeIDs[j].begin(), edgeIDs[j].end() );
 345             i = j = 4;
 346           }
 347     if ( connectedIDs.empty() ||                     // nothing
 348          allEdges.size() - connectedIDs.size() < 2 ) // too many
 349       return false;
 350
 351     // set shortEdges in the order as in allEdges
 352     if ( connectedIDs.count( 0 ) &&
 353          connectedIDs.count( allEdges.size()-1 ))
 354     {
 355       size_t iE = allEdges.size()-1;
 356       while ( connectedIDs.count( iE-1 ))
 357         --iE;
 358       for ( size_t i = 0; i < connectedIDs.size(); ++i )
 359       {
 360         shortEdges.push_back( allEdges[ iE ]);
 361         iE = ( iE + 1 ) % allEdges.size();
 362       }
 363     }
 364     else
 365     {
 366       set< size_t >::iterator i = connectedIDs.begin();
 367       for ( ; i != connectedIDs.end(); ++i )
 368         shortEdges.push_back( allEdges[ *i ]);
 369     }
 370     return true;
 371   }
 372
 373   //================================================================================
 374   /*!
 375    * \brief Find EDGEs to discretize using projection from MA
 376    *  \param [in,out] theSinuFace - the FACE to be meshed
 377    *  \return bool - OK or not
 378    *
 379    * It separates all EDGEs into four sides of a quadrangle connected in the order:
 380    * theSinuEdges[0], theShortEdges[0], theSinuEdges[1], theShortEdges[1]
 381    */
 382   //================================================================================
 383
 384   bool getSinuousEdges( SMESH_MesherHelper& theHelper,
 385                         SinuousFace&        theSinuFace)
 386   {
 387     vector<TopoDS_Edge> * theSinuEdges  = & theSinuFace._sinuSide [0];
 388     vector<TopoDS_Edge> * theShortEdges = & theSinuFace._shortSide[0];
 389     theSinuEdges[0].clear();
 390     theSinuEdges[1].clear();
 391     theShortEdges[0].clear();
 392     theShortEdges[1].clear();
 393
 394     vector<TopoDS_Edge> & allEdges = theSinuFace._edges;
 395     const size_t nbEdges = allEdges.size();
 396     if ( nbEdges < 4 && theSinuFace._nbWires == 1 )
 397       return false;
 398
 399     if ( theSinuFace._nbWires == 2 ) // ring
 400     {
 401       size_t nbOutEdges = theSinuFace._nbEdgesInWire.front();
 402       theSinuEdges[0].assign ( allEdges.begin(), allEdges.begin() + nbOutEdges );
 403       theSinuEdges[1].assign ( allEdges.begin() + nbOutEdges, allEdges.end() );
 404       return true;
 405     }
 406     if ( theSinuFace._nbWires > 2 )
 407       return false;
 408
 409     // create MedialAxis to find short edges by analyzing MA branches
 410     double minSegLen = getMinSegLen( theHelper, allEdges );
 411     SMESH_MAT2d::MedialAxis ma( theSinuFace.Face(), allEdges, minSegLen * 3 );
 412
 413     // in an initial request case, theFace represents a part of a river with almost parallel banks
 414     // so there should be two branch points
 415     using SMESH_MAT2d::BranchEnd;
 416     using SMESH_MAT2d::Branch;
 417     const vector< const BranchEnd* >& braPoints = ma.getBranchPoints();
 418     if ( braPoints.size() < 2 )
 419       return false;
 420     TopTools_MapOfShape shortMap;
 421     size_t nbBranchPoints = 0;
 422     for ( size_t i = 0; i < braPoints.size(); ++i )
 423     {
 424       vector< const Branch* > vertBranches; // branches with an end on VERTEX
 425       for ( size_t ib = 0; ib < braPoints[i]->_branches.size(); ++ib )
 426       {
 427         const Branch* branch = braPoints[i]->_branches[ ib ];
 428         if ( branch->hasEndOfType( SMESH_MAT2d::BE_ON_VERTEX ))
 429           vertBranches.push_back( branch );
 430       }
 431       if ( vertBranches.size() != 2 || braPoints[i]->_branches.size() != 3)
 432         continue;
 433
 434       // get common EDGEs of two branches
 435       if ( !getConnectedEdges( vertBranches[0], vertBranches[1],
 436                                allEdges, theShortEdges[ nbBranchPoints > 0 ] ))
 437         return false;
 438
 439       for ( size_t iS = 0; iS < theShortEdges[ nbBranchPoints ].size(); ++iS )
 440         shortMap.Add( theShortEdges[ nbBranchPoints ][ iS ]);
 441
 442       ++nbBranchPoints;
 443     }
 444
 445     if ( nbBranchPoints != 2 )
 446       return false;
 447
 448     // add to theSinuEdges all edges that are not theShortEdges
 449     vector< vector<TopoDS_Edge> > sinuEdges(1);
 450     TopoDS_Vertex vCommon;
 451     for ( size_t i = 0; i < allEdges.size(); ++i )
 452     {
 453       if ( !shortMap.Contains( allEdges[i] ))
 454       {
 455         if ( !sinuEdges.back().empty() )
 456           if ( !TopExp::CommonVertex( sinuEdges.back().back(), allEdges[ i ], vCommon ))
 457             sinuEdges.resize( sinuEdges.size() + 1 );
 458
 459         sinuEdges.back().push_back( allEdges[i] );
 460       }
 461     }
 462     if ( sinuEdges.size() == 3 )
 463     {
 464       if ( !TopExp::CommonVertex( sinuEdges.back().back(), sinuEdges[0][0], vCommon ))
 465         return false;
 466       vector<TopoDS_Edge>& last = sinuEdges.back();
 467       last.insert( last.end(), sinuEdges[0].begin(), sinuEdges[0].end() );
 468       sinuEdges[0].swap( last );
 469       sinuEdges.resize( 2 );
 470     }
 471     if ( sinuEdges.size() != 2 )
 472       return false;
 473
 474     theSinuEdges[0].swap( sinuEdges[0] );
 475     theSinuEdges[1].swap( sinuEdges[1] );
 476
 477     if ( !TopExp::CommonVertex( theSinuEdges[0].back(), theShortEdges[0][0], vCommon ) ||
 478          !vCommon.IsSame( theHelper.IthVertex( 1, theSinuEdges[0].back() )))
 479       theShortEdges[0].swap( theShortEdges[1] );
 480
 481     theSinuFace._sinuEdges = theSinuEdges[0];
 482     theSinuFace._sinuEdges.insert( theSinuFace._sinuEdges.end(),
 483                                    theSinuEdges[1].begin(), theSinuEdges[1].end() );
 484
 485     return ( theShortEdges[0].size() > 0 && theShortEdges[1].size() > 0 &&
 486              theSinuEdges [0].size() > 0 && theSinuEdges [1].size() > 0 );
 487
 488     // the sinuous EDGEs can be composite and C0 continuous,
 489     // therefor we use a complex criterion to find TWO short non-sinuous EDGEs
 490     // and the rest EDGEs will be treated as sinuous.
 491     // A short edge should have the following features:
 492     // a) straight
 493     // b) short
 494     // c) with convex corners at ends
 495     // d) far from the other short EDGE
 496
 497     // vector< double > isStraightEdge( nbEdges, 0 ); // criterion value
 498
 499     // // a0) evaluate continuity
 500     // const double contiWgt = 0.5; // weight of continuity in the criterion
 501     // multimap< int, TopoDS_Edge > continuity;
 502     // for ( size_t i = 0; i < nbEdges; ++I )
 503     // {
 504     //   BRepAdaptor_Curve curve( allEdges[i] );
 505     //   GeomAbs_Shape C = GeomAbs_CN;
 506     //   try:
 507     //     C = curve.Continuity(); // C0, G1, C1, G2, C2, C3, CN
 508     //   catch ( Standard_Failure ) {}
 509     //   continuity.insert( make_pair( C, allEdges[i] ));
 510     //   isStraight[i] += double( C ) / double( CN ) * contiWgt;
 511     // }
 512
 513     // // try to choose by continuity
 514     // int mostStraight = (int) continuity.rbegin()->first;
 515     // int lessStraight = (int) continuity.begin()->first;
 516     // if ( mostStraight != lessStraight )
 517     // {
 518     //   int nbStraight = continuity.count( mostStraight );
 519     //   if ( nbStraight == 2 )
 520     //   {
 521     //     getTwo( /*least=*/false, continuity, theShortEdges, theSinuEdges );
 522     //   }
 523     //   else if ( nbStraight == 3 && nbEdges == 4 )
 524     //   {
 525     //     theSinuEdges.push_back( continuity.begin()->second );
 526     //     vector<TopoDS_Edge>::iterator it =
 527     //       std::find( allEdges.begin(), allEdges.end(), theSinuEdges[0] );
 528     //     int i = std::distance( allEdges.begin(), it );
 529     //     theSinuEdges .push_back( allEdges[( i+2 )%4 ]);
 530     //     theShortEdges.push_back( allEdges[( i+1 )%4 ]);
 531     //     theShortEdges.push_back( allEdges[( i+3 )%4 ]);
 532     //   }
 533     //   if ( theShortEdges.size() == 2 )
 534     //     return true;
 535     // }
 536
 537     // // a) curvature; evaluate aspect ratio
 538     // {
 539     //   const double curvWgt = 0.5;
 540     //   for ( size_t i = 0; i < nbEdges; ++I )
 541     //   {
 542     //     BRepAdaptor_Curve curve( allEdges[i] );
 543     //     double curvature = 1;
 544     //     if ( !curve.IsClosed() )
 545     //     {
 546     //       const double f = curve.FirstParameter(), l = curve.LastParameter();
 547     //       gp_Pnt pf = curve.Value( f ), pl = curve.Value( l );
 548     //       gp_Lin line( pf, pl.XYZ() - pf.XYZ() );
 549     //       double distMax = 0;
 550     //       for ( double u = f; u < l; u += (l-f)/30. )
 551     //         distMax = Max( distMax, line.SquareDistance( curve.Value( u )));
 552     //       curvature = Sqrt( distMax ) / ( pf.Distance( pl ));
 553     //     }
 554     //     isStraight[i] += curvWgt / (              curvature + 1e-20 );
 555     //   }
 556     // }
 557     // // b) length
 558     // {
 559     //   const double lenWgt = 0.5;
 560     //   for ( size_t i = 0; i < nbEdges; ++I )
 561     //   {
 562     //     double length = SMESH_Algo::Length( allEdges[i] );
 563     //     if ( length > 0 )
 564     //       isStraight[i] += lenWgt / length;
 565     //   }
 566     // }
 567     // // c) with convex corners at ends
 568     // {
 569     //   const double cornerWgt = 0.25;
 570     //   for ( size_t i = 0; i < nbEdges; ++I )
 571     //   {
 572     //     double convex = 0;
 573     //     int iPrev = SMESH_MesherHelper::WrapIndex( int(i)-1, nbEdges );
 574     //     int iNext = SMESH_MesherHelper::WrapIndex( int(i)+1, nbEdges );
 575     //     TopoDS_Vertex v = helper.IthVertex( 0, allEdges[i] );
 576     //     double angle = SMESH_MesherHelper::GetAngle( allEdges[iPrev], allEdges[i], theFace, v );
 577     //     if ( angle < M_PI ) // [-PI; PI]
 578     //       convex += ( angle + M_PI ) / M_PI / M_PI;
 579     //     v = helper.IthVertex( 1, allEdges[i] );
 580     //     angle = SMESH_MesherHelper::GetAngle( allEdges[iNext], allEdges[i], theFace, v );
 581     //     if ( angle < M_PI ) // [-PI; PI]
 582     //       convex += ( angle + M_PI ) / M_PI / M_PI;
 583     //     isStraight[i] += cornerWgt * convex;
 584     //   }
 585     // }
 586   }
 587
 588   //================================================================================
 589   /*!
 590    * \brief Creates an EDGE from a sole branch of MA
 591    */
 592   //================================================================================
 593
 594   TopoDS_Edge makeEdgeFromMA( SMESH_MesherHelper&            theHelper,
 595                               const SMESH_MAT2d::MedialAxis& theMA,
 596                               const double                   theMinSegLen)
 597   {
 598     if ( theMA.nbBranches() != 1 )
 599       return TopoDS_Edge();
 600
 601     vector< gp_XY > uv;
 602     theMA.getPoints( theMA.getBranch(0), uv );
 603     if ( uv.size() < 2 )
 604       return TopoDS_Edge();
 605
 606     TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() );
 607     Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
 608
 609     vector< gp_Pnt > pnt;
 610     pnt.reserve( uv.size() * 2 );
 611     pnt.push_back( surface->Value( uv[0].X(), uv[0].Y() ));
 612     for ( size_t i = 1; i < uv.size(); ++i )
 613     {
 614       gp_Pnt p = surface->Value( uv[i].X(), uv[i].Y() );
 615       int nbDiv = int( p.Distance( pnt.back() ) / theMinSegLen );
 616       for ( int iD = 1; iD < nbDiv; ++iD )
 617       {
 618         double  R = iD / double( nbDiv );
 619         gp_XY uvR = uv[i-1] * (1 - R) + uv[i] * R;
 620         pnt.push_back( surface->Value( uvR.X(), uvR.Y() ));
 621       }
 622       pnt.push_back( p );
 623     }
 624
 625     // cout << "from salome.geom import geomBuilder" << endl;
 626     // cout << "geompy = geomBuilder.New(salome.myStudy)" << endl;
 627     Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, pnt.size());
 628     for ( size_t i = 0; i < pnt.size(); ++i )
 629     {
 630       gp_Pnt& p = pnt[i];
 631       points->SetValue( i+1, p );
 632       // cout << "geompy.MakeVertex( "<< p.X()<<", " << p.Y()<<", " << p.Z()
 633       //      <<" theName = 'p_" << i << "')" << endl;
 634     }
 635
 636     GeomAPI_Interpolate interpol( points, /*isClosed=*/false, gp::Resolution());
 637     interpol.Perform();
 638     if ( !interpol.IsDone())
 639       return TopoDS_Edge();
 640
 641     TopoDS_Edge branchEdge = BRepBuilderAPI_MakeEdge(interpol.Curve());
 642     return branchEdge;
 643   }
 644
 645   //================================================================================
 646   /*!
 647    * \brief Returns a type of shape, to which a hypothesis used to mesh a given edge is assigned
 648    */
 649   //================================================================================
 650
 651   TopAbs_ShapeEnum getHypShape( SMESH_Mesh* mesh, const TopoDS_Shape& edge )
 652   {
 653     TopAbs_ShapeEnum shapeType = TopAbs_SHAPE;
 654
 655     SMESH_subMesh* sm = mesh->GetSubMesh( edge );
 656     SMESH_Algo*  algo = sm->GetAlgo();
 657     if ( !algo ) return shapeType;
 658
 659     const list <const SMESHDS_Hypothesis *> & hyps =
 660       algo->GetUsedHypothesis( *mesh, edge, /*ignoreAuxiliary=*/true );
 661     if ( hyps.empty() ) return shapeType;
 662
 663     TopoDS_Shape shapeOfHyp =
 664       SMESH_MesherHelper::GetShapeOfHypothesis( hyps.front(), edge, mesh);
 665
 666     return SMESH_MesherHelper::GetGroupType( shapeOfHyp, /*woCompound=*/true);
 667   }
 668
 669   //================================================================================
 670   /*!
 671    * \brief Discretize a sole branch of MA an returns parameters of divisions on MA
 672    */
 673   //================================================================================
 674
 675   bool divideMA( SMESH_MesherHelper&            theHelper,
 676                  const SMESH_MAT2d::MedialAxis& theMA,
 677                  const SinuousFace&             theSinuFace,
 678                  SMESH_Algo*                    the1dAlgo,
 679                  const double                   theMinSegLen,
 680                  vector<double>&                theMAParams )
 681   {
 682     // check if all EDGEs of one size are meshed, then MA discretization is not needed
 683     SMESH_Mesh* mesh = theHelper.GetMesh();
 684     size_t nbComputedEdges[2] = { 0, 0 };
 685     for ( size_t iS = 0; iS < 2; ++iS )
 686       for ( size_t i = 0; i < theSinuFace._sinuSide[iS].size(); ++i )
 687       {
 688         bool isComputed = ( ! mesh->GetSubMesh( theSinuFace._sinuSide[iS][i] )->IsEmpty() );
 689         nbComputedEdges[ iS ] += isComputed;
 690       }
 691     if ( nbComputedEdges[0] == theSinuFace._sinuSide[0].size() ||
 692          nbComputedEdges[1] == theSinuFace._sinuSide[1].size() )
 693       return true; // discretization is not needed
 694
 695
 696     TopoDS_Edge branchEdge = makeEdgeFromMA( theHelper, theMA, theMinSegLen );
 697     if ( branchEdge.IsNull() )
 698       return false;
 699
 700     // const char* file = "/misc/dn25/salome/eap/salome/misc/tmp/MAedge.brep";
 701     // BRepTools::Write( branchEdge, file);
 702     // cout << "Write " << file << endl;
 703
 704     // look for a most local hyps assigned to theSinuEdges
 705     TopoDS_Edge edge = theSinuFace._sinuEdges[0];
 706     int mostSimpleShape = (int) getHypShape( mesh, edge );
 707     for ( size_t i = 1; i < theSinuFace._sinuEdges.size(); ++i )
 708     {
 709       int shapeType = (int) getHypShape( mesh, theSinuFace._sinuEdges[i] );
 710       if ( shapeType > mostSimpleShape )
 711         edge = theSinuFace._sinuEdges[i];
 712     }
 713
 714     SMESH_Algo* algo = the1dAlgo;
 715     if ( mostSimpleShape != TopAbs_SHAPE )
 716     {
 717       algo = mesh->GetSubMesh( edge )->GetAlgo();
 718       SMESH_Hypothesis::Hypothesis_Status status;
 719       if ( !algo->CheckHypothesis( *mesh, edge, status ))
 720         algo = the1dAlgo;
 721     }
 722
 723     TmpMesh tmpMesh;
 724     tmpMesh.ShapeToMesh( branchEdge );
 725     try {
 726       mesh->GetGen()->Compute( tmpMesh, branchEdge, true, true ); // make nodes on VERTEXes
 727       if ( !algo->Compute( tmpMesh, branchEdge ))
 728         return false;
 729     }
 730     catch (...) {
 731       return false;
 732     }
 733     return SMESH_Algo::GetNodeParamOnEdge( tmpMesh.GetMeshDS(), branchEdge, theMAParams );
 734   }
 735
 736   //================================================================================
 737   /*!
 738    * \brief Select division parameters on MA and make them coincide at ends with
 739    *        projections of VERTEXes to MA for a given pair of opposite EDGEs
 740    *  \param [in] theEdgePairInd - index of the EDGE pair
 741    *  \param [in] theDivPoints - the BranchPoint's dividing MA into parts each
 742    *         corresponding to a unique pair of opposite EDGEs
 743    *  \param [in] theMAParams - the MA division parameters
 744    *  \param [out] theSelectedMAParams - the selected MA parameters
 745    *  \return bool - is OK
 746    */
 747   //================================================================================
 748
 749   bool getParamsForEdgePair( const size_t                              theEdgePairInd,
 750                              const vector< SMESH_MAT2d::BranchPoint >& theDivPoints,
 751                              const vector<double>&                     theMAParams,
 752                              vector<double>&                           theSelectedMAParams)
 753   {
 754     if ( theDivPoints.empty() )
 755     {
 756       theSelectedMAParams = theMAParams;
 757       return true;
 758     }
 759     if ( theEdgePairInd > theDivPoints.size() || theMAParams.empty() )
 760       return false;
 761
 762     // find a range of params to copy
 763
 764     double par1 = 0;
 765     size_t iPar1 = 0;
 766     if ( theEdgePairInd > 0 )
 767     {
 768       const SMESH_MAT2d::BranchPoint& bp = theDivPoints[ theEdgePairInd-1 ];
 769       bp._branch->getParameter( bp, par1 );
 770       while ( theMAParams[ iPar1 ] < par1 ) ++iPar1;
 771       if ( par1 - theMAParams[ iPar1-1 ] < theMAParams[ iPar1 ] - par1 )
 772         --iPar1;
 773     }
 774
 775     double par2 = 1;
 776     size_t iPar2 = theMAParams.size() - 1;
 777     if ( theEdgePairInd < theDivPoints.size() )
 778     {
 779       const SMESH_MAT2d::BranchPoint& bp = theDivPoints[ theEdgePairInd ];
 780       bp._branch->getParameter( bp, par2 );
 781       iPar2 = iPar1;
 782       while ( theMAParams[ iPar2 ] < par2 ) ++iPar2;
 783       if ( par2 - theMAParams[ iPar2-1 ] < theMAParams[ iPar2 ] - par2 )
 784         --iPar2;
 785     }
 786
 787     theSelectedMAParams.assign( theMAParams.begin() + iPar1,
 788                                 theMAParams.begin() + iPar2 + 1 );
 789
 790     // adjust theSelectedMAParams to fit between par1 and par2
 791
 792     double d = par1 - theSelectedMAParams[0];
 793     double f = ( par2 - par1 ) / ( theSelectedMAParams.back() - theSelectedMAParams[0] );
 794
 795     for ( size_t i = 0; i < theSelectedMAParams.size(); ++i )
 796     {
 797       theSelectedMAParams[i] += d;
 798       theSelectedMAParams[i] = par1 + ( theSelectedMAParams[i] - par1 ) * f;
 799     }
 800
 801     return true;
 802   }
 803
 804   //--------------------------------------------------------------------------------
 805   // node or node parameter on EDGE
 806   struct NodePoint
 807   {
 808     const SMDS_MeshNode* _node;
 809     double               _u;
 810     int                  _edgeInd; // index in theSinuEdges vector
 811
 812     NodePoint(): _node(0), _u(0), _edgeInd(-1) {}
 813     NodePoint(const SMDS_MeshNode* n, double u, size_t iEdge ): _node(n), _u(u), _edgeInd(iEdge) {}
 814     NodePoint(double u, size_t iEdge) : _node(0), _u(u), _edgeInd(iEdge) {}
 815     NodePoint(const SMESH_MAT2d::BoundaryPoint& p) : _node(0), _u(p._param), _edgeInd(p._edgeIndex) {}
 816     gp_Pnt Point(const vector< Handle(Geom_Curve) >& curves) const
 817     {
 818       return curves[ _edgeInd ]->Value( _u );
 819     }
 820   };
 821
 822   //================================================================================
 823   /*!
 824    * \brief Finds a VERTEX corresponding to a point on EDGE, which is also filled
 825    *        with a node on the VERTEX, present or created
 826    *  \param [in,out] theNodePnt - the node position on the EDGE
 827    *  \param [in] theSinuEdges - the sinuous EDGEs
 828    *  \param [in] theMeshDS - the mesh
 829    *  \return bool - true if the \a theBndPnt is on VERTEX
 830    */
 831   //================================================================================
 832
 833   bool findVertex( NodePoint&                  theNodePnt,
 834                    const vector<TopoDS_Edge>&  theSinuEdges,
 835                    size_t                      theEdgeIndPrev,
 836                    size_t                      theEdgeIndNext,
 837                    SMESHDS_Mesh*               theMeshDS)
 838   {
 839     if ( theNodePnt._edgeInd >= theSinuEdges.size() )
 840       return false;
 841
 842     double f,l;
 843     BRep_Tool::Range( theSinuEdges[ theNodePnt._edgeInd ], f,l );
 844     const double tol = 1e-3 * ( l - f );
 845
 846     TopoDS_Vertex V;
 847     if      ( Abs( f - theNodePnt._u ) < tol )
 848       V = SMESH_MesherHelper::IthVertex( 0, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false);
 849     else if ( Abs( l - theNodePnt._u ) < tol )
 850       V = SMESH_MesherHelper::IthVertex( 1, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false);
 851     else if ( theEdgeIndPrev != theEdgeIndNext )
 852       TopExp::CommonVertex( theSinuEdges[theEdgeIndPrev], theSinuEdges[theEdgeIndNext], V );
 853
 854     if ( !V.IsNull() )
 855     {
 856       theNodePnt._node = SMESH_Algo::VertexNode( V, theMeshDS );
 857       if ( !theNodePnt._node )
 858       {
 859         gp_Pnt p = BRep_Tool::Pnt( V );
 860         theNodePnt._node = theMeshDS->AddNode( p.X(), p.Y(), p.Z() );
 861         theMeshDS->SetNodeOnVertex( theNodePnt._node, V );
 862       }
 863       return true;
 864     }
 865     return false;
 866   }
 867
 868   //================================================================================
 869   /*!
 870    * \brief Add to the map of NodePoint's those on VERTEXes
 871    *  \param [in,out] theHelper - the helper
 872    *  \param [in] theMA - Medial Axis
 873    *  \param [in] theMinSegLen - minimal segment length
 874    *  \param [in] theDivPoints - projections of VERTEXes to MA
 875    *  \param [in] theSinuEdges - the sinuous EDGEs
 876    *  \param [in] theSideEdgeIDs - indices of sinuous EDGEs per side
 877    *  \param [in] theIsEdgeComputed - is sinuous EGDE is meshed
 878    *  \param [in,out] thePointsOnE - the map to fill
 879    *  \param [out] theNodes2Merge - the map of nodes to merge
 880    */
 881   //================================================================================
 882
 883   bool projectVertices( SMESH_MesherHelper&                           theHelper,
 884                         //const double                                  theMinSegLen,
 885                         const SMESH_MAT2d::MedialAxis&                theMA,
 886                         const vector< SMESH_MAT2d::BranchPoint >&     theDivPoints,
 887                         const vector< std::size_t > &                 theEdgeIDs1,
 888                         const vector< std::size_t > &                 theEdgeIDs2,
 889                         const vector<TopoDS_Edge>&                    theSinuEdges,
 890                         const vector< Handle(Geom_Curve) >&           theCurves,
 891                         const vector< bool >&                         theIsEdgeComputed,
 892                         map< double, pair< NodePoint, NodePoint > > & thePointsOnE,
 893                         TMergeMap&                                    theNodes2Merge)
 894   {
 895     if ( theDivPoints.empty() )
 896       return true;
 897
 898     SMESHDS_Mesh* meshDS = theHelper.GetMeshDS();
 899
 900     double uMA;
 901     SMESH_MAT2d::BoundaryPoint bp[2];
 902     const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0);
 903
 904     // fill a map holding NodePoint's of ends of theSinuEdges
 905     map< double, pair< NodePoint, NodePoint > > extremaNP;
 906     map< double, pair< NodePoint, NodePoint > >::iterator u2NP0, u2NP1;
 907     if ( !branch.getBoundaryPoints( 0., bp[0], bp[1] ) ||
 908          !theMA.getBoundary().moveToClosestEdgeEnd( bp[0] ) ||
 909          !theMA.getBoundary().moveToClosestEdgeEnd( bp[1] )) return false;
 910     u2NP0 = extremaNP.insert
 911       ( make_pair( 0., make_pair( NodePoint( bp[0]), NodePoint( bp[1])))).first;
 912     if ( !branch.getBoundaryPoints( 1., bp[0], bp[1] ) ||
 913          !theMA.getBoundary().moveToClosestEdgeEnd( bp[0] ) ||
 914          !theMA.getBoundary().moveToClosestEdgeEnd( bp[1] )) return false;
 915     u2NP1 = extremaNP.insert
 916       ( make_pair( 1., make_pair( NodePoint( bp[0]), NodePoint( bp[1])))).first;
 917
 918     // project theDivPoints
 919     for ( size_t i = 0; i < theDivPoints.size(); ++i )
 920     {
 921       if ( !branch.getParameter( theDivPoints[i], uMA ))
 922         return false;
 923       if ( !branch.getBoundaryPoints( theDivPoints[i], bp[0], bp[1] ))
 924         return false;
 925
 926       NodePoint  np[2] = {
 927         NodePoint( bp[0] ),
 928         NodePoint( bp[1] )
 929       };
 930       bool isVertex[2] = {
 931         findVertex( np[0], theSinuEdges, theEdgeIDs1[i], theEdgeIDs1[i+1], meshDS ),
 932         findVertex( np[1], theSinuEdges, theEdgeIDs2[i], theEdgeIDs2[i+1], meshDS )
 933       };
 934
 935       map< double, pair< NodePoint, NodePoint > >::iterator u2NP =
 936         thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1]))).first;
 937
 938       if ( !isVertex[0] && !isVertex[1] ) return false; // error
 939       if ( isVertex[0] && isVertex[1] )
 940         continue;
 941       const size_t iVert = isVertex[0] ? 0 : 1;
 942       const size_t iNode   = 1 - iVert;
 943
 944       bool isOppComputed = theIsEdgeComputed[ np[ iNode ]._edgeInd ];
 945       if ( !isOppComputed )
 946         continue;
 947
 948       // a VERTEX is projected on a meshed EDGE; there are two options:
 949       // 1) a projected point is joined with a closet node if a strip between this and neighbor
 950       // projection is WIDE enough; joining is done by creating a node coincident with the
 951       //  existing node which will be merged together after all;
 952       // 2) a neighbor projection is merged with this one if it is TOO CLOSE; a node of deleted
 953       // projection is set to the BoundaryPoint of this projection
 954
 955       // evaluate distance to neighbor projections
 956       const double rShort = 0.2;
 957       bool isShortPrev[2], isShortNext[2];
 958       map< double, pair< NodePoint, NodePoint > >::iterator u2NPPrev = u2NP, u2NPNext = u2NP;
 959       --u2NPPrev; ++u2NPNext;
 960       bool hasPrev = ( u2NP     != thePointsOnE.begin() );
 961       bool hasNext = ( u2NPNext != thePointsOnE.end() );
 962       if ( !hasPrev ) u2NPPrev = u2NP0;
 963       if ( !hasNext ) u2NPNext = u2NP1;
 964       for ( int iS = 0; iS < 2; ++iS ) // side with Vertex and side with Nodes
 965       {
 966         NodePoint np     = get( u2NP->second,     iS );
 967         NodePoint npPrev = get( u2NPPrev->second, iS );
 968         NodePoint npNext = get( u2NPNext->second, iS );
 969         gp_Pnt     p     = np    .Point( theCurves );
 970         gp_Pnt     pPrev = npPrev.Point( theCurves );
 971         gp_Pnt     pNext = npNext.Point( theCurves );
 972         double  distPrev = p.Distance( pPrev );
 973         double  distNext = p.Distance( pNext );
 974         double         r = distPrev / ( distPrev + distNext );
 975         isShortPrev[iS] = ( r < rShort );
 976         isShortNext[iS] = (( 1 - r ) > ( 1 - rShort ));
 977       }
 978       // if ( !hasPrev ) isShortPrev[0] = isShortPrev[1] = false;
 979       // if ( !hasNext ) isShortNext[0] = isShortNext[1] = false;
 980
 981       map< double, pair< NodePoint, NodePoint > >::iterator u2NPClose;
 982
 983       if (( isShortPrev[0] && isShortPrev[1] ) || // option 2) -> remove a too close projection
 984           ( isShortNext[0] && isShortNext[1] ))
 985       {
 986         u2NPClose = isShortPrev[0] ? u2NPPrev : u2NPNext;
 987         NodePoint& npProj  = get( u2NP->second,      iNode ); // NP of VERTEX projection
 988         NodePoint npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj
 989         NodePoint npCloseV = get( u2NPClose->second, iVert ); // NP close to VERTEX
 990         if ( !npCloseV._node )
 991         {
 992           npProj = npCloseN;
 993           thePointsOnE.erase( isShortPrev[0] ? u2NPPrev : u2NPNext );
 994           continue;
 995         }
 996         else
 997         {
 998           // can't remove the neighbor projection as it is also from VERTEX, -> option 1)
 999         }
1000       }
1001       // else: option 1) - wide enough -> "duplicate" existing node
1002       {
1003         u2NPClose = isShortPrev[ iNode ] ? u2NPPrev : u2NPNext;
1004         NodePoint& npProj   = get( u2NP->second,      iNode ); // NP of VERTEX projection
1005         NodePoint& npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj
1006         // npProj._edgeInd = npCloseN._edgeInd;
1007         // npProj._u       = npCloseN._u + 1e-3 * Abs( get( u2NPPrev->second, iNode )._u -
1008         //                                             get( u2NPNext->second, iNode )._u );
1009         gp_Pnt        p = npProj.Point( theCurves );
1010         npProj._node    = meshDS->AddNode( p.X(), p.Y(), p.Z() );
1011         meshDS->SetNodeOnEdge( npProj._node, theSinuEdges[ npProj._edgeInd ], npProj._u  );
1012
1013         theNodes2Merge[ npCloseN._node ].push_back( npProj._node );
1014       }
1015     }
1016     return true;
1017   }
1018
1019   //================================================================================
1020   /*!
1021    * \brief Divide the sinuous EDGEs by projecting the division point of Medial
1022    *        Axis to the EGDEs
1023    *  \param [in] theHelper - the helper
1024    *  \param [in] theMinSegLen - minimal segment length
1025    *  \param [in] theMA - the Medial Axis
1026    *  \param [in] theMAParams - parameters of division points of \a theMA
1027    *  \param [in] theSinuEdges - the EDGEs to make nodes on
1028    *  \param [in] theSinuSide0Size - the number of EDGEs in the 1st sinuous side
1029    *  \return bool - is OK or not
1030    */
1031   //================================================================================
1032
1033   bool computeSinuEdges( SMESH_MesherHelper&        theHelper,
1034                          double                     /*theMinSegLen*/,
1035                          SMESH_MAT2d::MedialAxis&   theMA,
1036                          vector<double>&            theMAParams,
1037                          SinuousFace&               theSinuFace)
1038   {
1039     if ( theMA.nbBranches() != 1 )
1040       return false;
1041
1042     // normalize theMAParams
1043     for ( size_t i = 0; i < theMAParams.size(); ++i )
1044       theMAParams[i] /= theMAParams.back();
1045
1046
1047     SMESH_Mesh*     mesh = theHelper.GetMesh();
1048     SMESHDS_Mesh* meshDS = theHelper.GetMeshDS();
1049     double f,l;
1050
1051     const vector< TopoDS_Edge >& theSinuEdges = theSinuFace._sinuEdges;
1052     vector< Handle(Geom_Curve) > curves ( theSinuEdges.size() );
1053     vector< int >                edgeIDs( theSinuEdges.size() );
1054     vector< bool >            isComputed( theSinuEdges.size() );
1055     //bool hasComputed = false;
1056     for ( size_t i = 0; i < theSinuEdges.size(); ++i )
1057     {
1058       curves[i] = BRep_Tool::Curve( theSinuEdges[i], f,l );
1059       if ( !curves[i] )
1060         return false;
1061       SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] );
1062       edgeIDs   [i] = sm->GetId();
1063       isComputed[i] = ( !sm->IsEmpty() );
1064       if ( isComputed[i] )
1065       {
1066         TopAbs_ShapeEnum shape = getHypShape( mesh, theSinuEdges[i] );
1067         if ( shape == TopAbs_SHAPE || shape <= TopAbs_FACE )
1068         {
1069           // EDGE computed using global hypothesis -> clear it
1070           bool hasComputedFace = false;
1071           PShapeIteratorPtr faceIt = theHelper.GetAncestors( theSinuEdges[i], *mesh, TopAbs_FACE );
1072           while ( const TopoDS_Shape* face = faceIt->next() )
1073             if (( !face->IsSame( theSinuFace.Face())) &&
1074                 ( hasComputedFace = !mesh->GetSubMesh( *face )->IsEmpty() ))
1075               break;
1076           if ( !hasComputedFace )
1077             sm->ComputeStateEngine( SMESH_subMesh::CLEAN );
1078           isComputed[i] = false;
1079         }
1080       }
1081     }
1082
1083     const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0);
1084     SMESH_MAT2d::BoundaryPoint bp[2];
1085
1086     vector< std::size_t > edgeIDs1, edgeIDs2;
1087     vector< SMESH_MAT2d::BranchPoint > divPoints;
1088     branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints );
1089     for ( size_t i = 0; i < edgeIDs1.size(); ++i )
1090       if ( isComputed[ edgeIDs1[i]] &&
1091            isComputed[ edgeIDs2[i]])
1092         return false;
1093
1094     // map param on MA to parameters of nodes on a pair of theSinuEdges
1095     typedef map< double, pair< NodePoint, NodePoint > > TMAPar2NPoints;
1096     TMAPar2NPoints pointsOnE;
1097     vector<double> maParams;
1098
1099     // compute params of nodes on EDGEs by projecting division points from MA
1100     //const double tol = 1e-5 * theMAParams.back();
1101     size_t iEdgePair = 0;
1102     while ( iEdgePair < edgeIDs1.size() )
1103     {
1104       if ( isComputed[ edgeIDs1[ iEdgePair ]] ||
1105            isComputed[ edgeIDs2[ iEdgePair ]])
1106       {
1107         // "projection" from one side to the other
1108
1109         size_t iEdgeComputed = edgeIDs1[iEdgePair], iSideComputed = 0;
1110         if ( !isComputed[ iEdgeComputed ])
1111           ++iSideComputed, iEdgeComputed = edgeIDs2[iEdgePair];
1112
1113         map< double, const SMDS_MeshNode* > nodeParams; // params of existing nodes
1114         if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iEdgeComputed ], /*skipMedium=*/true, nodeParams ))
1115           return false;
1116
1117         SMESH_MAT2d::BoundaryPoint& bndPnt = bp[ 1-iSideComputed ];
1118         SMESH_MAT2d::BranchPoint brp;
1119         NodePoint npN, npB;
1120         NodePoint& np0 = iSideComputed ? npB : npN;
1121         NodePoint& np1 = iSideComputed ? npN : npB;
1122
1123         double maParam1st, maParamLast, maParam;
1124         if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.begin()->first, brp ))
1125             return false;
1126         branch.getParameter( brp, maParam1st );
1127         if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.rbegin()->first, brp ))
1128             return false;
1129         branch.getParameter( brp, maParamLast );
1130
1131         map< double, const SMDS_MeshNode* >::iterator u2n = nodeParams.begin(), u2nEnd = --nodeParams.end();
1132         TMAPar2NPoints::iterator end = pointsOnE.end(), pos = end;
1133         TMAPar2NPoints::iterator & hint = (maParamLast > maParam1st) ? end : pos;
1134         for ( ++u2n; u2n != u2nEnd; ++u2n )
1135         {
1136           if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, u2n->first, brp ))
1137             return false;
1138           if ( !branch.getBoundaryPoints( brp, bp[0], bp[1] ))
1139             return false;
1140           if ( !branch.getParameter( brp, maParam ))
1141             return false;
1142
1143           npN = NodePoint( u2n->second, u2n->first, iEdgeComputed );
1144           npB = NodePoint( bndPnt );
1145           pos = pointsOnE.insert( hint, make_pair( maParam, make_pair( np0, np1 )));
1146         }
1147
1148         // move iEdgePair forward
1149         while ( iEdgePair < edgeIDs1.size() )
1150           if ( edgeIDs1[ iEdgePair ] == bp[0]._edgeIndex &&
1151                edgeIDs2[ iEdgePair ] == bp[1]._edgeIndex )
1152             break;
1153           else
1154             ++iEdgePair;
1155       }
1156       else
1157       {
1158         // projection from MA
1159         maParams.clear();
1160         if ( !getParamsForEdgePair( iEdgePair, divPoints, theMAParams, maParams ))
1161           return false;
1162
1163         for ( size_t i = 1; i < maParams.size()-1; ++i )
1164         {
1165           if ( !branch.getBoundaryPoints( maParams[i], bp[0], bp[1] ))
1166             return false;
1167
1168           pointsOnE.insert( pointsOnE.end(), make_pair( maParams[i], make_pair( NodePoint(bp[0]),
1169                                                                                 NodePoint(bp[1]))));
1170         }
1171       }
1172       ++iEdgePair;
1173     }
1174
1175     if ( !projectVertices( theHelper, theMA, divPoints, edgeIDs1, edgeIDs2, theSinuEdges,
1176                            curves, isComputed, pointsOnE, theSinuFace._nodesToMerge ))
1177       return false;
1178
1179     // create nodes
1180     TMAPar2NPoints::iterator u2np = pointsOnE.begin();
1181     for ( ; u2np != pointsOnE.end(); ++u2np )
1182     {
1183       NodePoint* np[2] = { & u2np->second.first, & u2np->second.second };
1184       for ( int iSide = 0; iSide < 2; ++iSide )
1185       {
1186         if ( np[ iSide ]->_node ) continue;
1187         size_t       iEdge = np[ iSide ]->_edgeInd;
1188         double           u = np[ iSide ]->_u;
1189         gp_Pnt           p = curves[ iEdge ]->Value( u );
1190         np[ iSide ]->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() );
1191         meshDS->SetNodeOnEdge( np[ iSide ]->_node, edgeIDs[ iEdge ], u );
1192       }
1193     }
1194
1195     // create mesh segments on EDGEs
1196     theHelper.SetElementsOnShape( false );
1197     TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() );
1198     for ( size_t i = 0; i < theSinuEdges.size(); ++i )
1199     {
1200       SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] );
1201       if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 )
1202         continue;
1203
1204       StdMeshers_FaceSide side( face, theSinuEdges[i], mesh,
1205                                 /*isFwd=*/true, /*skipMediumNodes=*/true );
1206       vector<const SMDS_MeshNode*> nodes = side.GetOrderedNodes();
1207       for ( size_t in = 1; in < nodes.size(); ++in )
1208       {
1209         const SMDS_MeshElement* seg = theHelper.AddEdge( nodes[in-1], nodes[in], 0, false );
1210         meshDS->SetMeshElementOnShape( seg, edgeIDs[ i ] );
1211       }
1212     }
1213
1214     // update sub-meshes on VERTEXes
1215     for ( size_t i = 0; i < theSinuEdges.size(); ++i )
1216     {
1217       mesh->GetSubMesh( theHelper.IthVertex( 0, theSinuEdges[i] ))
1218         ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1219       mesh->GetSubMesh( theHelper.IthVertex( 1, theSinuEdges[i] ))
1220         ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1221     }
1222
1223     return true;
1224   }
1225
1226   //================================================================================
1227   /*!
1228    * \brief Mesh short EDGEs
1229    */
1230   //================================================================================
1231
1232   bool computeShortEdges( SMESH_MesherHelper&        theHelper,
1233                           const vector<TopoDS_Edge>& theShortEdges,
1234                           SMESH_Algo*                the1dAlgo )
1235   {
1236     for ( size_t i = 0; i < theShortEdges.size(); ++i )
1237     {
1238       theHelper.GetGen()->Compute( *theHelper.GetMesh(), theShortEdges[i], true, true );
1239
1240       SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh(theShortEdges[i] );
1241       if ( sm->IsEmpty() )
1242       {
1243         try {
1244           if ( !the1dAlgo->Compute( *theHelper.GetMesh(), theShortEdges[i] ))
1245             return false;
1246         }
1247         catch (...) {
1248           return false;
1249         }
1250         sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1251         if ( sm->IsEmpty() )
1252           return false;
1253       }
1254     }
1255     return true;
1256   }
1257
1258   inline double area( const UVPtStruct& p1, const UVPtStruct& p2, const UVPtStruct& p3 )
1259   {
1260     gp_XY v1 = p2.UV() - p1.UV();
1261     gp_XY v2 = p3.UV() - p1.UV();
1262     return v2 ^ v1;
1263   }
1264
1265   bool ellipticSmooth( FaceQuadStruct::Ptr quad, int nbLoops )
1266   {
1267     //nbLoops = 10;
1268     if ( quad->uv_grid.empty() )
1269       return true;
1270
1271     int nbhoriz  = quad->iSize;
1272     int nbvertic = quad->jSize;
1273
1274     const double dksi = 0.5, deta = 0.5;
1275     const double  dksi2 = dksi*dksi, deta2 = deta*deta;
1276     double err = 0., g11, g22, g12;
1277     int nbErr = 0;
1278
1279     FaceQuadStruct& q = *quad;
1280     UVPtStruct pNew;
1281
1282     double refArea = area( q.UVPt(0,0), q.UVPt(1,0), q.UVPt(1,1) );
1283
1284     for ( int iLoop = 0; iLoop < nbLoops; ++iLoop )
1285     {
1286       err = 0;
1287       for ( int i = 1; i < nbhoriz - 1; i++ )
1288         for ( int j = 1; j < nbvertic - 1; j++ )
1289         {
1290           g11 = ( (q.U(i,j+1) - q.U(i,j-1))*(q.U(i,j+1) - q.U(i,j-1))/dksi2 +
1291                   (q.V(i,j+1) - q.V(i,j-1))*(q.V(i,j+1) - q.V(i,j-1))/deta2 )/4;
1292
1293           g22 = ( (q.U(i+1,j) - q.U(i-1,j))*(q.U(i+1,j) - q.U(i-1,j))/dksi2 +
1294                   (q.V(i+1,j) - q.V(i-1,j))*(q.V(i+1,j) - q.V(i-1,j))/deta2 )/4;
1295
1296           g12 = ( (q.U(i+1,j) - q.U(i-1,j))*(q.U(i,j+1) - q.U(i,j-1))/dksi2 +
1297                   (q.V(i+1,j) - q.V(i-1,j))*(q.V(i,j+1) - q.V(i,j-1))/deta2 )/(4*dksi*deta);
1298
1299           pNew.u = dksi2/(2*(g11+g22)) * (g11*(q.U(i+1,j) + q.U(i-1,j))/dksi2 +
1300                                           g22*(q.U(i,j+1) + q.U(i,j-1))/dksi2
1301                                           - 0.5*g12*q.U(i+1,j+1) + 0.5*g12*q.U(i-1,j+1) +
1302                                           - 0.5*g12*q.U(i-1,j-1) + 0.5*g12*q.U(i+1,j-1));
1303
1304           pNew.v = deta2/(2*(g11+g22)) * (g11*(q.V(i+1,j) + q.V(i-1,j))/deta2 +
1305                                           g22*(q.V(i,j+1) + q.V(i,j-1))/deta2
1306                                           - 0.5*g12*q.V(i+1,j+1) + 0.5*g12*q.V(i-1,j+1) +
1307                                           - 0.5*g12*q.V(i-1,j-1) + 0.5*g12*q.V(i+1,j-1));
1308
1309           // if (( refArea * area( q.UVPt(i-1,j-1), q.UVPt(i,j-1), pNew ) > 0 ) &&
1310           //     ( refArea * area( q.UVPt(i+1,j-1), q.UVPt(i+1,j), pNew ) > 0 ) &&
1311           //     ( refArea * area( q.UVPt(i+1,j+1), q.UVPt(i,j+1), pNew ) > 0 ) &&
1312           //     ( refArea * area( q.UVPt(i-1,j), q.UVPt(i-1,j-1), pNew ) > 0 ))
1313           {
1314             err += sqrt(( q.U(i,j) - pNew.u ) * ( q.U(i,j) - pNew.u ) +
1315                         ( q.V(i,j) - pNew.v ) * ( q.V(i,j) - pNew.v ));
1316             q.U(i,j) = pNew.u;
1317             q.V(i,j) = pNew.v;
1318           }
1319           // else if ( ++nbErr < 10 )
1320           // {
1321           //   cout << i << ", " << j << endl;
1322           //   cout << "x = ["
1323           //        << "[ " << q.U(i-1,j-1) << ", " <<q.U(i,j-1) << ", " << q.U(i+1,j-1) << " ],"
1324           //        << "[ " << q.U(i-1,j-0) << ", " <<q.U(i,j-0) << ", " << q.U(i+1,j-0) << " ],"
1325           //        << "[ " << q.U(i-1,j+1) << ", " <<q.U(i,j+1) << ", " << q.U(i+1,j+1) << " ]]" << endl;
1326           //   cout << "y = ["
1327           //        << "[ " << q.V(i-1,j-1) << ", " <<q.V(i,j-1) << ", " << q.V(i+1,j-1) << " ],"
1328           //        << "[ " << q.V(i-1,j-0) << ", " <<q.V(i,j-0) << ", " << q.V(i+1,j-0) << " ],"
1329           //        << "[ " << q.V(i-1,j+1) << ", " <<q.V(i,j+1) << ", " << q.V(i+1,j+1) << " ]]" << endl<<endl;
1330           // }
1331         }
1332
1333       if ( err / ( nbhoriz - 2 ) / ( nbvertic - 2 ) < 1e-6 )
1334         break;
1335     }
1336     //cout << " ERR " << err / ( nbhoriz - 2 ) / ( nbvertic - 2 ) << endl;
1337
1338     return true;
1339   }
1340
1341   //================================================================================
1342   /*!
1343    * \brief Remove temporary node
1344    */
1345   //================================================================================
1346
1347   void mergeNodes( SMESH_MesherHelper& theHelper,
1348                    SinuousFace&        theSinuFace )
1349   {
1350     SMESH_MeshEditor editor( theHelper.GetMesh() );
1351     SMESH_MeshEditor::TListOfListOfNodes nodesGroups;
1352
1353     TMergeMap::iterator n2nn = theSinuFace._nodesToMerge.begin();
1354     for ( ; n2nn != theSinuFace._nodesToMerge.end(); ++n2nn )
1355     {
1356       nodesGroups.push_back( list< const SMDS_MeshNode* >() );
1357       list< const SMDS_MeshNode* > & group = nodesGroups.back();
1358
1359       group.push_back( n2nn->first );
1360       group.splice( group.end(), n2nn->second );
1361     }
1362     editor.MergeNodes( nodesGroups );
1363   }
1364
1365 } // namespace
1366
1367 //================================================================================
1368 /*!
1369  * \brief Create quadrangle elements
1370  *  \param [in] theHelper - the helper
1371  *  \param [in] theFace - the face to mesh
1372  *  \param [in] theSinuEdges - the sinuous EDGEs
1373  *  \param [in] theShortEdges - the short EDGEs
1374  *  \return bool - is OK or not
1375  */
1376 //================================================================================
1377
1378 bool StdMeshers_QuadFromMedialAxis_1D2D::computeQuads( SMESH_MesherHelper&       theHelper,
1379                                                        const TopoDS_Face&        theFace,
1380                                                        const vector<TopoDS_Edge> theSinuEdges[2],
1381                                                        const vector<TopoDS_Edge> theShortEdges[2])
1382 {
1383   SMESH_Mesh* mesh = theHelper.GetMesh();
1384   SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( *mesh, theFace );
1385   if ( !proxyMesh )
1386     return false;
1387
1388   StdMeshers_Quadrangle_2D::myProxyMesh  = proxyMesh;
1389   StdMeshers_Quadrangle_2D::myHelper     = &theHelper;
1390   StdMeshers_Quadrangle_2D::myNeedSmooth = false;
1391   StdMeshers_Quadrangle_2D::myCheckOri   = false;
1392   StdMeshers_Quadrangle_2D::myQuadList.clear();
1393
1394   // fill FaceQuadStruct
1395
1396   list< TopoDS_Edge > side[4];
1397   side[0].insert( side[0].end(), theShortEdges[0].begin(), theShortEdges[0].end() );
1398   side[1].insert( side[1].end(), theSinuEdges[1].begin(),  theSinuEdges[1].end() );
1399   side[2].insert( side[2].end(), theShortEdges[1].begin(), theShortEdges[1].end() );
1400   side[3].insert( side[3].end(), theSinuEdges[0].begin(),  theSinuEdges[0].end() );
1401
1402   FaceQuadStruct::Ptr quad( new FaceQuadStruct );
1403   quad->side.resize( 4 );
1404   quad->face = theFace;
1405   for ( int i = 0; i < 4; ++i )
1406   {
1407     quad->side[i] = StdMeshers_FaceSide::New( theFace, side[i], mesh, i < QUAD_TOP_SIDE,
1408                                               /*skipMediumNodes=*/true, proxyMesh );
1409   }
1410   int nbNodesShort0 = quad->side[0].NbPoints();
1411   int nbNodesShort1 = quad->side[2].NbPoints();
1412
1413   // compute UV of internal points
1414   myQuadList.push_back( quad );
1415   if ( !StdMeshers_Quadrangle_2D::setNormalizedGrid( quad ))
1416     return false;
1417
1418   // elliptic smooth of internal points to get boundary cell normal to the boundary
1419   ellipticSmooth( quad, 1 );
1420
1421   // create quadrangles
1422   bool ok;
1423   if ( nbNodesShort0 == nbNodesShort1 )
1424     ok = StdMeshers_Quadrangle_2D::computeQuadDominant( *mesh, theFace, quad );
1425   else
1426     ok = StdMeshers_Quadrangle_2D::computeTriangles( *mesh, theFace, quad );
1427
1428   StdMeshers_Quadrangle_2D::myProxyMesh.reset();
1429   StdMeshers_Quadrangle_2D::myHelper = 0;
1430
1431   return ok;
1432 }
1433
1434 //================================================================================
1435 /*!
1436  * \brief Generate quadrangle mesh
1437  */
1438 //================================================================================
1439
1440 bool StdMeshers_QuadFromMedialAxis_1D2D::Compute(SMESH_Mesh&         theMesh,
1441                                                  const TopoDS_Shape& theShape)
1442 {
1443   SMESH_MesherHelper helper( theMesh );
1444   helper.SetSubShape( theShape );
1445
1446   TopoDS_Face F = TopoDS::Face( theShape );
1447   if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD );
1448
1449   SinuousFace sinuFace( F );
1450
1451   _progress = 0.01;
1452
1453   if ( getSinuousEdges( helper, sinuFace ))
1454   {
1455     _progress = 0.2;
1456
1457     // if ( sinuFace._sinuEdges.size() > 2 )
1458     //   return error(COMPERR_BAD_SHAPE, "Not yet supported case" );
1459
1460     double minSegLen = getMinSegLen( helper, sinuFace._sinuEdges );
1461     SMESH_MAT2d::MedialAxis ma( F, sinuFace._sinuEdges, minSegLen, /*ignoreCorners=*/true );
1462
1463     if ( !_regular1D )
1464       _regular1D = new Algo1D( _studyId, _gen );
1465     _regular1D->SetSegmentLength( minSegLen );
1466
1467     vector<double> maParams;
1468     if ( ! divideMA( helper, ma, sinuFace, _regular1D, minSegLen, maParams ))
1469       return error(COMPERR_BAD_SHAPE);
1470
1471     _progress = 0.4;
1472
1473     if ( !computeShortEdges( helper, sinuFace._shortSide[0], _regular1D ) ||
1474          !computeShortEdges( helper, sinuFace._shortSide[1], _regular1D ))
1475       return error("Failed to mesh short edges");
1476
1477     _progress = 0.6;
1478
1479     if ( !computeSinuEdges( helper, minSegLen, ma, maParams, sinuFace ))
1480       return error("Failed to mesh sinuous edges");
1481
1482     _progress = 0.8;
1483
1484     bool ok = computeQuads( helper, F, sinuFace._sinuSide, sinuFace._shortSide );
1485
1486     if ( ok )
1487       mergeNodes( helper, sinuFace );
1488
1489     _progress = 1.;
1490
1491     return ok;
1492   }
1493
1494   return error(COMPERR_BAD_SHAPE, "Not implemented so far");
1495 }
1496
1497 //================================================================================
1498 /*!
1499  * \brief Predict nb of elements
1500  */
1501 //================================================================================
1502
1503 bool StdMeshers_QuadFromMedialAxis_1D2D::Evaluate(SMESH_Mesh &         theMesh,
1504                                                   const TopoDS_Shape & theShape,
1505                                                   MapShapeNbElems&     theResMap)
1506 {
1507   return StdMeshers_Quadrangle_2D::Evaluate(theMesh,theShape,theResMap);
1508 }
1509