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           {
 499             if ( paths.size() - nbPaths >= 2 ) // use a face non-parallel to the plane
 500             {
 501               const SMDS_MeshElement* goodFace = 0;
 502               for ( size_t j = nbPaths; j < paths.size(); ++j )
 503               {
 504                 path = paths[j];
 505                 if ( path.Extend( plnNorm, plnOrig ))
 506                   goodFace = paths[j].myFace;
 507                 else
 508                   paths[j].myFace = 0;
 509               }
 510               if ( !goodFace )
 511                 throw SALOME_Exception ( SMESH_Comment("Cant move from point ") << iP+1
 512                                          << " of a PolySegment " << iSeg );
 513               for ( size_t j = nbPaths; j < paths.size(); ++j )
 514                 if ( !paths[j].myFace )
 515                 {
 516                   paths[j].myFace = goodFace;
 517                   paths[j].myNodeInd1 = goodFace->GetNodeIndex( paths[j].myNode1.Node() );
 518                   paths[j].myNodeInd2 = goodFace->GetNodeIndex( paths[j].myNode2.Node() );
 519                 }
 520             }
 521             else // use the sole found face
 522             {
 523               path = paths.back();
 524               std::swap( path.myNode1,    path.myNode2 );
 525               std::swap( path.myNodeInd1, path.myNodeInd2 );
 526               paths.push_back( path );
 527             }
 528           }
 529         }
 530
 531         else if ( polySeg.myNode1[ iP ] ) // the end point is at a node
 532         {
 533           path.myFace = 0;
 534           path.SetCutAtCorner( polySeg.myNode1[ iP ], plnNorm, plnOrig, &paths );
 535         }
 536
 537
 538         // look for a one-segment path
 539         for ( size_t i = 0; i < nbPaths; ++i )
 540           for ( size_t j = nbPaths; j < paths.size(); ++j )
 541             if ( paths[i].myFace == paths[j].myFace )
 542             {
 543               myPaths[ iSeg ].myPoints.push_back( paths[i].myPoints[0] );
 544               myPaths[ iSeg ].myPoints.push_back( paths[j].myPoints[0] );
 545               paths.clear();
 546             }
 547
 548       }  // loop on the polySeg end points to initialize all possible paths
 549
 550
 551       // 2) extend paths and compose the shortest one connecting the two points
 552
 553       myPaths[ iSeg ].myLength = 1e100;
 554
 555       while ( paths.size() >= 2 )
 556       {
 557         for ( size_t i = 0; i < paths.size(); ++i )
 558         {
 559           Path& path = paths[ i ];
 560           if ( !path.Extend( plnNorm, plnOrig, &paths ) || // path reached a mesh boundary
 561                path.myLength > myPaths[ iSeg ].myLength )  // path is longer than others
 562           {
 563             Path::Remove( paths, i );
 564             continue;
 565           }
 566
 567           // join paths that reach same point
 568           for ( size_t j = 0; j < paths.size(); ++j )
 569           {
 570             if ( i != j && paths[i].ReachSamePoint( paths[j] ))
 571             {
 572               double   distLast = ( paths[i].myPoints.back() - paths[j].myPoints.back() ).Modulus();
 573               double fullLength = ( paths[i].myLength + paths[j].myLength + distLast );
 574               if ( fullLength < myPaths[ iSeg ].myLength )
 575               {
 576                 myPaths[ iSeg ].myLength = fullLength;
 577                 std::vector< gp_XYZ > & allPoints = myPaths[ iSeg ].myPoints;
 578                 allPoints.swap( paths[i].myPoints );
 579                 allPoints.insert( allPoints.end(),
 580                                   paths[j].myPoints.rbegin(),
 581                                   paths[j].myPoints.rend() );
 582               }
 583               if ( i < j ) std::swap( i, j );
 584               Path::Remove( paths, i );
 585               Path::Remove( paths, j );
 586               break;
 587             }
 588           }
 589         }
 590         if ( !paths.empty() && (int) paths[0].myPoints.size() > myMesh->NbFaces() )
 591           throw SALOME_Exception(LOCALIZED( "Infinite loop in MakePolyLine()"));
 592       }
 593
 594       if ( myPaths[ iSeg ].myPoints.empty() )
 595         throw SALOME_Exception( SMESH_Comment("Can't find a full path for PolySegment #") << iSeg );
 596
 597       // reverse the path
 598       double d00 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.front() ).SquareModulus();
 599       double d01 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.back()  ).SquareModulus();
 600       if ( d00 > d01 )
 601         std::reverse( myPaths[ iSeg ].myPoints.begin(), myPaths[ iSeg ].myPoints.end() );
 602
 603     } // PolyPathCompute::Compute()
 604
 605   }; // struct PolyPathCompute
 606
 607 } // namespace
 608
 609 //=======================================================================
 610 //function : MakePolyLine
 611 //purpose  : Create a polyline consisting of 1D mesh elements each lying on a 2D element of
 612 //           the initial mesh
 613 //=======================================================================
 614
 615 void SMESH_MeshAlgos::MakePolyLine( SMDS_Mesh*                            theMesh,
 616                                     TListOfPolySegments&                  theSegments,
 617                                     std::vector<const SMDS_MeshElement*>& theNewEdges,
 618                                     std::vector< const SMDS_MeshNode* >&  theNewNodes,
 619                                     SMDS_MeshGroup*                       theGroup,
 620                                     SMESH_ElementSearcher*                theSearcher)
 621 {
 622   std::vector< Path > segPaths( theSegments.size() ); // path of each of segments
 623
 624   SMESH_ElementSearcher* searcher = theSearcher;
 625   SMESHUtils::Deleter<SMESH_ElementSearcher> delSearcher;
 626   if ( !searcher )
 627   {
 628     searcher = SMESH_MeshAlgos::GetElementSearcher( *theMesh );
 629     delSearcher._obj = searcher;
 630   }
 631
 632   // get cutting planes
 633
 634   std::vector< bool > isVectorOK( theSegments.size(), true );
 635   const double planarCoef = 0.333; // plane height in planar case
 636
 637   for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
 638   {
 639     PolySegment& polySeg = theSegments[ iSeg ];
 640
 641     gp_XYZ p1 = SMESH_NodeXYZ( polySeg.myNode1[0] );
 642     gp_XYZ p2 = SMESH_NodeXYZ( polySeg.myNode1[1] );
 643     if ( polySeg.myNode2[0] ) p1 = 0.5 * ( p1 + SMESH_NodeXYZ( polySeg.myNode2[0] ));
 644     if ( polySeg.myNode2[1] ) p2 = 0.5 * ( p2 + SMESH_NodeXYZ( polySeg.myNode2[1] ));
 645
 646     polySeg.myFace[0] = polySeg.myFace[1] = 0;
 647     if ( !polySeg.myNode1[0] && !polySeg.myNode2[0] )
 648     {
 649       p1 = searcher->Project( polySeg.myXYZ[0], SMDSAbs_Face, &polySeg.myFace[0] );
 650     }
 651     if ( !polySeg.myNode1[1] && !polySeg.myNode2[1] )
 652     {
 653       p2 = searcher->Project( polySeg.myXYZ[1], SMDSAbs_Face, &polySeg.myFace[1] );
 654     }
 655     polySeg.myXYZ[0] = p1;
 656     polySeg.myXYZ[1] = p2;
 657
 658     gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
 659
 660     isVectorOK[ iSeg ] = ( plnNorm.Modulus() > std::numeric_limits<double>::min() );
 661     if ( !isVectorOK[ iSeg ] && ( p1 - p2 ).SquareModulus() > 0. )
 662     {
 663       gp_XYZ pMid = 0.5 * ( p1 + p2 );
 664       const SMDS_MeshElement* face;
 665       polySeg.myMidProjPoint = searcher->Project( pMid, SMDSAbs_Face, &face );
 666       polySeg.myVector       = polySeg.myMidProjPoint.XYZ() - pMid;
 667
 668       gp_XYZ faceNorm;
 669       SMESH_MeshAlgos::FaceNormal( face, faceNorm, /*normalized=*/false );
 670
 671       const double tol = Precision::Confusion();
 672       if ( polySeg.myVector.Magnitude() < tol || polySeg.myVector * faceNorm  < tol )
 673       {
 674         polySeg.myVector = faceNorm;
 675         polySeg.myMidProjPoint = pMid + faceNorm * ( p1 - p2 ).Modulus() * planarCoef;
 676       }
 677       plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
 678       if ( plnNorm.SquareModulus() == 0 ) // p1-p2 perpendicular to mesh
 679       {
 680         double radius = faceNorm.Modulus();
 681         std::vector< const SMDS_MeshElement* > foundElems;
 682         while ( plnNorm.SquareModulus() == 0  &&  radius < 1e200 )
 683         {
 684           foundElems.clear();
 685           searcher->GetElementsInSphere( p1, radius, SMDSAbs_Face, foundElems );
 686           searcher->GetElementsInSphere( p2, radius, SMDSAbs_Face, foundElems );
 687           radius *= 2;
 688           polySeg.myVector.SetCoord( 0,0,0 );
 689           for ( size_t i = 0; i < foundElems.size(); ++i )
 690           {
 691             SMESH_MeshAlgos::FaceNormal( foundElems[i], faceNorm );
 692             polySeg.myVector += faceNorm / foundElems.size();
 693           }
 694           plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
 695         }
 696       }
 697
 698     }
 699     else
 700     {
 701       polySeg.myVector = plnNorm ^ ( p1 - p2 );
 702     }
 703   }
 704
 705   // assure that inverse elements are constructed, avoid their concurrent building in threads
 706   theMesh->nodesIterator()->next()->NbInverseElements();
 707
 708   // find paths
 709
 710   PolyPathCompute algo( theSegments, segPaths, theMesh );
 711   OSD_Parallel::For( 0, theSegments.size(), algo, theSegments.size() == 1 );
 712
 713   for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
 714     if ( !algo.myErrors[ iSeg ].empty() )
 715       throw SALOME_Exception( algo.myErrors[ iSeg ].c_str() );
 716
 717   // create an 1D mesh
 718
 719   const SMDS_MeshNode *n, *nPrev = 0;
 720
 721   for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
 722   {
 723     const Path& path = segPaths[iSeg];
 724     if ( path.myPoints.size() < 2 )
 725       continue;
 726
 727     double tol = path.myLength / path.myPoints.size() / 1000.;
 728     if ( !nPrev || ( SMESH_NodeXYZ( nPrev ) - path.myPoints[0] ).SquareModulus() > tol*tol )
 729     {
 730       nPrev = theMesh->AddNode( path.myPoints[0].X(), path.myPoints[0].Y(), path.myPoints[0].Z() );
 731       theNewNodes.push_back( nPrev );
 732     }
 733     for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
 734     {
 735       n = theMesh->AddNode( path.myPoints[iP].X(), path.myPoints[iP].Y(), path.myPoints[iP].Z() );
 736       theNewNodes.push_back( n );
 737
 738       const SMDS_MeshElement* elem = theMesh->AddEdge( nPrev, n );
 739       theNewEdges.push_back( elem );
 740       if ( theGroup )
 741         theGroup->Add( elem );
 742
 743       nPrev = n;
 744     }
 745
 746     // return a vector
 747
 748     gp_XYZ pMid = 0.5 * ( path.myPoints[0] + path.myPoints.back() );
 749     if ( isVectorOK[ iSeg ])
 750     {
 751       // find the most distant point of a path
 752       double maxDist = 0;
 753       for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
 754       {
 755         double dist = Abs( theSegments[iSeg].myVector * ( path.myPoints[iP] - path.myPoints[0] ));
 756         if ( dist > maxDist )
 757         {
 758           maxDist = dist;
 759           theSegments[iSeg].myMidProjPoint = path.myPoints[iP];
 760         }
 761       }
 762       if ( maxDist < Precision::Confusion() ) // planar case
 763         theSegments[iSeg].myMidProjPoint =
 764           pMid + theSegments[iSeg].myVector.XYZ().Normalized() * path.myLength * planarCoef;
 765     }
 766     theSegments[iSeg].myVector = gp_Vec( pMid, theSegments[iSeg].myMidProjPoint );
 767   }
 768
 769   return;
 770 }