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(const SMDS_MeshElement* face=0, int srcInd=-1):
  56       myLength(0.0), myFace(face), mySrcPntInd( srcInd ) {}
  57
  58     void CopyNoPoints( const Path& other );
  59
  60     bool Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig, std::vector< Path > * paths = 0 );
  61
  62     bool SetCutAtCorner( const SMESH_NodeXYZ&    cornerNode,
  63                          const gp_XYZ&           plnNorm,
  64                          const gp_XYZ&           plnOrig,
  65                          std::vector< Path >*    paths);
  66
  67     bool SetCutAtCorner( const SMESH_NodeXYZ&    cornerNode,
  68                          const SMDS_MeshElement* face,
  69                          const gp_XYZ&           plnNorm,
  70                          const gp_XYZ&           plnOrig );
  71
  72     void AddPoint( const gp_XYZ& p );
  73
  74     bool ReachSamePoint( const Path& other );
  75
  76     static void Remove( std::vector< Path > & paths, size_t& i );
  77   };
  78
  79   //================================================================================
  80   /*!
  81    * \brief Return true if this Path meats another
  82    */
  83   //================================================================================
  84
  85   bool Path::ReachSamePoint( const Path& other )
  86   {
  87     return ( mySrcPntInd != other.mySrcPntInd &&
  88              myFace == other.myFace );
  89   }
  90
  91   //================================================================================
  92   /*!
  93    * \brief Copy data except points
  94    */
  95   //================================================================================
  96
  97   void Path::CopyNoPoints( const Path& other )
  98   {
  99     myLength    = other.myLength;
 100     mySrcPntInd = other.mySrcPntInd;
 101     myFace      = other.myFace;
 102     myNode1     = other.myNode1;
 103     myNode2     = other.myNode2;
 104     myNodeInd1  = other.myNodeInd1;
 105     myNodeInd2  = other.myNodeInd2;
 106     myDot1      = other.myDot1;
 107     myDot2      = other.myDot2;
 108   }
 109
 110   //================================================================================
 111   /*!
 112    * \brief Remove a path from a vector
 113    */
 114   //================================================================================
 115
 116   void Path::Remove( std::vector< Path > & paths, size_t& i )
 117   {
 118     if ( paths.size() > 1 )
 119     {
 120       size_t j = paths.size() - 1; // last item to be removed
 121       if ( i < j )
 122       {
 123         paths[ i ].CopyNoPoints ( paths[ j ]);
 124         paths[ i ].myPoints.swap( paths[ j ].myPoints );
 125       }
 126     }
 127     paths.pop_back();
 128     if ( i > 0 )
 129       --i;
 130   }
 131
 132   //================================================================================
 133   /*!
 134    * \brief Try to extend self by a point located at a node.
 135    *        Return a success flag.
 136    */
 137   //================================================================================
 138
 139   bool Path::SetCutAtCorner( const SMESH_NodeXYZ&  cornerNode,
 140                              const gp_XYZ&         plnNorm,
 141                              const gp_XYZ&         plnOrig,
 142                              std::vector< Path > * paths )
 143   {
 144     bool ok = false;
 145     const bool isContinuation = myFace; // extend this path or find all possible paths?
 146     const SMDS_MeshElement* lastFace = myFace;
 147
 148     myAvoidSet.clear();
 149
 150     SMDS_ElemIteratorPtr fIt = cornerNode->GetInverseElementIterator(SMDSAbs_Face);
 151     while ( fIt->more() )
 152     {
 153       Path path( lastFace, mySrcPntInd );
 154       if ( !path.SetCutAtCorner( cornerNode, fIt->next(), plnNorm, plnOrig ))
 155         continue;
 156
 157       if ( path.myDot1 == 0 &&
 158            !myAvoidSet.insert( path.myNode1.Node() ).second )
 159         continue;
 160       if ( path.myDot2 == 0 &&
 161            !myAvoidSet.insert( path.myNode2.Node() ).second )
 162         continue;
 163
 164       if ( isContinuation )
 165       {
 166         if ( ok ) // non-manifold continuation
 167         {
 168           path.myPoints = myPoints;
 169           path.myLength = myLength;
 170           path.AddPoint( cornerNode );
 171           paths->push_back( path );
 172         }
 173         else
 174         {
 175           double len = myLength;
 176           this->CopyNoPoints( path );
 177           this->myLength = len;
 178           this->AddPoint( path.myPoints.back() );
 179         }
 180       }
 181       else
 182       {
 183         paths->push_back( path );
 184       }
 185       ok = true;
 186     }
 187
 188     return ok;
 189   }
 190
 191   //================================================================================
 192   /*!
 193    * \brief Store a point that is at a node of a face if the face is intersected by plane.
 194    *        Return false if the node is a sole intersection point of the face and the plane
 195    */
 196   //================================================================================
 197
 198   bool Path::SetCutAtCorner( const SMESH_NodeXYZ&    cornerNode,
 199                              const SMDS_MeshElement* face,
 200                              const gp_XYZ&           plnNorm,
 201                              const gp_XYZ&           plnOrig )
 202   {
 203     if ( face == myFace )
 204       return false;
 205     myNodeInd1 = face->GetNodeIndex( cornerNode._node );
 206     myNodeInd2 = ( myNodeInd1 + 1 ) % face->NbCornerNodes();
 207     int   ind3 = ( myNodeInd1 + 2 ) % face->NbCornerNodes();
 208     myNode1.Set( face->GetNode( ind3 ));
 209     myNode2.Set( face->GetNode( myNodeInd2 ));
 210
 211     myDot1 = plnNorm * ( myNode1 - plnOrig );
 212     myDot2 = plnNorm * ( myNode2 - plnOrig );
 213
 214     bool ok = ( myDot1 * myDot2 < 0 );
 215     if ( !ok && myDot1 * myDot2 == 0 )
 216     {
 217       ok = ( myDot1 != myDot2 );
 218       if ( ok && myFace )
 219         ok = ( myFace->GetNodeIndex(( myDot1 == 0 ? myNode1 : myNode2 )._node ) < 0 );
 220     }
 221     if ( ok )
 222     {
 223       myFace = face;
 224       myDot1 = 0;
 225       AddPoint( cornerNode );
 226     }
 227     return ok;
 228   }
 229
 230   //================================================================================
 231   /*!
 232    * \brief Store a point and update myLength
 233    */
 234   //================================================================================
 235
 236   void Path::AddPoint( const gp_XYZ& p )
 237   {
 238     if ( !myPoints.empty() )
 239       myLength += ( p - myPoints.back() ).Modulus();
 240     else
 241       myLength = 0;
 242     myPoints.push_back( p );
 243   }
 244
 245   //================================================================================
 246   /*!
 247    * \brief Try to find the next point
 248    *  \param [in] plnNorm - cutting plane normal
 249    *  \param [in] plnOrig - cutting plane origin
 250    *  \param [in] paths - all paths
 251    */
 252   //================================================================================
 253
 254   bool Path::Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig, std::vector< Path > * paths )
 255   {
 256     bool ok = false;
 257
 258     int nodeInd3 = ( myNodeInd1 + 1 ) % myFace->NbCornerNodes();
 259     if ( myNodeInd2 == nodeInd3 )
 260       nodeInd3 = ( myNodeInd1 + 2 ) % myFace->NbCornerNodes();
 261
 262     SMESH_NodeXYZ node3 = myFace->GetNode( nodeInd3 );
 263     double         dot3 = plnNorm * ( node3 - plnOrig );
 264
 265     if ( dot3 * myDot1 < 0. )
 266     {
 267       myNode2    = node3;
 268       myNodeInd2 = nodeInd3;
 269       myDot2     = dot3;
 270     }
 271     else if ( dot3 * myDot2 < 0. )
 272     {
 273       myNode1    = node3;
 274       myNodeInd1 = nodeInd3;
 275       myDot1     = dot3;
 276     }
 277     else if ( dot3 == 0. )
 278     {
 279       ok = SetCutAtCorner( node3, plnNorm, plnOrig, paths );
 280       return ok;
 281     }
 282     else if ( myDot2 == 0. )
 283     {
 284       ok = SetCutAtCorner( myNode2, plnNorm, plnOrig, paths );
 285       return ok;
 286     }
 287
 288     double r = Abs( myDot1 / ( myDot2 - myDot1 ));
 289     AddPoint( myNode1 * ( 1 - r ) + myNode2 * r );
 290
 291     myAvoidSet.clear();
 292     myAvoidSet.insert( myFace );
 293     const SMDS_MeshElement* nextFace;
 294     int ind1, ind2;
 295     while (( nextFace = SMESH_MeshAlgos::FindFaceInSet( myNode1._node, myNode2._node,
 296                                                         myElemSet, myAvoidSet,
 297                                                         &ind1, &ind2 )))
 298     {
 299       if ( ok ) // non-manifold continuation
 300       {
 301         paths->push_back( *this );
 302         paths->back().myFace = nextFace;
 303         paths->back().myNodeInd1 = ind1;
 304         paths->back().myNodeInd2 = ind2;
 305       }
 306       else
 307       {
 308         myFace = nextFace;
 309         myNodeInd1 = ind1;
 310         myNodeInd2 = ind2;
 311       }
 312       ok = true;
 313       if ( !paths )
 314         break;
 315       myAvoidSet.insert( nextFace );
 316     }
 317
 318     return ok;
 319   }
 320
 321   //================================================================================
 322   /*!
 323    * \brief Compute a path between two points of PolySegment
 324    */
 325   struct PolyPathCompute
 326   {
 327     SMESH_MeshAlgos::TListOfPolySegments& mySegments; //!< inout PolySegment's
 328     std::vector< Path >&                  myPaths;    //!< path of each of segments to compute
 329     SMDS_Mesh*                            myMesh;
 330     mutable std::vector< std::string >    myErrors;
 331
 332     PolyPathCompute( SMESH_MeshAlgos::TListOfPolySegments& theSegments,
 333                      std::vector< Path >&                  thePaths,
 334                      SMDS_Mesh*                            theMesh):
 335       mySegments( theSegments ),
 336       myPaths( thePaths ),
 337       myMesh( theMesh ),
 338       myErrors( theSegments.size() )
 339     {
 340     }
 341
 342 #undef SMESH_CAUGHT
 343 #define SMESH_CAUGHT myErrors[i] =
 344     void operator() ( const int i ) const
 345     {
 346       SMESH_TRY;
 347       const_cast< PolyPathCompute* >( this )->Compute( i );
 348       SMESH_CATCH( SMESH::returnError );
 349     }
 350 #undef SMESH_CAUGHT
 351
 352     //================================================================================
 353     /*!
 354      * \brief Compute a path of a given segment
 355      */
 356     //================================================================================
 357
 358     void Compute( const int iSeg )
 359     {
 360       SMESH_MeshAlgos::PolySegment& polySeg = mySegments[ iSeg ];
 361
 362       if ( ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ).SquareModulus() == 0 )
 363       {
 364         myPaths[ iSeg ].AddPoint( polySeg.myXYZ[0] );
 365         myPaths[ iSeg ].AddPoint( polySeg.myXYZ[1] );
 366         return;
 367       }
 368
 369       // the cutting plane
 370       gp_XYZ plnNorm = ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ) ^ polySeg.myVector.XYZ();
 371       gp_XYZ plnOrig = polySeg.myXYZ[1];
 372
 373       // Find paths connecting the 2 end points of polySeg
 374
 375       std::vector< Path > paths; paths.reserve(10);
 376
 377       // 1) initialize paths; two paths starts at each end point
 378
 379       for ( int iP = 0; iP < 2; ++iP ) // loop on the polySeg end points
 380       {
 381         Path path( 0, iP );
 382         size_t nbPaths = paths.size();
 383
 384         if ( polySeg.myFace[ iP ]) // the end point lies on polySeg.myFace[ iP ]
 385         {
 386           // check coincidence of polySeg.myXYZ[ iP ] with nodes
 387           const double tol = 1e-17;
 388           SMESH_NodeXYZ nodes[4];
 389           for ( int i = 0; i < 3 &&  !polySeg.myNode1[ iP ]; ++i )
 390           {
 391             nodes[ i ] = polySeg.myFace[ iP ]->GetNode( i );
 392             if (( nodes[ i ] - polySeg.myXYZ[ iP ]).SquareModulus() < tol*tol )
 393               polySeg.myNode1[ iP ] = nodes[ i ].Node();
 394           }
 395           nodes[ 3 ] = nodes[ 0 ];
 396
 397           double dot[ 4 ];
 398           for ( int i = 0; i < 3; ++i )
 399             dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
 400           dot[ 3 ] = dot[ 0 ];
 401
 402           // check coincidence of polySeg.myXYZ[ iP ] with edges
 403           for ( int i = 0; i < 3 &&  !polySeg.myNode1[ iP ]; ++i )
 404           {
 405             SMDS_LinearEdge edge( nodes[i].Node(), nodes[i+1].Node() );
 406             if ( SMESH_MeshAlgos::GetDistance( &edge, polySeg.myXYZ[ iP ]) < tol )
 407             {
 408               polySeg.myNode1[ iP ] = nodes[ i     ].Node();
 409               polySeg.myNode2[ iP ] = nodes[ i + 1 ].Node();
 410
 411               int i3 = ( i + 2 ) % 3;
 412               if ( dot[ i   ] * dot [ i3 ] > 0 &&
 413                    dot[ i+1 ] * dot [ i3 ] > 0 ) // point iP is inside a neighbor triangle
 414               {
 415                 path.myAvoidSet.insert( polySeg.myFace[ iP ]);
 416                 const SMDS_MeshElement* face2 =
 417                   SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
 418                                                   polySeg.myNode2[ iP ],
 419                                                   path.myElemSet,
 420                                                   path.myAvoidSet );
 421                 if ( face2 )
 422                   polySeg.myFace[ iP ] = face2;
 423                 else
 424                   ;// ??
 425                 for ( int i = 0; i < 3; ++i )
 426                 {
 427                   nodes[ i ] = polySeg.myFace[ iP ]->GetNode( i );
 428                   dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
 429                 }
 430                 dot[ 3 ] = dot[ 0 ];
 431                 polySeg.myNode1[ iP ] = polySeg.myNode2[ iP ] = 0;
 432                 break;
 433               }
 434             }
 435           }
 436
 437           if ( !polySeg.myNode1[ iP ] ) // polySeg.myXYZ[ iP ] is within polySeg.myFace[ iP ]
 438           {
 439             int iCut = 0; // index of a cut edge
 440             if      ( dot[ 1 ] * dot[ 2 ] < 0. ) iCut = 1;
 441             else if ( dot[ 2 ] * dot[ 3 ] < 0. ) iCut = 2;
 442
 443             // initialize path so as if it entered the face via iCut-th edge
 444             path.myFace = polySeg.myFace[ iP ];
 445             path.myNodeInd1 = iCut;
 446             path.myNodeInd2 = iCut + 1;
 447             path.myNode1.Set( nodes[ iCut     ].Node() );
 448             path.myNode2.Set( nodes[ iCut + 1 ].Node() );
 449             path.myDot1 = dot[ iCut ];
 450             path.myDot2 = dot[ iCut + 1 ];
 451             path.myPoints.clear();
 452             path.AddPoint( polySeg.myXYZ[ iP ]);
 453             paths.push_back( path );
 454
 455             path.Extend( plnNorm, plnOrig ); // to get another edge cut
 456             path.myFace = polySeg.myFace[ iP ];
 457             if ( path.myDot1 == 0. ) // cut at a node
 458             {
 459               path.myNodeInd1 = ( iCut + 2 ) % 3;
 460               path.myNodeInd2 = ( iCut + 3 ) % 3;
 461               path.myNode2.Set( path.myFace->GetNode( path.myNodeInd2 ));
 462               path.myDot2 = dot[ path.myNodeInd2 ];
 463             }
 464             else
 465             {
 466               path.myNodeInd1 = path.myFace->GetNodeIndex( path.myNode1.Node() );
 467               path.myNodeInd2 = path.myFace->GetNodeIndex( path.myNode2.Node() );
 468             }
 469             path.myPoints.clear();
 470             path.AddPoint( polySeg.myXYZ[ iP ]);
 471             paths.push_back( path );
 472           }
 473         }
 474
 475         if ( polySeg.myNode2[ iP ] && polySeg.myNode2[ iP ] != polySeg.myNode1[ iP ] )
 476         {
 477           // the end point is on an edge
 478           while (( path.myFace = SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
 479                                                                  polySeg.myNode2[ iP ],
 480                                                                  path.myElemSet,
 481                                                                  path.myAvoidSet,
 482                                                                  &path.myNodeInd1,
 483                                                                  &path.myNodeInd2 )))
 484           {
 485             path.myNode1.Set( polySeg.myNode1[ iP ]);
 486             path.myNode2.Set( polySeg.myNode2[ iP ]);
 487             path.myDot1 = plnNorm * ( path.myNode1 - plnOrig );
 488             path.myDot2 = plnNorm * ( path.myNode2 - plnOrig );
 489             path.myPoints.clear();
 490             path.AddPoint( polySeg.myXYZ[ iP ]);
 491             path.myAvoidSet.insert( path.myFace );
 492             paths.push_back( path );
 493             std::swap( polySeg.myNode1[ iP ], polySeg.myNode2[ iP ]);
 494           }
 495           if ( nbPaths == paths.size() )
 496             throw SALOME_Exception ( SMESH_Comment("No face edge found by point ") << iP+1
 497                                      << " in a PolySegment " << iSeg );
 498
 499           if ( path.myDot1 == 0. &&
 500                path.myDot2 == 0. )
 501           {
 502             if ( paths.size() - nbPaths >= 2 ) // use a face non-parallel to the plane
 503             {
 504               const SMDS_MeshElement* goodFace = 0;
 505               for ( size_t j = nbPaths; j < paths.size(); ++j )
 506               {
 507                 path = paths[j];
 508                 if ( path.Extend( plnNorm, plnOrig ))
 509                   goodFace = paths[j].myFace;
 510                 else
 511                   paths[j].myFace = 0;
 512               }
 513               if ( !goodFace )
 514                 throw SALOME_Exception ( SMESH_Comment("Cant move from point ") << iP+1
 515                                          << " of a PolySegment " << iSeg );
 516               for ( size_t j = nbPaths; j < paths.size(); ++j )
 517                 if ( !paths[j].myFace )
 518                 {
 519                   paths[j].myFace = goodFace;
 520                   paths[j].myNodeInd1 = goodFace->GetNodeIndex( paths[j].myNode1.Node() );
 521                   paths[j].myNodeInd2 = goodFace->GetNodeIndex( paths[j].myNode2.Node() );
 522                 }
 523             }
 524             else // use the sole found face
 525             {
 526               path = paths.back();
 527               std::swap( path.myNode1,    path.myNode2 );
 528               std::swap( path.myNodeInd1, path.myNodeInd2 );
 529               paths.push_back( path );
 530             }
 531           }
 532         }
 533
 534         else if ( polySeg.myNode1[ iP ] ) // the end point is at a node
 535         {
 536           path.myFace = 0;
 537           path.SetCutAtCorner( polySeg.myNode1[ iP ], plnNorm, plnOrig, &paths );
 538         }
 539
 540
 541         // look for a one-segment path
 542         for ( size_t i = 0; i < nbPaths; ++i )
 543           for ( size_t j = nbPaths; j < paths.size(); ++j )
 544             if ( paths[i].myFace == paths[j].myFace )
 545             {
 546               myPaths[ iSeg ].myPoints.push_back( paths[i].myPoints[0] );
 547               myPaths[ iSeg ].myPoints.push_back( paths[j].myPoints[0] );
 548               paths.clear();
 549             }
 550
 551       }  // loop on the polySeg end points to initialize all possible paths
 552
 553
 554       // 2) extend paths and compose the shortest one connecting the two points
 555
 556       myPaths[ iSeg ].myLength = 1e100;
 557
 558       while ( paths.size() >= 2 )
 559       {
 560         for ( size_t i = 0; i < paths.size(); ++i )
 561         {
 562           Path& path = paths[ i ];
 563           if ( !path.Extend( plnNorm, plnOrig, &paths ) || // path reached a mesh boundary
 564                path.myLength > myPaths[ iSeg ].myLength )  // path is longer than others
 565           {
 566             Path::Remove( paths, i );
 567             continue;
 568           }
 569
 570           // join paths that reach same point
 571           for ( size_t j = 0; j < paths.size(); ++j )
 572           {
 573             if ( i != j && paths[i].ReachSamePoint( paths[j] ))
 574             {
 575               double   distLast = ( paths[i].myPoints.back() - paths[j].myPoints.back() ).Modulus();
 576               double fullLength = ( paths[i].myLength + paths[j].myLength + distLast );
 577               if ( fullLength < myPaths[ iSeg ].myLength )
 578               {
 579                 myPaths[ iSeg ].myLength = fullLength;
 580                 std::vector< gp_XYZ > & allPoints = myPaths[ iSeg ].myPoints;
 581                 allPoints.swap( paths[i].myPoints );
 582                 allPoints.insert( allPoints.end(),
 583                                   paths[j].myPoints.rbegin(),
 584                                   paths[j].myPoints.rend() );
 585               }
 586               if ( i < j ) std::swap( i, j );
 587               Path::Remove( paths, i );
 588               Path::Remove( paths, j );
 589               break;
 590             }
 591           }
 592         }
 593         if ( !paths.empty() && (int) paths[0].myPoints.size() > myMesh->NbFaces() )
 594           throw SALOME_Exception(LOCALIZED( "Infinite loop in MakePolyLine()"));
 595       }
 596
 597       if ( myPaths[ iSeg ].myPoints.empty() )
 598         throw SALOME_Exception( SMESH_Comment("Can't find a full path for PolySegment #") << iSeg );
 599
 600       // reverse the path
 601       double d00 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.front() ).SquareModulus();
 602       double d01 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.back()  ).SquareModulus();
 603       if ( d00 > d01 )
 604         std::reverse( myPaths[ iSeg ].myPoints.begin(), myPaths[ iSeg ].myPoints.end() );
 605
 606     } // PolyPathCompute::Compute()
 607
 608   }; // struct PolyPathCompute
 609
 610 } // namespace
 611
 612 //=======================================================================
 613 //function : MakePolyLine
 614 //purpose  : Create a polyline consisting of 1D mesh elements each lying on a 2D element of
 615 //           the initial mesh
 616 //=======================================================================
 617
 618 void SMESH_MeshAlgos::MakePolyLine( SMDS_Mesh*                            theMesh,
 619                                     TListOfPolySegments&                  theSegments,
 620                                     std::vector<const SMDS_MeshElement*>& theNewEdges,
 621                                     std::vector< const SMDS_MeshNode* >&  theNewNodes,
 622                                     SMDS_MeshGroup*                       theGroup,
 623                                     SMESH_ElementSearcher*                theSearcher)
 624 {
 625   std::vector< Path > segPaths( theSegments.size() ); // path of each of segments
 626
 627   SMESH_ElementSearcher* searcher = theSearcher;
 628   SMESHUtils::Deleter<SMESH_ElementSearcher> delSearcher;
 629   if ( !searcher )
 630   {
 631     searcher = SMESH_MeshAlgos::GetElementSearcher( *theMesh );
 632     delSearcher._obj = searcher;
 633   }
 634
 635   // get cutting planes
 636
 637   std::vector< bool > isVectorOK( theSegments.size(), true );
 638   const double planarCoef = 0.333; // plane height in planar case
 639
 640   for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
 641   {
 642     PolySegment& polySeg = theSegments[ iSeg ];
 643
 644     gp_XYZ p1 = SMESH_NodeXYZ( polySeg.myNode1[0] );
 645     gp_XYZ p2 = SMESH_NodeXYZ( polySeg.myNode1[1] );
 646     if ( polySeg.myNode2[0] ) p1 = 0.5 * ( p1 + SMESH_NodeXYZ( polySeg.myNode2[0] ));
 647     if ( polySeg.myNode2[1] ) p2 = 0.5 * ( p2 + SMESH_NodeXYZ( polySeg.myNode2[1] ));
 648
 649     polySeg.myFace[0] = polySeg.myFace[1] = 0;
 650     if ( !polySeg.myNode1[0] && !polySeg.myNode2[0] )
 651     {
 652       p1 = searcher->Project( polySeg.myXYZ[0], SMDSAbs_Face, &polySeg.myFace[0] );
 653     }
 654     if ( !polySeg.myNode1[1] && !polySeg.myNode2[1] )
 655     {
 656       p2 = searcher->Project( polySeg.myXYZ[1], SMDSAbs_Face, &polySeg.myFace[1] );
 657     }
 658     polySeg.myXYZ[0] = p1;
 659     polySeg.myXYZ[1] = p2;
 660
 661     gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
 662
 663     isVectorOK[ iSeg ] = ( plnNorm.Modulus() > std::numeric_limits<double>::min() );
 664     if ( !isVectorOK[ iSeg ] && ( p1 - p2 ).SquareModulus() > 0. )
 665     {
 666       gp_XYZ pMid = 0.5 * ( p1 + p2 );
 667       const SMDS_MeshElement* face;
 668       polySeg.myMidProjPoint = searcher->Project( pMid, SMDSAbs_Face, &face );
 669       polySeg.myVector       = polySeg.myMidProjPoint.XYZ() - pMid;
 670
 671       gp_XYZ faceNorm;
 672       SMESH_MeshAlgos::FaceNormal( face, faceNorm, /*normalized=*/false );
 673
 674       const double tol = Precision::Confusion();
 675       if ( polySeg.myVector.Magnitude() < tol || polySeg.myVector * faceNorm  < tol )
 676       {
 677         polySeg.myVector = faceNorm;
 678         polySeg.myMidProjPoint = pMid + faceNorm * ( p1 - p2 ).Modulus() * planarCoef;
 679       }
 680       plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
 681       if ( plnNorm.SquareModulus() == 0 ) // p1-p2 perpendicular to mesh
 682       {
 683         double radius = faceNorm.Modulus();
 684         std::vector< const SMDS_MeshElement* > foundElems;
 685         while ( plnNorm.SquareModulus() == 0  &&  radius < 1e200 )
 686         {
 687           foundElems.clear();
 688           searcher->GetElementsInSphere( p1, radius, SMDSAbs_Face, foundElems );
 689           searcher->GetElementsInSphere( p2, radius, SMDSAbs_Face, foundElems );
 690           radius *= 2;
 691           polySeg.myVector.SetCoord( 0,0,0 );
 692           for ( size_t i = 0; i < foundElems.size(); ++i )
 693           {
 694             SMESH_MeshAlgos::FaceNormal( foundElems[i], faceNorm );
 695             polySeg.myVector += faceNorm / foundElems.size();
 696           }
 697           plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
 698         }
 699       }
 700
 701     }
 702     else
 703     {
 704       polySeg.myVector = plnNorm ^ ( p1 - p2 );
 705     }
 706   }
 707
 708   // assure that inverse elements are constructed, avoid their concurrent building in threads
 709   theMesh->nodesIterator()->next()->NbInverseElements();
 710
 711   // find paths
 712
 713   PolyPathCompute algo( theSegments, segPaths, theMesh );
 714   OSD_Parallel::For( 0, theSegments.size(), algo, theSegments.size() == 1 );
 715
 716   for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
 717     if ( !algo.myErrors[ iSeg ].empty() )
 718       throw SALOME_Exception( algo.myErrors[ iSeg ].c_str() );
 719
 720   // create an 1D mesh
 721
 722   const SMDS_MeshNode *n, *nPrev = 0;
 723
 724   for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
 725   {
 726     const Path& path = segPaths[iSeg];
 727     if ( path.myPoints.size() < 2 )
 728       continue;
 729
 730     double tol = path.myLength / path.myPoints.size() / 1000.;
 731     if ( !nPrev || ( SMESH_NodeXYZ( nPrev ) - path.myPoints[0] ).SquareModulus() > tol*tol )
 732     {
 733       nPrev = theMesh->AddNode( path.myPoints[0].X(), path.myPoints[0].Y(), path.myPoints[0].Z() );
 734       theNewNodes.push_back( nPrev );
 735     }
 736     for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
 737     {
 738       n = theMesh->AddNode( path.myPoints[iP].X(), path.myPoints[iP].Y(), path.myPoints[iP].Z() );
 739       theNewNodes.push_back( n );
 740
 741       const SMDS_MeshElement* elem = theMesh->AddEdge( nPrev, n );
 742       theNewEdges.push_back( elem );
 743       if ( theGroup )
 744         theGroup->Add( elem );
 745
 746       nPrev = n;
 747     }
 748
 749     // return a vector
 750
 751     gp_XYZ pMid = 0.5 * ( path.myPoints[0] + path.myPoints.back() );
 752     if ( isVectorOK[ iSeg ])
 753     {
 754       // find the most distant point of a path
 755       double maxDist = 0;
 756       for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
 757       {
 758         double dist = Abs( theSegments[iSeg].myVector * ( path.myPoints[iP] - path.myPoints[0] ));
 759         if ( dist > maxDist )
 760         {
 761           maxDist = dist;
 762           theSegments[iSeg].myMidProjPoint = path.myPoints[iP];
 763         }
 764       }
 765       if ( maxDist < Precision::Confusion() ) // planar case
 766         theSegments[iSeg].myMidProjPoint =
 767           pMid + theSegments[iSeg].myVector.XYZ().Normalized() * path.myLength * planarCoef;
 768     }
 769     theSegments[iSeg].myVector = gp_Vec( pMid, theSegments[iSeg].myMidProjPoint );
 770   }
 771
 772   return;
 773 }