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