src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx

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