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