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