src/StdMeshers/StdMeshers_FaceSide.cxx

   1 // Copyright (C) 2007-2014  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
  23 // File      : StdMeshers_FaceSide.hxx
  24 // Created   : Wed Jan 31 18:41:25 2007
  25 // Author    : Edward AGAPOV (eap)
  26 // Module    : SMESH
  27 //
  28 #include "StdMeshers_FaceSide.hxx"
  29
  30 #include "SMDS_EdgePosition.hxx"
  31 #include "SMDS_MeshNode.hxx"
  32 #include "SMESHDS_Mesh.hxx"
  33 #include "SMESHDS_SubMesh.hxx"
  34 #include "SMESH_Algo.hxx"
  35 #include "SMESH_Mesh.hxx"
  36 #include "SMESH_MesherHelper.hxx"
  37 #include "SMESH_ComputeError.hxx"
  38 #include "SMESH_Block.hxx"
  39
  40 #include <Adaptor2d_Curve2d.hxx>
  41 #include <BRepAdaptor_CompCurve.hxx>
  42 #include <BRep_Builder.hxx>
  43 #include <BRep_Tool.hxx>
  44 #include <GCPnts_AbscissaPoint.hxx>
  45 #include <Geom2dAdaptor_Curve.hxx>
  46 #include <TopExp.hxx>
  47 #include <TopExp_Explorer.hxx>
  48 #include <TopoDS.hxx>
  49 #include <TopoDS_Face.hxx>
  50 #include <TopoDS_Vertex.hxx>
  51 #include <TopoDS_Wire.hxx>
  52
  53 #include <map>
  54 #include <limits>
  55
  56 #include "utilities.h"
  57
  58 //================================================================================
  59 /*!
  60  * \brief Constructor of a side of one edge
  61   * \param theFace - the face
  62   * \param theEdge - the edge
  63  */
  64 //================================================================================
  65
  66 StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face&   theFace,
  67                                          const TopoDS_Edge&   theEdge,
  68                                          SMESH_Mesh*          theMesh,
  69                                          const bool           theIsForward,
  70                                          const bool           theIgnoreMediumNodes,
  71                                          SMESH_ProxyMesh::Ptr theProxyMesh)
  72 {
  73   list<TopoDS_Edge> edges(1,theEdge);
  74   *this = StdMeshers_FaceSide( theFace, edges, theMesh, theIsForward,
  75                                theIgnoreMediumNodes, theProxyMesh );
  76 }
  77
  78 //================================================================================
  79 /*!
  80  * \brief Constructor of a side of several edges
  81   * \param theFace - the face
  82   * \param theEdge - the edge
  83  */
  84 //================================================================================
  85
  86 StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face&   theFace,
  87                                          list<TopoDS_Edge>&   theEdges,
  88                                          SMESH_Mesh*          theMesh,
  89                                          const bool           theIsForward,
  90                                          const bool           theIgnoreMediumNodes,
  91                                          SMESH_ProxyMesh::Ptr theProxyMesh)
  92 {
  93   int nbEdges = theEdges.size();
  94   myEdge.resize      ( nbEdges );
  95   myEdgeID.resize    ( nbEdges );
  96   myC2d.resize       ( nbEdges );
  97   myC3dAdaptor.resize( nbEdges );
  98   myFirst.resize     ( nbEdges );
  99   myLast.resize      ( nbEdges );
 100   myNormPar.resize   ( nbEdges );
 101   myEdgeLength.resize( nbEdges );
 102   myIsUniform.resize ( nbEdges, true );
 103   myLength             = 0;
 104   myNbPonits           = myNbSegments = 0;
 105   myProxyMesh          = theProxyMesh;
 106   myMissingVertexNodes = false;
 107   myIgnoreMediumNodes  = theIgnoreMediumNodes;
 108   myDefaultPnt2d       = gp_Pnt2d( 1e+100, 1e+100 );
 109   if ( !myProxyMesh ) myProxyMesh.reset( new SMESH_ProxyMesh( *theMesh ));
 110   if ( nbEdges == 0 ) return;
 111
 112   SMESHDS_Mesh* meshDS = myProxyMesh->GetMeshDS();
 113
 114   int nbDegen = 0;
 115   list<TopoDS_Edge>::iterator edge = theEdges.begin();
 116   TopoDS_Iterator vExp;
 117   for ( int index = 0; edge != theEdges.end(); ++index, ++edge )
 118   {
 119     int i = theIsForward ? index : nbEdges-index-1;
 120     myEdgeLength[i] = SMESH_Algo::EdgeLength( *edge );
 121     if ( myEdgeLength[i] < DBL_MIN ) nbDegen++;
 122     myLength += myEdgeLength[i];
 123     myEdge[i] = *edge;
 124     myEdgeID[i] = meshDS->ShapeToIndex( *edge );
 125     if ( !theIsForward ) myEdge[i].Reverse();
 126
 127     if ( theFace.IsNull() )
 128       BRep_Tool::Range( *edge, myFirst[i], myLast[i] );
 129     else
 130       myC2d[i] = BRep_Tool::CurveOnSurface( *edge, theFace, myFirst[i], myLast[i] );
 131     if ( myEdge[i].Orientation() == TopAbs_REVERSED )
 132       std::swap( myFirst[i], myLast[i] );
 133
 134     if ( const SMESHDS_SubMesh* sm = myProxyMesh->GetSubMesh( *edge )) {
 135       int nbN = sm->NbNodes();
 136       if ( theIgnoreMediumNodes ) {
 137         SMDS_ElemIteratorPtr elemIt = sm->GetElements();
 138         if ( elemIt->more() && elemIt->next()->IsQuadratic() )
 139           nbN -= sm->NbElements();
 140       }
 141       myNbPonits += nbN;
 142       myNbSegments += sm->NbElements();
 143     }
 144
 145     // TopExp::FirstVertex() and TopExp::LastVertex() return NULL from INTERNAL edge
 146     vExp.Initialize( *edge );
 147     if ( vExp.Value().Orientation() == TopAbs_REVERSED ) vExp.Next();
 148     if ( SMESH_Algo::VertexNode( TopoDS::Vertex( vExp.Value()), meshDS ))
 149       myNbPonits += 1; // for the first end
 150     else
 151       myMissingVertexNodes = true;
 152
 153     // check if the edge has a non-uniform parametrization (issue 0020705)
 154     if ( !myC2d[i].IsNull() && myEdgeLength[i] > DBL_MIN)
 155     {
 156       Geom2dAdaptor_Curve A2dC( myC2d[i],
 157                                 std::min( myFirst[i], myLast[i] ),
 158                                 std::max( myFirst[i], myLast[i] ));
 159       double p2 = myFirst[i]+(myLast[i]-myFirst[i])/2., p4 = myFirst[i]+(myLast[i]-myFirst[i])/4.;
 160       double d2 = GCPnts_AbscissaPoint::Length( A2dC, myFirst[i], p2 );
 161       double d4 = GCPnts_AbscissaPoint::Length( A2dC, myFirst[i], p4 );
 162       //cout<<"len = "<<len<<"  d2 = "<<d2<<"  fabs(2*d2/len-1.0) = "<<fabs(2*d2/len-1.0)<<endl;
 163       myIsUniform[i] = !( fabs(2*d2/myEdgeLength[i]-1.0) > 0.01 || fabs(2*d4/d2-1.0) > 0.01 );
 164       //if ( !myIsUniform[i] ) to implement Value3d(u)
 165       {
 166         double fp,lp;
 167         Handle(Geom_Curve) C3d = BRep_Tool::Curve(myEdge[i],fp,lp);
 168         myC3dAdaptor[i].Load( C3d, fp,lp );
 169       }
 170     }
 171     // reverse a proxy submesh
 172     if ( !theIsForward )
 173       reverseProxySubmesh( myEdge[i] );
 174
 175   } // loop on edges
 176
 177   vExp.Initialize( theEdges.back() );
 178   if ( vExp.Value().Orientation() != TopAbs_REVERSED ) vExp.Next();
 179   if ( vExp.More() )
 180   {
 181     if ( SMESH_Algo::VertexNode( TopoDS::Vertex( vExp.Value()), meshDS ))
 182       myNbPonits++; // for the last end
 183     else
 184       myMissingVertexNodes = true;
 185   }
 186   if ( nbEdges > 1 && myLength > DBL_MIN ) {
 187     const double degenNormLen = 1.e-5;
 188     double totLength = myLength;
 189     if ( nbDegen )
 190       totLength += myLength * degenNormLen * nbDegen;
 191     double prevNormPar = 0;
 192     for ( int i = 0; i < nbEdges; ++i ) {
 193       if ( myEdgeLength[ i ] < DBL_MIN )
 194         myEdgeLength[ i ] = myLength * degenNormLen;
 195       myNormPar[ i ] = prevNormPar + myEdgeLength[i]/totLength;
 196       prevNormPar = myNormPar[ i ];
 197     }
 198   }
 199   myNormPar[nbEdges-1] = 1.;
 200   //dump();
 201 }
 202
 203 //================================================================================
 204 /*!
 205  * \brief Constructor of a side for vertex using data from other FaceSide
 206  *  \param theVertex - the vertex
 207  *  \param theSide - the side
 208  */
 209 //================================================================================
 210
 211 StdMeshers_FaceSide::StdMeshers_FaceSide(const StdMeshers_FaceSide*  theSide,
 212                                          const SMDS_MeshNode*        theNode,
 213                                          const gp_Pnt2d*             thePnt2d1,
 214                                          const gp_Pnt2d*             thePnt2d2,
 215                                          const Handle(Geom2d_Curve)& theC2d,
 216                                          const double                theUFirst,
 217                                          const double                theULast)
 218 {
 219   myC2d.push_back      ( theC2d );
 220   myFirst.push_back    ( theUFirst );
 221   myLast.push_back     ( theULast );
 222   myNormPar.push_back  ( 1. );
 223   myIsUniform.push_back( true );
 224   myEdgeID.push_back   ( 0 );
 225   myLength       = 0;
 226   myProxyMesh    = theSide->myProxyMesh;
 227   myDefaultPnt2d = *thePnt2d1;
 228   myPoints       = theSide->GetUVPtStruct();
 229   myNbPonits     = myPoints.size();
 230   myNbSegments   = theSide->myNbSegments;
 231   if ( thePnt2d2 )
 232     for ( size_t i = 0; i < myPoints.size(); ++i )
 233     {
 234       double r = i / ( myPoints.size() - 1. );
 235       myPoints[i].u = (1-r) * thePnt2d1->X() + r * thePnt2d2->X();
 236       myPoints[i].v = (1-r) * thePnt2d1->Y() + r * thePnt2d2->Y();
 237       myPoints[i].node = theNode;
 238     }
 239   else
 240     for ( size_t i = 0; i < myPoints.size(); ++i )
 241     {
 242       myPoints[i].u = thePnt2d1->X();
 243       myPoints[i].v = thePnt2d1->Y();
 244       myPoints[i].node = theNode;
 245     }
 246 }
 247
 248 //================================================================================
 249 /*
 250  * Create a side from an UVPtStructVec
 251  */
 252 //================================================================================
 253
 254 StdMeshers_FaceSide::StdMeshers_FaceSide(UVPtStructVec&     theSideNodes,
 255                                          const TopoDS_Face& theFace)
 256 {
 257   myEdge.resize( 1 );
 258   myEdgeID.resize( 1, -1 );
 259   myC2d.resize( 1 );
 260   myC3dAdaptor.resize( 1 );
 261   myFirst.resize( 1, 0. );
 262   myLast.resize( 1, 1. );
 263   myNormPar.resize( 1, 1. );
 264   myIsUniform.resize( 1, 1 );
 265   myMissingVertexNodes = myIgnoreMediumNodes = false;
 266   myDefaultPnt2d.SetCoord( 1e100, 1e100 );
 267
 268   myPoints     = theSideNodes;
 269   myNbPonits   = myPoints.size();
 270   myNbSegments = myNbPonits + 1;
 271
 272   myLength = 0;
 273   if ( !myPoints.empty() )
 274   {
 275     myPoints[0].normParam = 0;
 276     if ( myPoints[0].node &&
 277          myPoints.back().node &&
 278          myPoints[ myNbPonits/2 ].node )
 279     {
 280       gp_Pnt pPrev = SMESH_TNodeXYZ( myPoints[0].node );
 281       for ( size_t i = 1; i < myPoints.size(); ++i )
 282       {
 283         gp_Pnt p = SMESH_TNodeXYZ( myPoints[i].node );
 284         myLength += p.Distance( pPrev );
 285         myPoints[i].normParam = myLength;
 286         pPrev = p;
 287       }
 288     }
 289     else if ( !theFace.IsNull() )
 290     {
 291       TopLoc_Location loc;
 292       Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace, loc );
 293       gp_Pnt pPrev = surf->Value( myPoints[0].u, myPoints[0].v );
 294       for ( size_t i = 1; i < myPoints.size(); ++i )
 295       {
 296         gp_Pnt p = surf->Value( myPoints[i].u, myPoints[i].v );
 297         myLength += p.Distance( pPrev );
 298         myPoints[i].normParam = myLength;
 299         pPrev = p;
 300       }
 301     }
 302     else
 303     {
 304       gp_Pnt2d pPrev = myPoints[0].UV();
 305       for ( size_t i = 1; i < myPoints.size(); ++i )
 306       {
 307         gp_Pnt2d p = myPoints[i].UV();
 308         myLength += p.Distance( pPrev );
 309         myPoints[i].normParam = myLength;
 310         pPrev = p;
 311       }
 312     }
 313     if ( myLength > std::numeric_limits<double>::min() )
 314       for ( size_t i = 1; i < myPoints.size(); ++i )
 315         myPoints[i].normParam /= myLength;
 316   }
 317   myEdgeLength.resize( 1, myLength );
 318 }
 319
 320 //================================================================================
 321 /*
 322  * Return info on nodes on the side
 323  */
 324 //================================================================================
 325
 326 const vector<UVPtStruct>& StdMeshers_FaceSide::GetUVPtStruct(bool   isXConst,
 327                                                              double constValue) const
 328 {
 329   if ( myPoints.empty() ) {
 330
 331     if ( NbEdges() == 0 ) return myPoints;
 332
 333     SMESHDS_Mesh* meshDS = myProxyMesh->GetMeshDS();
 334     SMESH_MesherHelper helper(*myProxyMesh->GetMesh());
 335     bool paramOK;
 336     double eps = 1e-100;
 337
 338     // sort nodes of all edges putting them into a map
 339
 340     map< double, const SMDS_MeshNode*> u2node;
 341     vector< const SMESH_ProxyMesh::SubMesh* > proxySubMesh( myEdge.size());
 342     int nbProxyNodes = 0;
 343     for ( size_t iE = 0; iE < myEdge.size(); ++iE )
 344     {
 345       proxySubMesh[iE] = myProxyMesh->GetProxySubMesh( myEdge[iE] );
 346       if ( proxySubMesh[iE] )
 347       {
 348         if ( proxySubMesh[iE]->GetUVPtStructVec().empty() ) {
 349           proxySubMesh[iE] = 0;
 350         }
 351         else {
 352           nbProxyNodes += proxySubMesh[iE]->GetUVPtStructVec().size() - 1;
 353           if ( iE+1 == myEdge.size() )
 354             ++nbProxyNodes;
 355           continue;
 356         }
 357       }
 358       // Put 1st vertex node of a current edge
 359       TopoDS_Vertex VV[2]; // TopExp::FirstVertex() returns NULL for INTERNAL edge
 360       VV[0] = SMESH_MesherHelper::IthVertex( 0, myEdge[iE]);
 361       VV[1] = SMESH_MesherHelper::IthVertex( 1, myEdge[iE]);
 362       const SMDS_MeshNode* node = SMESH_Algo::VertexNode( VV[0], meshDS );
 363       double prevNormPar = ( iE == 0 ? 0 : myNormPar[ iE-1 ]); // normalized param
 364       if ( node ) { // nodes on internal vertices may be missing
 365         u2node.insert( u2node.end(), make_pair( prevNormPar, node ));
 366       }
 367       else if ( iE == 0 ) {
 368         MESSAGE(" NO NODE on VERTEX" );
 369         return myPoints;
 370       }
 371
 372       // Put internal nodes
 373       if ( const SMESHDS_SubMesh* sm = myProxyMesh->GetSubMesh( myEdge[iE] ))
 374       {
 375         vector< pair< double, const SMDS_MeshNode*> > u2nodeVec;
 376         u2nodeVec.reserve( sm->NbNodes() );
 377         SMDS_NodeIteratorPtr nItr = sm->GetNodes();
 378         double paramSize = myLast[iE] - myFirst[iE];
 379         double r         = myNormPar[iE] - prevNormPar;
 380         helper.SetSubShape( myEdge[iE] );
 381         helper.ToFixNodeParameters( true );
 382         if ( !myIsUniform[iE] )
 383           while ( nItr->more() )
 384           {
 385             const SMDS_MeshNode* node = nItr->next();
 386             if ( myIgnoreMediumNodes && SMESH_MeshEditor::IsMedium( node, SMDSAbs_Edge ))
 387               continue;
 388             double u = helper.GetNodeU( myEdge[iE], node, 0, &paramOK );
 389             double aLenU = GCPnts_AbscissaPoint::Length
 390               ( const_cast<GeomAdaptor_Curve&>( myC3dAdaptor[iE]), myFirst[iE], u );
 391             if ( myEdgeLength[iE] < aLenU ) // nonregression test "3D_mesh_NETGEN/G6"
 392             {
 393               u2nodeVec.clear();
 394               break;
 395             }
 396             double normPar = prevNormPar + r*aLenU/myEdgeLength[iE];
 397             u2nodeVec.push_back( make_pair( normPar, node ));
 398           }
 399         nItr = sm->GetNodes();
 400         if ( u2nodeVec.empty() )
 401           while ( nItr->more() )
 402           {
 403             const SMDS_MeshNode* node = nItr->next();
 404             if ( myIgnoreMediumNodes && SMESH_MeshEditor::IsMedium( node, SMDSAbs_Edge ))
 405               continue;
 406             double u = helper.GetNodeU( myEdge[iE], node, 0, &paramOK );
 407
 408             // paramSize is signed so orientation is taken into account
 409             double normPar = prevNormPar + r * ( u - myFirst[iE] ) / paramSize;
 410             u2nodeVec.push_back( make_pair( normPar, node ));
 411           }
 412         for ( size_t j = 0; j < u2nodeVec.size(); ++j )
 413           u2node.insert( u2node.end(), u2nodeVec[j] );
 414       }
 415
 416       // Put 2nd vertex node for a last edge
 417       if ( iE+1 == myEdge.size() ) {
 418         node = SMESH_Algo::VertexNode( VV[1], meshDS );
 419         if ( !node ) {
 420           MESSAGE(" NO NODE on VERTEX" );
 421           return myPoints;
 422         }
 423         u2node.insert( u2node.end(), make_pair( 1., node ));
 424       }
 425     } // loop on myEdge's
 426
 427     if ( u2node.size() + nbProxyNodes != myNbPonits &&
 428          u2node.size() + nbProxyNodes != NbPoints( /*update=*/true ))
 429     {
 430       MESSAGE("Wrong node parameters on edges, u2node.size():"
 431               <<u2node.size()<<" !=  myNbPonits:"<<myNbPonits);
 432       return myPoints;
 433     }
 434
 435     // fill array of UVPtStruct
 436
 437     UVPtStructVec& points = const_cast< UVPtStructVec& >( myPoints );
 438     points.resize( myNbPonits );
 439
 440     int iPt = 0;
 441     double prevNormPar = 0, paramSize = myNormPar[ 0 ];
 442     map< double, const SMDS_MeshNode*>::iterator u_node = u2node.begin();
 443     for ( size_t iE = 0; iE < myEdge.size(); ++iE )
 444     {
 445       if ( proxySubMesh[ iE ] ) // copy data from a proxy sub-mesh
 446       {
 447         const UVPtStructVec& edgeUVPtStruct = proxySubMesh[iE]->GetUVPtStructVec();
 448         std::copy( edgeUVPtStruct.begin(), edgeUVPtStruct.end(), & points[iPt] );
 449         // check orientation
 450         double du1 = edgeUVPtStruct.back().param - edgeUVPtStruct[0].param;
 451         double du2 = myLast[iE] - myFirst[iE];
 452         if ( du1 * du2 < 0 )
 453         {
 454           std::reverse( & points[iPt], & points[iPt + edgeUVPtStruct.size()]);
 455           for ( size_t i = 0; i < edgeUVPtStruct.size(); ++i )
 456             points[iPt+i].normParam = 1. - points[iPt+i].normParam;
 457         }
 458         // update normalized params
 459         if ( myEdge.size() > 1 ) {
 460           for ( size_t i = 0; i < edgeUVPtStruct.size(); ++i, ++iPt )
 461           {
 462             UVPtStruct & uvPt = points[iPt];
 463             uvPt.normParam    = prevNormPar + uvPt.normParam * paramSize;
 464             uvPt.x = uvPt.y   = uvPt.normParam;
 465           }
 466           --iPt; // to point to the 1st VERTEX of the next EDGE
 467         }
 468       }
 469       else
 470       {
 471         for ( ; u_node != u2node.end(); ++u_node, ++iPt )
 472         {
 473           if ( myNormPar[ iE ]-eps < u_node->first )
 474             break; // u_node is at VERTEX of the next EDGE
 475
 476           UVPtStruct & uvPt = points[iPt];
 477           uvPt.node       = u_node->second;
 478           // -- normParam, x, y --------------------------------
 479           uvPt.normParam  = u_node->first;
 480           uvPt.x = uvPt.y = uvPt.normParam;
 481           // -- U ----------------------------------------------
 482           const SMDS_EdgePosition* epos =
 483             dynamic_cast<const SMDS_EdgePosition*>(uvPt.node->GetPosition());
 484           if ( epos ) {
 485             uvPt.param = epos->GetUParameter();
 486           }
 487           else {
 488             double r = ( uvPt.normParam - prevNormPar )/ paramSize;
 489             uvPt.param = ( r > 0.5 ? myLast[iE] : myFirst[iE] );
 490           }
 491           // -- UV ---------------------------------------------
 492           if ( !myC2d[ iE ].IsNull() ) {
 493             gp_Pnt2d p = myC2d[ iE ]->Value( uvPt.param );
 494             uvPt.u = p.X();
 495             uvPt.v = p.Y();
 496           }
 497           else {
 498             uvPt.u = uvPt.v = 1e+100;
 499           }
 500         }
 501       }
 502       // prepare for the next EDGE
 503       if ( iE+1 < myEdge.size() )
 504       {
 505         prevNormPar = myNormPar[ iE ];
 506         paramSize   = myNormPar[ iE+1 ] - prevNormPar;
 507       }
 508     } // loop on myEdge's
 509
 510     // set <constValue>
 511     if ( isXConst )
 512       for ( iPt = 0; iPt < points.size(); ++iPt ) points[ iPt ].x = constValue;
 513     else
 514       for ( iPt = 0; iPt < points.size(); ++iPt ) points[ iPt ].y = constValue;
 515
 516   } // if ( myPoints.empty())
 517
 518   return myPoints;
 519 }
 520
 521 //================================================================================
 522 /*!
 523  * \brief Falsificate info on nodes
 524  * \param nbSeg - nb of segments on the side
 525  * \retval UVPtStruct* - array of data structures
 526  */
 527 //================================================================================
 528
 529 const vector<UVPtStruct>& StdMeshers_FaceSide::SimulateUVPtStruct(int    nbSeg,
 530                                                                   bool   isXConst,
 531                                                                   double constValue) const
 532 {
 533   if ( myFalsePoints.empty() ) {
 534
 535     if ( NbEdges() == 0 ) return myFalsePoints;
 536
 537     vector<uvPtStruct>* points = const_cast<vector<uvPtStruct>*>( &myFalsePoints );
 538     points->resize( nbSeg+1 );
 539
 540     int EdgeIndex = 0;
 541     double prevNormPar = 0, paramSize = myNormPar[ EdgeIndex ];
 542     for (int i = 0 ; i < myFalsePoints.size(); ++i ) {
 543       double normPar = double(i) / double(nbSeg);
 544       UVPtStruct & uvPt = (*points)[i];
 545       uvPt.node = 0;
 546       uvPt.x = uvPt.y = uvPt.param = uvPt.normParam = normPar;
 547       if ( isXConst ) uvPt.x = constValue;
 548       else            uvPt.y = constValue;
 549       if ( myNormPar[ EdgeIndex ] < normPar ) {
 550         prevNormPar = myNormPar[ EdgeIndex ];
 551         ++EdgeIndex;
 552         paramSize = myNormPar[ EdgeIndex ] - prevNormPar;
 553       }
 554       double r = ( normPar - prevNormPar )/ paramSize;
 555       uvPt.param = myFirst[EdgeIndex] * ( 1 - r ) + myLast[EdgeIndex] * r;
 556       if ( !myC2d[ EdgeIndex ].IsNull() ) {
 557         gp_Pnt2d p = myC2d[ EdgeIndex ]->Value( uvPt.param );
 558         uvPt.u = p.X();
 559         uvPt.v = p.Y();
 560       }
 561       else {
 562         uvPt.u = uvPt.v = 1e+100;
 563       }
 564     }
 565   }
 566   return myFalsePoints;
 567 }
 568
 569 //=======================================================================
 570 //function : GetOrderedNodes
 571 //purpose  : Return nodes in the order they encounter while walking along the side
 572 //=======================================================================
 573
 574 std::vector<const SMDS_MeshNode*> StdMeshers_FaceSide::GetOrderedNodes() const
 575 {
 576   vector<const SMDS_MeshNode*> resultNodes;
 577   if ( myPoints.empty() )
 578   {
 579     if ( NbEdges() == 0 ) return resultNodes;
 580
 581     SMESHDS_Mesh* meshDS = myProxyMesh->GetMeshDS();
 582     SMESH_MesherHelper helper(*myProxyMesh->GetMesh());
 583     bool paramOK;
 584
 585     // Sort nodes of all edges putting them into a map
 586
 587     map< double, const SMDS_MeshNode*> u2node;
 588     for ( int i = 0; i < myEdge.size(); ++i )
 589     {
 590       // Put 1st vertex node of a current edge
 591       TopoDS_Vertex VV[2]; // TopExp::FirstVertex() returns NULL for INTERNAL edge
 592       VV[0] = SMESH_MesherHelper::IthVertex( 0, myEdge[i]);
 593       VV[1] = SMESH_MesherHelper::IthVertex( 1, myEdge[i]);
 594       const SMDS_MeshNode* node = SMESH_Algo::VertexNode( VV[0], meshDS );
 595       double prevNormPar = ( i == 0 ? 0 : myNormPar[ i-1 ]); // normalized param
 596       if ( node ) { // nodes on internal vertices may be missing
 597         u2node.insert( make_pair( prevNormPar, node ));
 598       }
 599       else if ( i == 0 ) {
 600         MESSAGE(" NO NODE on VERTEX" );
 601         return resultNodes;
 602       }
 603
 604       // Put internal nodes
 605       if ( SMESHDS_SubMesh* sm = meshDS->MeshElements( myEdge[i] ))
 606       {
 607         SMDS_NodeIteratorPtr nItr = sm->GetNodes();
 608         double paramSize = myLast[i] - myFirst[i];
 609         double r         = myNormPar[i] - prevNormPar;
 610         helper.SetSubShape( myEdge[i] );
 611         helper.ToFixNodeParameters( true );
 612         while ( nItr->more() )
 613         {
 614           const SMDS_MeshNode* node = nItr->next();
 615           if ( myIgnoreMediumNodes && SMESH_MeshEditor::IsMedium( node, SMDSAbs_Edge ))
 616             continue;
 617           double u = helper.GetNodeU( myEdge[i], node, 0, &paramOK );
 618
 619           // paramSize is signed so orientation is taken into account
 620           double normPar = prevNormPar + r * ( u - myFirst[i] ) / paramSize;
 621           u2node.insert( u2node.end(), make_pair( normPar, node ));
 622         }
 623       }
 624
 625       // Put 2nd vertex node for a last edge
 626       if ( i+1 == myEdge.size() ) {
 627         node = SMESH_Algo::VertexNode( VV[1], meshDS );
 628         if ( !node ) {
 629           return resultNodes;
 630         }
 631         u2node.insert( u2node.end(), make_pair( 1., node ));
 632       }
 633     }
 634
 635     // Fill the result vector
 636
 637     if ( u2node.size() == myNbPonits )
 638     {
 639       resultNodes.reserve( u2node.size() );
 640       map< double, const SMDS_MeshNode*>::iterator u2n = u2node.begin();
 641       for ( ; u2n != u2node.end(); ++u2n )
 642         resultNodes.push_back( u2n->second );
 643     }
 644   }
 645   else
 646   {
 647     resultNodes.resize( myPoints.size() );
 648     for ( size_t i = 0; i < myPoints.size(); ++i )
 649       resultNodes[i] = myPoints[i].node;
 650   }
 651
 652   return resultNodes;
 653 }
 654
 655 //================================================================================
 656 /*!
 657  * \brief reverse order of vector elements
 658   * \param vec - vector to reverse
 659  */
 660 //================================================================================
 661
 662 template <typename T > void reverse(vector<T> & vec)
 663 {
 664   for ( int f=0, r=vec.size()-1; f < r; ++f, --r )
 665     std::swap( vec[f], vec[r] );
 666 }
 667
 668 //================================================================================
 669 /*!
 670  * \brief Change orientation of side geometry
 671  */
 672 //================================================================================
 673
 674 void StdMeshers_FaceSide::Reverse()
 675 {
 676   int nbEdges = myEdge.size();
 677   for ( int i = nbEdges-1; i >= 0; --i ) {
 678     std::swap( myFirst[i], myLast[i] );
 679     if ( !myEdge[i].IsNull() )
 680       myEdge[i].Reverse();
 681     if ( i > 0 ) // at the first loop 1. is overwritten
 682       myNormPar[i] = 1 - myNormPar[i-1];
 683   }
 684   if ( nbEdges > 1 ) {
 685     reverse( myEdge );
 686     reverse( myEdgeID );
 687     reverse( myC2d );
 688     //reverse( myC3dAdaptor );
 689     reverse( myFirst );
 690     reverse( myLast );
 691     reverse( myNormPar );
 692     reverse( myEdgeLength );
 693     reverse( myIsUniform );
 694   }
 695   if ( nbEdges > 0 )
 696   {
 697     myNormPar[nbEdges-1]=1.;
 698     if ( !myEdge[0].IsNull() )
 699     {
 700       for ( size_t i = 0; i < myEdge.size(); ++i )
 701         reverseProxySubmesh( myEdge[i] );
 702       myPoints.clear();
 703       myFalsePoints.clear();
 704     }
 705     else
 706     {
 707       for ( size_t i = 0; i < myPoints.size(); ++i )
 708       {
 709         UVPtStruct & uvPt = myPoints[i];
 710         uvPt.normParam = 1 - uvPt.normParam;
 711         uvPt.x         = 1 - uvPt.x;
 712         uvPt.y         = 1 - uvPt.y;
 713       }
 714       reverse( myPoints );
 715
 716       for ( size_t i = 0; i < myFalsePoints.size(); ++i )
 717       {
 718         UVPtStruct & uvPt = myFalsePoints[i];
 719         uvPt.normParam = 1 - uvPt.normParam;
 720         uvPt.x         = 1 - uvPt.x;
 721         uvPt.y         = 1 - uvPt.y;
 722       }
 723       reverse( myFalsePoints );
 724     }
 725   }
 726   for ( size_t i = 0; i < myEdge.size(); ++i )
 727   {
 728     if ( myEdge[i].IsNull() ) continue; // for a side on points only
 729     double fp,lp;
 730     Handle(Geom_Curve) C3d = BRep_Tool::Curve(myEdge[i],fp,lp);
 731     if ( !C3d.IsNull() )
 732       myC3dAdaptor[i].Load( C3d, fp,lp );
 733   }
 734 }
 735
 736 //=======================================================================
 737 //function : SetIgnoreMediumNodes
 738 //purpose  : Make ignore medium nodes
 739 //=======================================================================
 740
 741 void StdMeshers_FaceSide::SetIgnoreMediumNodes(bool toIgnore)
 742 {
 743   if ( myIgnoreMediumNodes != toIgnore )
 744   {
 745     myIgnoreMediumNodes = toIgnore;
 746
 747     if ( !myPoints.empty() )
 748     {
 749       UVPtStructVec newPoints;
 750       newPoints.reserve( myPoints.size()/2 + 1 );
 751       for ( size_t i = 0; i < myPoints.size(); i += 2 )
 752         newPoints.push_back( myPoints[i] );
 753
 754       myPoints.swap( newPoints );
 755     }
 756     else
 757     {
 758       NbPoints( /*update=*/true );
 759     }
 760   }
 761 }
 762
 763 //=======================================================================
 764 //function : NbPoints
 765 //purpose  : Return nb nodes on edges and vertices (+1 to be == GetUVPtStruct().size() )
 766 //           Call it with update == true if mesh of this side can be recomputed
 767 //           since creation of this side
 768 //=======================================================================
 769
 770 int StdMeshers_FaceSide::NbPoints(const bool update) const
 771 {
 772   if ( !myPoints.empty() )
 773     return myPoints.size();
 774
 775   // if ( !myFalsePoints.empty() )
 776   //   return myFalsePoints.size();
 777
 778   if ( update && myEdge.size() > 0 )
 779   {
 780     StdMeshers_FaceSide* me = (StdMeshers_FaceSide*) this;
 781     me->myNbPonits = 0;
 782     me->myNbSegments = 0;
 783     me->myMissingVertexNodes = false;
 784
 785     for ( int i = 0; i < NbEdges(); ++i )
 786     {
 787       TopoDS_Vertex v1 = SMESH_MesherHelper::IthVertex( 0, myEdge[i] );
 788       if ( SMESH_Algo::VertexNode( v1, myProxyMesh->GetMeshDS() ))
 789         me->myNbPonits += 1; // for the first end
 790       else
 791         me->myMissingVertexNodes = true;
 792
 793       if ( const SMESHDS_SubMesh* sm = myProxyMesh->GetSubMesh( Edge(i) )) {
 794         int nbN = sm->NbNodes();
 795         if ( myIgnoreMediumNodes ) {
 796           SMDS_ElemIteratorPtr elemIt = sm->GetElements();
 797           if ( elemIt->more() && elemIt->next()->IsQuadratic() )
 798             nbN -= sm->NbElements();
 799         }
 800         me->myNbPonits   += nbN;
 801         me->myNbSegments += sm->NbElements();
 802       }
 803     }
 804     TopoDS_Vertex v1 = SMESH_MesherHelper::IthVertex( 1, Edge( NbEdges()-1 ));
 805     if ( SMESH_Algo::VertexNode( v1, myProxyMesh->GetMeshDS() ))
 806       me->myNbPonits++; // for the last end
 807     else
 808       me->myMissingVertexNodes = true;
 809   }
 810   return myNbPonits;
 811 }
 812
 813 //=======================================================================
 814 //function : NbSegments
 815 //purpose  : Return nb edges
 816 //           Call it with update == true if mesh of this side can be recomputed
 817 //           since creation of this side
 818 //=======================================================================
 819
 820 int StdMeshers_FaceSide::NbSegments(const bool update) const
 821 {
 822   return NbPoints( update ), myNbSegments;
 823 }
 824
 825 //================================================================================
 826 /*!
 827  * \brief Reverse UVPtStructVec if a proxy sub-mesh of E
 828  */
 829 //================================================================================
 830
 831 void StdMeshers_FaceSide::reverseProxySubmesh( const TopoDS_Edge& E )
 832 {
 833   if ( !myProxyMesh ) return;
 834   if ( const SMESH_ProxyMesh::SubMesh* sm = myProxyMesh->GetProxySubMesh( E ))
 835   {
 836     UVPtStructVec& edgeUVPtStruct = (UVPtStructVec& ) sm->GetUVPtStructVec();
 837     for ( size_t i = 0; i < edgeUVPtStruct.size(); ++i )
 838     {
 839       UVPtStruct & uvPt = edgeUVPtStruct[i];
 840       uvPt.normParam = 1 - uvPt.normParam;
 841       uvPt.x         = 1 - uvPt.x;
 842       uvPt.y         = 1 - uvPt.y;
 843     }
 844     reverse( edgeUVPtStruct );
 845   }
 846 }
 847
 848 //================================================================================
 849 /*!
 850  * \brief Show side features
 851  */
 852 //================================================================================
 853
 854 void StdMeshers_FaceSide::dump(const char* msg) const
 855 {
 856 #ifdef _DEBUG_
 857   if (msg) MESSAGE ( std::endl << msg );
 858   MESSAGE_BEGIN ("NB EDGES: "<< myEdge.size() );
 859   MESSAGE_ADD ( "nbPoints: "<<myNbPonits<<" vecSize: " << myPoints.size()<<" "<<myFalsePoints.size() );
 860   for ( int i=0; i<myEdge.size(); ++i)
 861   {
 862     MESSAGE_ADD ( "\t"<<i+1 );
 863     MESSAGE_ADD ( "\tEDGE: " );
 864     if (myEdge[i].IsNull()) {
 865       MESSAGE_ADD ( "NULL" );
 866     }
 867     else {
 868       TopAbs::Print(myEdge[i].Orientation(),cout)<<" "<<myEdge[i].TShape().operator->()<<endl;
 869       MESSAGE_ADD ( "\tV1: " << TopExp::FirstVertex( myEdge[i], 1).TShape().operator->()
 870                  << "  V2: " << TopExp::LastVertex( myEdge[i], 1).TShape().operator->() );
 871     }
 872     MESSAGE_ADD ( "\tC2d: ");
 873
 874     if (myC2d[i].IsNull()) {
 875       MESSAGE_ADD ( "NULL" );
 876     }
 877     else {
 878       MESSAGE_ADD ( myC2d[i].operator->() );
 879     }
 880
 881     MESSAGE_ADD ( "\tF: "<<myFirst[i]<< " L: "<< myLast[i] );
 882     MESSAGE_END ( "\tnormPar: "<<myNormPar[i]<<endl );
 883   }
 884 #endif
 885 }
 886
 887 //================================================================================
 888 /*!
 889  * \brief Creates a Adaptor2d_Curve2d to be used in SMESH_Block
 890   * \retval Adaptor2d_Curve2d* -
 891  */
 892 //================================================================================
 893
 894 struct Adaptor2dCurve2d : public Adaptor2d_Curve2d
 895 {
 896   const StdMeshers_FaceSide* mySide;
 897   Adaptor2dCurve2d(const StdMeshers_FaceSide* faceSide):mySide(faceSide) {}
 898   gp_Pnt2d Value(const Standard_Real U) const { return mySide->Value2d( U ); }
 899   Standard_Real FirstParameter() const { return 0; }
 900   Standard_Real LastParameter() const { return 1; }
 901 };
 902
 903 Adaptor2d_Curve2d* StdMeshers_FaceSide::GetCurve2d() const
 904 {
 905   return new Adaptor2dCurve2d( this );
 906 }
 907
 908 //================================================================================
 909 /*!
 910  * \brief Creates a fully functional Adaptor_Curve
 911  */
 912 //================================================================================
 913
 914 BRepAdaptor_CompCurve* StdMeshers_FaceSide::GetCurve3d() const
 915 {
 916   if ( myEdge.empty() )
 917     return 0;
 918
 919   TopoDS_Wire aWire;
 920   BRep_Builder aBuilder;
 921   aBuilder.MakeWire(aWire);
 922   for ( int i=0; i<myEdge.size(); ++i )
 923     aBuilder.Add( aWire, myEdge[i] );
 924
 925   if ( myEdge.size() == 2 && FirstVertex().IsSame( LastVertex() ))
 926     aWire.Closed(true); // issue 0021141
 927
 928   return new BRepAdaptor_CompCurve( aWire );
 929 }
 930
 931 //================================================================================
 932 /*!
 933  * \brief Return 2D point by normalized parameter
 934   * \param U - normalized parameter value
 935   * \retval gp_Pnt2d - point
 936  */
 937 //================================================================================
 938
 939 gp_Pnt2d StdMeshers_FaceSide::Value2d(double U) const
 940 {
 941   if ( !myC2d[0].IsNull() ) {
 942     int i = EdgeIndex( U );
 943     double prevU = i ? myNormPar[ i-1 ] : 0;
 944     double r = ( U - prevU )/ ( myNormPar[ i ] - prevU );
 945
 946     double par = myFirst[i] * ( 1 - r ) + myLast[i] * r;
 947
 948     // check parametrization of curve
 949     if( !myIsUniform[i] )
 950     {
 951       double aLen3dU = r * myEdgeLength[i] * ( myFirst[i]>myLast[i] ? -1. : 1.);
 952       GCPnts_AbscissaPoint AbPnt
 953         ( const_cast<GeomAdaptor_Curve&>( myC3dAdaptor[i]), aLen3dU, myFirst[i] );
 954       if( AbPnt.IsDone() ) {
 955         par = AbPnt.Parameter();
 956       }
 957     }
 958     return myC2d[ i ]->Value(par);
 959
 960   }
 961   else if ( !myPoints.empty() )
 962   {
 963     int i = U * double( myPoints.size()-1 );
 964     while ( i > 0 && myPoints[ i ].normParam > U )
 965       --i;
 966     while ( i+1 < myPoints.size() && myPoints[ i+1 ].normParam < U )
 967       ++i;
 968     double r = (( U - myPoints[ i ].normParam ) /
 969                 ( myPoints[ i+1 ].normParam - myPoints[ i ].normParam ));
 970     return ( myPoints[ i   ].UV() * ( 1 - r ) +
 971              myPoints[ i+1 ].UV() * r );
 972   }
 973   return myDefaultPnt2d;
 974 }
 975
 976 //================================================================================
 977 /*!
 978  * \brief Return XYZ by normalized parameter
 979   * \param U - normalized parameter value
 980   * \retval gp_Pnt - point
 981  */
 982 //================================================================================
 983
 984 gp_Pnt StdMeshers_FaceSide::Value3d(double U) const
 985 {
 986   int        i = EdgeIndex( U );
 987   double prevU = i ? myNormPar[ i-1 ] : 0;
 988   double     r = ( U - prevU )/ ( myNormPar[ i ] - prevU );
 989
 990   double par = myFirst[i] * ( 1 - r ) + myLast[i] * r;
 991
 992   // check parametrization of curve
 993   if( !myIsUniform[i] )
 994   {
 995     double aLen3dU = r * myEdgeLength[i] * ( myFirst[i]>myLast[i] ? -1. : 1.);
 996     GCPnts_AbscissaPoint AbPnt
 997       ( const_cast<GeomAdaptor_Curve&>( myC3dAdaptor[i]), aLen3dU, myFirst[i] );
 998     if( AbPnt.IsDone() ) {
 999       par = AbPnt.Parameter();
1000     }
1001   }
1002   return myC3dAdaptor[ i ].Value(par);
1003 }
1004
1005 //================================================================================
1006 /*!
1007  * \brief Return wires of a face as StdMeshers_FaceSide's
1008  */
1009 //================================================================================
1010
1011 TSideVector StdMeshers_FaceSide::GetFaceWires(const TopoDS_Face&   theFace,
1012                                               SMESH_Mesh &         theMesh,
1013                                               const bool           theIgnoreMediumNodes,
1014                                               TError &             theError,
1015                                               SMESH_ProxyMesh::Ptr theProxyMesh)
1016 {
1017   list< TopoDS_Edge > edges, internalEdges;
1018   list< int > nbEdgesInWires;
1019   int nbWires = SMESH_Block::GetOrderedEdges (theFace, edges, nbEdgesInWires);
1020
1021   // split list of all edges into separate wires
1022   TSideVector wires( nbWires );
1023   list< int >::iterator nbE = nbEdgesInWires.begin();
1024   list< TopoDS_Edge >::iterator from = edges.begin(), to = from;
1025   for ( int iW = 0; iW < nbWires; ++iW, ++nbE )
1026   {
1027     std::advance( to, *nbE );
1028     if ( *nbE == 0 ) // Issue 0020676
1029     {
1030       --nbWires;
1031       --iW;
1032       wires.resize( nbWires );
1033       continue;
1034     }
1035     list< TopoDS_Edge > wireEdges( from, to );
1036     // assure that there is a node on the first vertex
1037     // as StdMeshers_FaceSide::GetUVPtStruct() requires
1038     if ( wireEdges.front().Orientation() != TopAbs_INTERNAL ) // Issue 0020676
1039     {
1040       while ( !SMESH_Algo::VertexNode( TopExp::FirstVertex( wireEdges.front(), true),
1041                                        theMesh.GetMeshDS()))
1042       {
1043         wireEdges.splice(wireEdges.end(), wireEdges,
1044                          wireEdges.begin(), ++wireEdges.begin());
1045         if ( from->IsSame( wireEdges.front() )) {
1046           theError = TError
1047             ( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"No nodes on vertices"));
1048           return TSideVector(0);
1049         }
1050       }
1051     }
1052     else if ( *nbE > 1 ) // Issue 0020676 (Face_pb_netgen.brep) - several internal edges in a wire
1053     {
1054       internalEdges.splice( internalEdges.end(), wireEdges, ++wireEdges.begin(), wireEdges.end());
1055     }
1056
1057     StdMeshers_FaceSide* wire = new StdMeshers_FaceSide( theFace, wireEdges, &theMesh,
1058                                                          /*isForward=*/true, theIgnoreMediumNodes,
1059                                                          theProxyMesh );
1060     wires[ iW ] = StdMeshers_FaceSidePtr( wire );
1061     from = to;
1062   }
1063   while ( !internalEdges.empty() )
1064   {
1065     StdMeshers_FaceSide* wire = new StdMeshers_FaceSide( theFace, internalEdges.back(), &theMesh,
1066                                                          /*isForward=*/true, theIgnoreMediumNodes,
1067                                                          theProxyMesh );
1068     wires.push_back( StdMeshers_FaceSidePtr( wire ));
1069     internalEdges.pop_back();
1070   }
1071   return wires;
1072 }
1073
1074 //================================================================================
1075 /*!
1076  * \brief Return 1st vertex of the i-the edge
1077  */
1078 //================================================================================
1079
1080 TopoDS_Vertex StdMeshers_FaceSide::FirstVertex(int i) const
1081 {
1082   TopoDS_Vertex v;
1083   if ( i < NbEdges() )
1084   {
1085     v = myEdge[i].Orientation() <= TopAbs_REVERSED ? // FORWARD || REVERSED
1086         TopExp::FirstVertex( myEdge[i], 1 )        :
1087         TopoDS::Vertex( TopoDS_Iterator( myEdge[i] ).Value() );
1088   }
1089   return v;
1090 }
1091
1092 //================================================================================
1093 /*!
1094  * \brief Return last vertex of the i-the edge
1095  */
1096 //================================================================================
1097
1098 TopoDS_Vertex StdMeshers_FaceSide::LastVertex(int i) const
1099 {
1100   TopoDS_Vertex v;
1101   if ( i < NbEdges() )
1102   {
1103     const TopoDS_Edge& edge = i<0 ? myEdge[ NbEdges() + i ] : myEdge[i];
1104     if ( edge.Orientation() <= TopAbs_REVERSED ) // FORWARD || REVERSED
1105       v = TopExp::LastVertex( edge, 1 );
1106     else
1107       for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
1108         v = TopoDS::Vertex( vIt.Value() );
1109   }
1110   return v;
1111 }
1112