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