src/SMESHUtils/SMESH_PolyLine.cxx

   1 // Copyright (C) 2018-2019  OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // File      : SMESH_PolyLine.cxx
  20 // Created   : Thu Dec  6 17:33:26 2018
  21 // Author    : Edward AGAPOV (eap)
  22
  23 #include "SMESH_MeshAlgos.hxx"
  24
  25 #include "SMDS_MeshGroup.hxx"
  26 #include "SMDS_LinearEdge.hxx"
  27 #include "SMDS_Mesh.hxx"
  28 #include "SMESH_TryCatch.hxx"
  29
  30 #include <OSD_Parallel.hxx>
  31 #include <Precision.hxx>
  32
  33 namespace
  34 {
  35   //================================================================================
  36   /*!
  37    * \brief Sequence of found points and a current point data
  38    */
  39   struct Path
  40   {
  41     std::vector< gp_XYZ >   myPoints;
  42     double                  myLength;
  43
  44     const SMDS_MeshElement* myFace;
  45     SMESH_NodeXYZ           myNode1; // nodes of the edge the path entered myFace
  46     SMESH_NodeXYZ           myNode2;
  47     int                     myNodeInd1;
  48     int                     myNodeInd2;
  49     double                  myDot1;
  50     double                  myDot2;
  51
  52     int                     mySrcPntInd; //!< start point index
  53     TIDSortedElemSet        myElemSet, myAvoidSet;
  54
  55     Path(): myLength(0.0), myFace(0) {}
  56
  57     bool SetCutAtCorner( const SMESH_NodeXYZ&    cornerNode,
  58                          const SMDS_MeshElement* face,
  59                          const gp_XYZ&           plnNorm,
  60                          const gp_XYZ&           plnOrig );
  61
  62     void AddPoint( const gp_XYZ& p );
  63
  64     bool Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig );
  65
  66     bool ReachSamePoint( const Path& other );
  67
  68     static void Remove( std::vector< Path > & paths, size_t& i );
  69   };
  70
  71   //================================================================================
  72   /*!
  73    * \brief Return true if this Path meats another
  74    */
  75   //================================================================================
  76
  77   bool Path::ReachSamePoint( const Path& other )
  78   {
  79     return ( mySrcPntInd != other.mySrcPntInd &&
  80              myFace == other.myFace );
  81   }
  82
  83   //================================================================================
  84   /*!
  85    * \brief Remove a path from a vector
  86    */
  87   //================================================================================
  88
  89   void Path::Remove( std::vector< Path > & paths, size_t& i )
  90   {
  91     if ( paths.size() > 1 )
  92     {
  93       size_t j = paths.size() - 1; // last item to be removed
  94       if ( i < j )
  95       {
  96         paths[ i ].myPoints.swap( paths[ j ].myPoints );
  97         paths[ i ].myLength    = paths[ j ].myLength;
  98         paths[ i ].mySrcPntInd = paths[ j ].mySrcPntInd;
  99         paths[ i ].myFace      = paths[ j ].myFace;
 100         paths[ i ].myNode1     = paths[ j ].myNode1;
 101         paths[ i ].myNode2     = paths[ j ].myNode2;
 102         paths[ i ].myNodeInd1  = paths[ j ].myNodeInd1;
 103         paths[ i ].myNodeInd2  = paths[ j ].myNodeInd2;
 104         paths[ i ].myDot1      = paths[ j ].myDot1;
 105         paths[ i ].myDot2      = paths[ j ].myDot2;
 106       }
 107     }
 108     paths.pop_back();
 109     if ( i > 0 )
 110       --i;
 111   }
 112
 113   //================================================================================
 114   /*!
 115    * \brief Store a point that is at a node of a face if the face is intersected by plane.
 116    *        Return false if the node is a sole intersection point of the face and the plane
 117    */
 118   //================================================================================
 119
 120   bool Path::SetCutAtCorner( const SMESH_NodeXYZ&    cornerNode,
 121                              const SMDS_MeshElement* face,
 122                              const gp_XYZ&           plnNorm,
 123                              const gp_XYZ&           plnOrig )
 124   {
 125     if ( face == myFace )
 126       return false;
 127     myNodeInd1 = face->GetNodeIndex( cornerNode._node );
 128     myNodeInd2 = ( myNodeInd1 + 1 ) % face->NbCornerNodes();
 129     int   ind3 = ( myNodeInd1 + 2 ) % face->NbCornerNodes();
 130     myNode1.Set( face->GetNode( ind3 ));
 131     myNode2.Set( face->GetNode( myNodeInd2 ));
 132
 133     myDot1 = plnNorm * ( myNode1 - plnOrig );
 134     myDot2 = plnNorm * ( myNode2 - plnOrig );
 135
 136     bool ok = ( myDot1 * myDot2 < 0 );
 137     if ( !ok && myDot1 * myDot2 == 0 )
 138     {
 139       ok = ( myDot1 != myDot2 );
 140       if ( ok && myFace )
 141         ok = ( myFace->GetNodeIndex(( myDot1 == 0 ? myNode1 : myNode2 )._node ) < 0 );
 142     }
 143     if ( ok )
 144     {
 145       myFace = face;
 146       myDot1 = 0;
 147       AddPoint( cornerNode );
 148     }
 149     return ok;
 150   }
 151
 152   //================================================================================
 153   /*!
 154    * \brief Store a point and update myLength
 155    */
 156   //================================================================================
 157
 158   void Path::AddPoint( const gp_XYZ& p )
 159   {
 160     if ( !myPoints.empty() )
 161       myLength += ( p - myPoints.back() ).Modulus();
 162     else
 163       myLength = 0;
 164     myPoints.push_back( p );
 165   }
 166
 167   //================================================================================
 168   /*!
 169    * \brief Try to find the next point
 170    *  \param [in] plnNorm - cutting plane normal
 171    *  \param [in] plnOrig - cutting plane origin
 172    */
 173   //================================================================================
 174
 175   bool Path::Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig )
 176   {
 177     int nodeInd3 = ( myNodeInd1 + 1 ) % myFace->NbCornerNodes();
 178     if ( myNodeInd2 == nodeInd3 )
 179       nodeInd3 = ( myNodeInd1 + 2 ) % myFace->NbCornerNodes();
 180
 181     SMESH_NodeXYZ node3 = myFace->GetNode( nodeInd3 );
 182     double         dot3 = plnNorm * ( node3 - plnOrig );
 183
 184     if ( dot3 * myDot1 < 0. )
 185     {
 186       myNode2    = node3;
 187       myNodeInd2 = nodeInd3;
 188       myDot2     = dot3;
 189     }
 190     else if ( dot3 * myDot2 < 0. )
 191     {
 192       myNode1    = node3;
 193       myNodeInd1 = nodeInd3;
 194       myDot1     = dot3;
 195     }
 196     else if ( dot3 == 0. )
 197     {
 198       SMDS_ElemIteratorPtr fIt = node3._node->GetInverseElementIterator(SMDSAbs_Face);
 199       while ( fIt->more() )
 200         if ( SetCutAtCorner( node3, fIt->next(), plnNorm, plnOrig ))
 201           return true;
 202       return false;
 203     }
 204     else if ( myDot2 == 0. )
 205     {
 206       SMESH_NodeXYZ node2 = myNode2; // copy as myNode2 changes in SetCutAtCorner()
 207       SMDS_ElemIteratorPtr fIt = node2._node->GetInverseElementIterator(SMDSAbs_Face);
 208       while ( fIt->more() )
 209         if ( SetCutAtCorner( node2, fIt->next(), plnNorm, plnOrig ))
 210           return true;
 211       return false;
 212     }
 213
 214     double r = Abs( myDot1 / ( myDot2 - myDot1 ));
 215     AddPoint( myNode1 * ( 1 - r ) + myNode2 * r );
 216
 217     myAvoidSet.clear();
 218     myAvoidSet.insert( myFace );
 219     myFace = SMESH_MeshAlgos::FindFaceInSet( myNode1._node, myNode2._node,
 220                                              myElemSet,   myAvoidSet,
 221                                              &myNodeInd1, &myNodeInd2 );
 222     return myFace;
 223   }
 224
 225   //================================================================================
 226   /*!
 227    * \brief Compute a path between two points of PolySegment
 228    */
 229   struct PolyPathCompute
 230   {
 231     SMESH_MeshAlgos::TListOfPolySegments& mySegments; //!< inout PolySegment's
 232     std::vector< Path >&                  myPaths;    //!< path of each of segments to compute
 233     SMDS_Mesh*                            myMesh;
 234     mutable std::vector< std::string >    myErrors;
 235
 236     PolyPathCompute( SMESH_MeshAlgos::TListOfPolySegments& theSegments,
 237                      std::vector< Path >&                  thePaths,
 238                      SMDS_Mesh*                            theMesh):
 239       mySegments( theSegments ),
 240       myPaths( thePaths ),
 241       myMesh( theMesh ),
 242       myErrors( theSegments.size() )
 243     {
 244     }
 245
 246 #undef SMESH_CAUGHT
 247 #define SMESH_CAUGHT myErrors[i] =
 248     void operator() ( const int i ) const
 249     {
 250       SMESH_TRY;
 251       const_cast< PolyPathCompute* >( this )->Compute( i );
 252       SMESH_CATCH( SMESH::returnError );
 253     }
 254 #undef SMESH_CAUGHT
 255
 256     //================================================================================
 257     /*!
 258      * \brief Compute a path of a given segment
 259      */
 260     //================================================================================
 261
 262     void Compute( const int iSeg )
 263     {
 264       SMESH_MeshAlgos::PolySegment& polySeg = mySegments[ iSeg ];
 265
 266       // the cutting plane
 267       gp_XYZ plnNorm = ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ) ^ polySeg.myVector.XYZ();
 268       gp_XYZ plnOrig = polySeg.myXYZ[1];
 269
 270       // Find paths connecting the 2 end points of polySeg
 271
 272       std::vector< Path > paths; paths.reserve(10);
 273
 274       // 1) initialize paths; two paths starts at each end point
 275
 276       for ( int iP = 0; iP < 2; ++iP ) // loop on the polySeg end points
 277       {
 278         Path path;
 279         path.mySrcPntInd = iP;
 280         size_t nbPaths = paths.size();
 281
 282         if ( polySeg.myFace[ iP ]) // the end point lies on polySeg.myFace[ iP ]
 283         {
 284           // check coincidence of polySeg.myXYZ[ iP ] with nodes
 285           const double tol = 1e-20;
 286           SMESH_NodeXYZ nodes[4];
 287           for ( int i = 0; i < 3 &&  !polySeg.myNode1[ iP ]; ++i )
 288           {
 289             nodes[ i ] = polySeg.myFace[ iP ]->GetNode( i );
 290             if (( nodes[ i ] - polySeg.myXYZ[ iP ]).SquareModulus() < tol*tol )
 291               polySeg.myNode1[ iP ] = nodes[ i ].Node();
 292           }
 293           nodes[ 3 ] = nodes[ 0 ];
 294
 295           // check coincidence of polySeg.myXYZ[ iP ] with edges
 296           for ( int i = 0; i < 3 &&  !polySeg.myNode1[ iP ]; ++i )
 297           {
 298             SMDS_LinearEdge edge( nodes[i].Node(), nodes[i+1].Node() );
 299             if ( SMESH_MeshAlgos::GetDistance( &edge, polySeg.myXYZ[ iP ]) < tol )
 300             {
 301               polySeg.myNode1[ iP ] = nodes[ i     ].Node();
 302               polySeg.myNode2[ iP ] = nodes[ i + 1 ].Node();
 303             }
 304           }
 305
 306           if ( !polySeg.myNode1[ iP ] ) // polySeg.myXYZ[ iP ] is within polySeg.myFace[ iP ]
 307           {
 308             double dot[ 4 ];
 309             for ( int i = 0; i < 3; ++i )
 310               dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
 311             dot[ 3 ] = dot[ 0 ];
 312
 313             int iCut = 0; // index of a cut edge
 314             if      ( dot[ 1 ] * dot[ 2 ] < 0. ) iCut = 1;
 315             else if ( dot[ 2 ] * dot[ 3 ] < 0. ) iCut = 2;
 316
 317             // initialize path so as if it entered the face via iCut-th edge
 318             path.myFace = polySeg.myFace[ iP ];
 319             path.myNodeInd1 = iCut;
 320             path.myNodeInd2 = iCut + 1;
 321             path.myNode1.Set( nodes[ iCut     ].Node() );
 322             path.myNode2.Set( nodes[ iCut + 1 ].Node() );
 323             path.myDot1 = dot[ iCut ];
 324             path.myDot2 = dot[ iCut + 1 ];
 325             path.myPoints.clear();
 326             path.AddPoint( polySeg.myXYZ[ iP ]);
 327             paths.push_back( path );
 328
 329             path.Extend( plnNorm, plnOrig ); // to get another edge cut
 330             path.myFace = polySeg.myFace[ iP ];
 331             if ( path.myDot1 == 0. ) // cut at a node
 332             {
 333               path.myNodeInd1 = ( iCut + 2 ) % 3;
 334               path.myNodeInd2 = ( iCut + 3 ) % 3;
 335               path.myNode2.Set( path.myFace->GetNode( path.myNodeInd2 ));
 336               path.myDot2 = dot[ path.myNodeInd2 ];
 337             }
 338             else
 339             {
 340               path.myNodeInd1 = path.myFace->GetNodeIndex( path.myNode1.Node() );
 341               path.myNodeInd2 = path.myFace->GetNodeIndex( path.myNode2.Node() );
 342             }
 343             path.myPoints.clear();
 344             path.AddPoint( polySeg.myXYZ[ iP ]);
 345             paths.push_back( path );
 346           }
 347         }
 348
 349         if ( polySeg.myNode2[ iP ] && polySeg.myNode2[ iP ] != polySeg.myNode1[ iP ] )
 350         {
 351           // the end point is on an edge
 352           while (( path.myFace = SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
 353                                                                  polySeg.myNode2[ iP ],
 354                                                                  path.myElemSet,
 355                                                                  path.myAvoidSet,
 356                                                                  &path.myNodeInd1,
 357                                                                  &path.myNodeInd2 )))
 358           {
 359             path.myNode1.Set( polySeg.myNode1[ iP ]);
 360             path.myNode2.Set( polySeg.myNode2[ iP ]);
 361             path.myDot1 = plnNorm * ( path.myNode1 - plnOrig );
 362             path.myDot2 = plnNorm * ( path.myNode2 - plnOrig );
 363             path.myPoints.clear();
 364             path.AddPoint( polySeg.myXYZ[ iP ]);
 365             path.myAvoidSet.insert( path.myFace );
 366             paths.push_back( path );
 367           }
 368           if ( nbPaths == paths.size() )
 369             throw SALOME_Exception ( SMESH_Comment("No face edge found by point ") << iP+1
 370                                      << " in a PolySegment " << iSeg );
 371         }
 372         else if ( polySeg.myNode1[ iP ] ) // the end point is at a node
 373         {
 374           std::set<const SMDS_MeshNode* > nodes;
 375           SMDS_ElemIteratorPtr fIt = polySeg.myNode1[ iP ]->GetInverseElementIterator(SMDSAbs_Face);
 376           while ( fIt->more() )
 377           {
 378             path.myPoints.clear();
 379             if ( path.SetCutAtCorner( polySeg.myNode1[ iP ], fIt->next(), plnNorm, plnOrig ))
 380             {
 381               if (( path.myDot1 * path.myDot2 != 0 ) ||
 382                   ( nodes.insert( path.myDot1 == 0 ? path.myNode1._node : path.myNode2._node ).second ))
 383                 paths.push_back( path );
 384             }
 385           }
 386         }
 387
 388         // look for a one-segment path
 389         for ( size_t i = 0; i < nbPaths; ++i )
 390           for ( size_t j = nbPaths; j < paths.size(); ++j )
 391             if ( paths[i].myFace == paths[j].myFace )
 392             {
 393               myPaths[ iSeg ].myPoints.push_back( paths[i].myPoints[0] );
 394               myPaths[ iSeg ].myPoints.push_back( paths[j].myPoints[0] );
 395               paths.clear();
 396             }
 397       }
 398
 399       // 2) extend paths and compose the shortest one connecting the two points
 400
 401       myPaths[ iSeg ].myLength = 1e100;
 402
 403       while ( paths.size() >= 2 )
 404       {
 405         for ( size_t i = 0; i < paths.size(); ++i )
 406         {
 407           Path& path = paths[ i ];
 408           if ( !path.Extend( plnNorm, plnOrig ) ||        // path reached a mesh boundary
 409                path.myLength > myPaths[ iSeg ].myLength ) // path is longer than others
 410           {
 411             Path::Remove( paths, i );
 412             continue;
 413           }
 414
 415           // join paths that reach same point
 416           for ( size_t j = 0; j < paths.size(); ++j )
 417           {
 418             if ( i != j && paths[i].ReachSamePoint( paths[j] ))
 419             {
 420               double   distLast = ( paths[i].myPoints.back() - paths[j].myPoints.back() ).Modulus();
 421               double fullLength = ( paths[i].myLength + paths[j].myLength + distLast );
 422               if ( fullLength < myPaths[ iSeg ].myLength )
 423               {
 424                 myPaths[ iSeg ].myLength = fullLength;
 425                 std::vector< gp_XYZ > & allPoints = myPaths[ iSeg ].myPoints;
 426                 allPoints.swap( paths[i].myPoints );
 427                 allPoints.insert( allPoints.end(),
 428                                   paths[j].myPoints.rbegin(),
 429                                   paths[j].myPoints.rend() );
 430               }
 431               Path::Remove( paths, i );
 432               Path::Remove( paths, j );
 433             }
 434           }
 435         }
 436         if ( !paths.empty() && (int) paths[0].myPoints.size() > myMesh->NbFaces() )
 437           throw SALOME_Exception(LOCALIZED( "Infinite loop in MakePolyLine()"));
 438       }
 439
 440       if ( myPaths[ iSeg ].myPoints.empty() )
 441         throw SALOME_Exception( SMESH_Comment("Can't find a full path for PolySegment #") << iSeg );
 442
 443       // reverse the path
 444       double d00 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.front() ).SquareModulus();
 445       double d01 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.back()  ).SquareModulus();
 446       if ( d00 > d01 )
 447         std::reverse( myPaths[ iSeg ].myPoints.begin(), myPaths[ iSeg ].myPoints.end() );
 448
 449     } // PolyPathCompute::Compute()
 450
 451   }; // struct PolyPathCompute
 452
 453 } // namespace
 454
 455 //=======================================================================
 456 //function : MakePolyLine
 457 //purpose  : Create a polyline consisting of 1D mesh elements each lying on a 2D element of
 458 //           the initial mesh
 459 //=======================================================================
 460
 461 void SMESH_MeshAlgos::MakePolyLine( SMDS_Mesh*                            theMesh,
 462                                     TListOfPolySegments&                  theSegments,
 463                                     std::vector<const SMDS_MeshElement*>& theNewEdges,
 464                                     std::vector< const SMDS_MeshNode* >&  theNewNodes,
 465                                     SMDS_MeshGroup*                       theGroup,
 466                                     SMESH_ElementSearcher*                theSearcher)
 467 {
 468   std::vector< Path > segPaths( theSegments.size() ); // path of each of segments
 469
 470   SMESH_ElementSearcher* searcher = theSearcher;
 471   SMESHUtils::Deleter<SMESH_ElementSearcher> delSearcher;
 472   if ( !searcher )
 473   {
 474     searcher = SMESH_MeshAlgos::GetElementSearcher( *theMesh );
 475     delSearcher._obj = searcher;
 476   }
 477
 478   // get cutting planes
 479
 480   std::vector< bool > isVectorOK( theSegments.size(), true );
 481   const double planarCoef = 0.333; // plane height in planar case
 482
 483   for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
 484   {
 485     PolySegment& polySeg = theSegments[ iSeg ];
 486
 487     gp_XYZ p1 = SMESH_NodeXYZ( polySeg.myNode1[0] );
 488     gp_XYZ p2 = SMESH_NodeXYZ( polySeg.myNode1[1] );
 489     if ( polySeg.myNode2[0] ) p1 = 0.5 * ( p1 + SMESH_NodeXYZ( polySeg.myNode2[0] ));
 490     if ( polySeg.myNode2[1] ) p2 = 0.5 * ( p2 + SMESH_NodeXYZ( polySeg.myNode2[1] ));
 491
 492     polySeg.myFace[0] = polySeg.myFace[1] = 0;
 493     if ( !polySeg.myNode1[0] && !polySeg.myNode2[0] )
 494     {
 495       p1 = searcher->Project( polySeg.myXYZ[0], SMDSAbs_Face, &polySeg.myFace[0] );
 496     }
 497     if ( !polySeg.myNode1[1] && !polySeg.myNode2[1] )
 498     {
 499       p2 = searcher->Project( polySeg.myXYZ[1], SMDSAbs_Face, &polySeg.myFace[1] );
 500     }
 501     polySeg.myXYZ[0] = p1;
 502     polySeg.myXYZ[1] = p2;
 503
 504     gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
 505
 506     isVectorOK[ iSeg ] = ( plnNorm.Modulus() > std::numeric_limits<double>::min() );
 507     if ( !isVectorOK[ iSeg ])
 508     {
 509       gp_XYZ pMid = 0.5 * ( p1 + p2 );
 510       const SMDS_MeshElement* face;
 511       polySeg.myMidProjPoint = searcher->Project( pMid, SMDSAbs_Face, &face );
 512       polySeg.myVector       = polySeg.myMidProjPoint.XYZ() - pMid;
 513
 514       gp_XYZ faceNorm;
 515       SMESH_MeshAlgos::FaceNormal( face, faceNorm );
 516
 517       if ( polySeg.myVector.Magnitude() < Precision::Confusion() ||
 518            polySeg.myVector * faceNorm  < Precision::Confusion() )
 519       {
 520         polySeg.myVector = faceNorm;
 521         polySeg.myMidProjPoint = pMid + faceNorm * ( p1 - p2 ).Modulus() * planarCoef;
 522       }
 523     }
 524     else
 525     {
 526       polySeg.myVector = plnNorm ^ ( p1 - p2 );
 527     }
 528   }
 529
 530   // assure that inverse elements are constructed, avoid their concurrent building in threads
 531   theMesh->nodesIterator()->next()->NbInverseElements();
 532
 533   // find paths
 534
 535   PolyPathCompute algo( theSegments, segPaths, theMesh );
 536   OSD_Parallel::For( 0, theSegments.size(), algo, theSegments.size() == 1 );
 537
 538   for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
 539     if ( !algo.myErrors[ iSeg ].empty() )
 540       throw SALOME_Exception( algo.myErrors[ iSeg ].c_str() );
 541
 542   // create an 1D mesh
 543
 544   const SMDS_MeshNode *n, *nPrev = 0;
 545
 546   for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
 547   {
 548     const Path& path = segPaths[iSeg];
 549     if ( path.myPoints.size() < 2 )
 550       continue;
 551
 552     double tol = path.myLength / path.myPoints.size() / 1000.;
 553     if ( !nPrev || ( SMESH_NodeXYZ( nPrev ) - path.myPoints[0] ).SquareModulus() > tol*tol )
 554     {
 555       nPrev = theMesh->AddNode( path.myPoints[0].X(), path.myPoints[0].Y(), path.myPoints[0].Z() );
 556       theNewNodes.push_back( nPrev );
 557     }
 558     for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
 559     {
 560       n = theMesh->AddNode( path.myPoints[iP].X(), path.myPoints[iP].Y(), path.myPoints[iP].Z() );
 561       theNewNodes.push_back( n );
 562
 563       const SMDS_MeshElement* elem = theMesh->AddEdge( nPrev, n );
 564       theNewEdges.push_back( elem );
 565       if ( theGroup )
 566         theGroup->Add( elem );
 567
 568       nPrev = n;
 569     }
 570
 571     // return a vector
 572
 573     gp_XYZ pMid = 0.5 * ( path.myPoints[0] + path.myPoints.back() );
 574     if ( isVectorOK[ iSeg ])
 575     {
 576       // find the most distant point of a path
 577       double maxDist = 0;
 578       for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
 579       {
 580         double dist = Abs( theSegments[iSeg].myVector * ( path.myPoints[iP] - path.myPoints[0] ));
 581         if ( dist > maxDist )
 582         {
 583           maxDist = dist;
 584           theSegments[iSeg].myMidProjPoint = path.myPoints[iP];
 585         }
 586       }
 587       if ( maxDist < Precision::Confusion() ) // planar case
 588         theSegments[iSeg].myMidProjPoint =
 589           pMid + theSegments[iSeg].myVector.XYZ().Normalized() * path.myLength * planarCoef;
 590     }
 591     theSegments[iSeg].myVector = gp_Vec( pMid, theSegments[iSeg].myMidProjPoint );
 592   }
 593
 594   return;
 595 }