1 // Copyright (C) 2018-2020 CEA/DEN, EDF R&D, OPEN CASCADE
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // File : SMESH_PolyLine.cxx
20 // Created : Thu Dec 6 17:33:26 2018
21 // Author : Edward AGAPOV (eap)
23 #include "SMESH_MeshAlgos.hxx"
25 #include "SMDS_MeshGroup.hxx"
26 #include "SMDS_LinearEdge.hxx"
27 #include "SMDS_Mesh.hxx"
28 #include "SMESH_TryCatch.hxx"
30 #include <OSD_Parallel.hxx>
31 #include <Precision.hxx>
35 //================================================================================
37 * \brief Sequence of found points and a current point data
41 std::vector< gp_XYZ > myPoints;
44 const SMDS_MeshElement* myFace;
45 SMESH_NodeXYZ myNode1; // nodes of the edge the path entered myFace
46 SMESH_NodeXYZ myNode2;
52 int mySrcPntInd; //!< start point index
53 TIDSortedElemSet myElemSet, myAvoidSet;
55 Path(const SMDS_MeshElement* face=0, int srcInd=-1):
56 myLength(0.0), myFace(face), mySrcPntInd( srcInd ) {}
58 void CopyNoPoints( const Path& other );
60 bool Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig, std::vector< Path > * paths = 0 );
62 bool SetCutAtCorner( const SMESH_NodeXYZ& cornerNode,
63 const gp_XYZ& plnNorm,
64 const gp_XYZ& plnOrig,
65 std::vector< Path >* paths);
67 bool SetCutAtCorner( const SMESH_NodeXYZ& cornerNode,
68 const SMDS_MeshElement* face,
69 const gp_XYZ& plnNorm,
70 const gp_XYZ& plnOrig );
72 void AddPoint( const gp_XYZ& p );
74 bool ReachSamePoint( const Path& other );
76 static void Remove( std::vector< Path > & paths, size_t& i );
79 //================================================================================
81 * \brief Return true if this Path meats another
83 //================================================================================
85 bool Path::ReachSamePoint( const Path& other )
87 return ( mySrcPntInd != other.mySrcPntInd &&
88 myFace == other.myFace );
91 //================================================================================
93 * \brief Copy data except points
95 //================================================================================
97 void Path::CopyNoPoints( const Path& other )
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;
110 //================================================================================
112 * \brief Remove a path from a vector
114 //================================================================================
116 void Path::Remove( std::vector< Path > & paths, size_t& i )
118 if ( paths.size() > 1 )
120 size_t j = paths.size() - 1; // last item to be removed
123 paths[ i ].CopyNoPoints ( paths[ j ]);
124 paths[ i ].myPoints.swap( paths[ j ].myPoints );
132 //================================================================================
134 * \brief Try to extend self by a point located at a node.
135 * Return a success flag.
137 //================================================================================
139 bool Path::SetCutAtCorner( const SMESH_NodeXYZ& cornerNode,
140 const gp_XYZ& plnNorm,
141 const gp_XYZ& plnOrig,
142 std::vector< Path > * paths )
145 const bool isContinuation = myFace; // extend this path or find all possible paths?
146 const SMDS_MeshElement* lastFace = myFace;
150 SMDS_ElemIteratorPtr fIt = cornerNode->GetInverseElementIterator(SMDSAbs_Face);
151 while ( fIt->more() )
153 Path path( lastFace, mySrcPntInd );
154 if ( !path.SetCutAtCorner( cornerNode, fIt->next(), plnNorm, plnOrig ))
157 if ( path.myDot1 == 0 &&
158 !myAvoidSet.insert( path.myNode1.Node() ).second )
160 if ( path.myDot2 == 0 &&
161 !myAvoidSet.insert( path.myNode2.Node() ).second )
164 if ( isContinuation )
166 if ( ok ) // non-manifold continuation
168 path.myPoints = myPoints;
169 path.myLength = myLength;
170 path.AddPoint( cornerNode );
171 paths->push_back( path );
175 double len = myLength;
176 this->CopyNoPoints( path );
177 this->myLength = len;
178 this->AddPoint( path.myPoints.back() );
183 paths->push_back( path );
191 //================================================================================
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
196 //================================================================================
198 bool Path::SetCutAtCorner( const SMESH_NodeXYZ& cornerNode,
199 const SMDS_MeshElement* face,
200 const gp_XYZ& plnNorm,
201 const gp_XYZ& plnOrig )
203 if ( face == myFace )
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 ));
211 myDot1 = plnNorm * ( myNode1 - plnOrig );
212 myDot2 = plnNorm * ( myNode2 - plnOrig );
214 bool ok = ( myDot1 * myDot2 < 0 );
215 if ( !ok && myDot1 * myDot2 == 0 )
217 ok = ( myDot1 != myDot2 );
219 ok = ( myFace->GetNodeIndex(( myDot1 == 0 ? myNode1 : myNode2 )._node ) < 0 );
225 AddPoint( cornerNode );
230 //================================================================================
232 * \brief Store a point and update myLength
234 //================================================================================
236 void Path::AddPoint( const gp_XYZ& p )
238 if ( !myPoints.empty() )
239 myLength += ( p - myPoints.back() ).Modulus();
242 myPoints.push_back( p );
245 //================================================================================
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
252 //================================================================================
254 bool Path::Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig, std::vector< Path > * paths )
258 int nodeInd3 = ( myNodeInd1 + 1 ) % myFace->NbCornerNodes();
259 if ( myNodeInd2 == nodeInd3 )
260 nodeInd3 = ( myNodeInd1 + 2 ) % myFace->NbCornerNodes();
262 SMESH_NodeXYZ node3 = myFace->GetNode( nodeInd3 );
263 double dot3 = plnNorm * ( node3 - plnOrig );
265 if ( dot3 * myDot1 < 0. )
268 myNodeInd2 = nodeInd3;
271 else if ( dot3 * myDot2 < 0. )
274 myNodeInd1 = nodeInd3;
277 else if ( dot3 == 0. )
279 ok = SetCutAtCorner( node3, plnNorm, plnOrig, paths );
282 else if ( myDot2 == 0. )
284 ok = SetCutAtCorner( myNode2, plnNorm, plnOrig, paths );
288 double r = Abs( myDot1 / ( myDot2 - myDot1 ));
289 AddPoint( myNode1 * ( 1 - r ) + myNode2 * r );
292 myAvoidSet.insert( myFace );
293 const SMDS_MeshElement* nextFace;
295 while (( nextFace = SMESH_MeshAlgos::FindFaceInSet( myNode1._node, myNode2._node,
296 myElemSet, myAvoidSet,
299 if ( ok ) // non-manifold continuation
301 paths->push_back( *this );
302 paths->back().myFace = nextFace;
303 paths->back().myNodeInd1 = ind1;
304 paths->back().myNodeInd2 = ind2;
315 myAvoidSet.insert( nextFace );
321 //================================================================================
323 * \brief Compute a path between two points of PolySegment
325 struct PolyPathCompute
327 SMESH_MeshAlgos::TListOfPolySegments& mySegments; //!< inout PolySegment's
328 std::vector< Path >& myPaths; //!< path of each of segments to compute
330 mutable std::vector< std::string > myErrors;
332 PolyPathCompute( SMESH_MeshAlgos::TListOfPolySegments& theSegments,
333 std::vector< Path >& thePaths,
335 mySegments( theSegments ),
338 myErrors( theSegments.size() )
343 #define SMESH_CAUGHT myErrors[i] =
344 void operator() ( const int i ) const
347 const_cast< PolyPathCompute* >( this )->Compute( i );
348 SMESH_CATCH( SMESH::returnError );
352 //================================================================================
354 * \brief Compute a path of a given segment
356 //================================================================================
358 void Compute( const int iSeg )
360 SMESH_MeshAlgos::PolySegment& polySeg = mySegments[ iSeg ];
362 if (( polySeg.myXYZ[0] - polySeg.myXYZ[1] ).SquareModulus() == 0 )
364 myPaths[ iSeg ].AddPoint( polySeg.myXYZ[0] );
365 myPaths[ iSeg ].AddPoint( polySeg.myXYZ[1] );
370 gp_XYZ plnNorm = ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ) ^ polySeg.myVector.XYZ();
371 gp_XYZ plnOrig = polySeg.myXYZ[1];
373 // Find paths connecting the 2 end points of polySeg
375 std::vector< Path > paths; paths.reserve(10);
377 // 1) initialize paths; two paths starts at each end point
379 for ( int iP = 0; iP < 2; ++iP ) // loop on the polySeg end points
382 size_t nbPaths = paths.size();
384 if ( polySeg.myFace[ iP ]) // the end point lies on polySeg.myFace[ iP ]
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 )
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();
395 nodes[ 3 ] = nodes[ 0 ];
398 for ( int i = 0; i < 3; ++i )
399 dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
402 // check coincidence of polySeg.myXYZ[ iP ] with edges
403 for ( int i = 0; i < 3 && !polySeg.myNode1[ iP ]; ++i )
405 SMDS_LinearEdge edge( nodes[i].Node(), nodes[i+1].Node() );
406 if ( SMESH_MeshAlgos::GetDistance( &edge, polySeg.myXYZ[ iP ]) < tol )
408 polySeg.myNode1[ iP ] = nodes[ i ].Node();
409 polySeg.myNode2[ iP ] = nodes[ i + 1 ].Node();
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
415 path.myAvoidSet.insert( polySeg.myFace[ iP ]);
416 const SMDS_MeshElement* face2 =
417 SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
418 polySeg.myNode2[ iP ],
422 polySeg.myFace[ iP ] = face2;
425 for ( int i = 0; i < 3; ++i )
427 nodes[ i ] = polySeg.myFace[ iP ]->GetNode( i );
428 dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
431 polySeg.myNode1[ iP ] = polySeg.myNode2[ iP ] = 0;
437 if ( !polySeg.myNode1[ iP ] ) // polySeg.myXYZ[ iP ] is within polySeg.myFace[ iP ]
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;
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 );
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
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 ];
466 path.myNodeInd1 = path.myFace->GetNodeIndex( path.myNode1.Node() );
467 path.myNodeInd2 = path.myFace->GetNodeIndex( path.myNode2.Node() );
469 path.myPoints.clear();
470 path.AddPoint( polySeg.myXYZ[ iP ]);
471 paths.push_back( path );
475 if ( polySeg.myNode2[ iP ] && polySeg.myNode2[ iP ] != polySeg.myNode1[ iP ] )
477 // the end point is on an edge
478 while (( path.myFace = SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
479 polySeg.myNode2[ iP ],
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 ]);
495 if ( nbPaths == paths.size() )
496 throw SALOME_Exception ( SMESH_Comment("No face edge found by point ") << iP+1
497 << " in a PolySegment " << iSeg );
499 if ( path.myDot1 == 0. &&
502 if ( paths.size() - nbPaths >= 2 ) // use a face non-parallel to the plane
504 const SMDS_MeshElement* goodFace = 0;
505 for ( size_t j = nbPaths; j < paths.size(); ++j )
508 if ( path.Extend( plnNorm, plnOrig ))
509 goodFace = paths[j].myFace;
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 )
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() );
524 else // use the sole found face
527 std::swap( path.myNode1, path.myNode2 );
528 std::swap( path.myNodeInd1, path.myNodeInd2 );
529 paths.push_back( path );
534 else if ( polySeg.myNode1[ iP ] ) // the end point is at a node
537 path.SetCutAtCorner( polySeg.myNode1[ iP ], plnNorm, plnOrig, &paths );
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 )
546 myPaths[ iSeg ].myPoints.push_back( paths[i].myPoints[0] );
547 myPaths[ iSeg ].myPoints.push_back( paths[j].myPoints[0] );
551 } // loop on the polySeg end points to initialize all possible paths
554 // 2) extend paths and compose the shortest one connecting the two points
556 myPaths[ iSeg ].myLength = 1e100;
558 while ( paths.size() >= 2 )
560 for ( size_t i = 0; i < paths.size(); ++i )
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
566 Path::Remove( paths, i );
570 // join paths that reach same point
571 for ( size_t j = 0; j < paths.size(); ++j )
573 if ( i != j && paths[i].ReachSamePoint( paths[j] ))
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 )
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() );
586 if ( i < j ) std::swap( i, j );
587 Path::Remove( paths, i );
588 Path::Remove( paths, j );
593 if ( !paths.empty() && (int) paths[0].myPoints.size() > myMesh->NbFaces() )
594 throw SALOME_Exception(LOCALIZED( "Infinite loop in MakePolyLine()"));
597 if ( myPaths[ iSeg ].myPoints.empty() )
598 throw SALOME_Exception( SMESH_Comment("Can't find a full path for PolySegment #") << iSeg );
601 double d00 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.front() ).SquareModulus();
602 double d01 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.back() ).SquareModulus();
604 std::reverse( myPaths[ iSeg ].myPoints.begin(), myPaths[ iSeg ].myPoints.end() );
606 } // PolyPathCompute::Compute()
608 }; // struct PolyPathCompute
612 //=======================================================================
613 //function : MakePolyLine
614 //purpose : Create a polyline consisting of 1D mesh elements each lying on a 2D element of
616 //=======================================================================
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)
625 std::vector< Path > segPaths( theSegments.size() ); // path of each of segments
627 SMESH_ElementSearcher* searcher = theSearcher;
628 SMESHUtils::Deleter<SMESH_ElementSearcher> delSearcher;
631 searcher = SMESH_MeshAlgos::GetElementSearcher( *theMesh );
632 delSearcher._obj = searcher;
635 // get cutting planes
637 std::vector< bool > isVectorOK( theSegments.size(), true );
638 const double planarCoef = 0.333; // plane height in planar case
640 for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
642 PolySegment& polySeg = theSegments[ iSeg ];
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] ));
649 polySeg.myFace[0] = polySeg.myFace[1] = 0;
650 if ( !polySeg.myNode1[0] && !polySeg.myNode2[0] )
652 p1 = searcher->Project( polySeg.myXYZ[0], SMDSAbs_Face, &polySeg.myFace[0] );
654 if ( !polySeg.myNode1[1] && !polySeg.myNode2[1] )
656 p2 = searcher->Project( polySeg.myXYZ[1], SMDSAbs_Face, &polySeg.myFace[1] );
658 polySeg.myXYZ[0] = p1;
659 polySeg.myXYZ[1] = p2;
661 gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
663 isVectorOK[ iSeg ] = ( plnNorm.Modulus() > std::numeric_limits<double>::min() );
664 if ( !isVectorOK[ iSeg ] && ( p1 - p2 ).SquareModulus() > 0. )
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;
672 SMESH_MeshAlgos::FaceNormal( face, faceNorm, /*normalized=*/false );
674 const double tol = Precision::Confusion();
675 if ( polySeg.myVector.Magnitude() < tol || polySeg.myVector * faceNorm < tol )
677 polySeg.myVector = faceNorm;
678 polySeg.myMidProjPoint = pMid + faceNorm * ( p1 - p2 ).Modulus() * planarCoef;
680 plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
681 if ( plnNorm.SquareModulus() == 0 ) // p1-p2 perpendicular to mesh
683 double radius = faceNorm.Modulus();
684 std::vector< const SMDS_MeshElement* > foundElems;
685 while ( plnNorm.SquareModulus() == 0 && radius < 1e200 )
688 searcher->GetElementsInSphere( p1, radius, SMDSAbs_Face, foundElems );
689 searcher->GetElementsInSphere( p2, radius, SMDSAbs_Face, foundElems );
691 polySeg.myVector.SetCoord( 0,0,0 );
692 for ( size_t i = 0; i < foundElems.size(); ++i )
694 SMESH_MeshAlgos::FaceNormal( foundElems[i], faceNorm );
695 polySeg.myVector += faceNorm / foundElems.size();
697 plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
704 polySeg.myVector = plnNorm ^ ( p1 - p2 );
708 // assure that inverse elements are constructed, avoid their concurrent building in threads
709 theMesh->nodesIterator()->next()->NbInverseElements();
713 PolyPathCompute algo( theSegments, segPaths, theMesh );
714 OSD_Parallel::For( 0, theSegments.size(), algo, theSegments.size() == 1 );
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() );
722 const SMDS_MeshNode *n, *nPrev = 0;
724 for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
726 const Path& path = segPaths[iSeg];
727 if ( path.myPoints.size() < 2 )
730 double tol = path.myLength / path.myPoints.size() / 1000.;
731 if ( !nPrev || ( SMESH_NodeXYZ( nPrev ) - path.myPoints[0] ).SquareModulus() > tol*tol )
733 nPrev = theMesh->AddNode( path.myPoints[0].X(), path.myPoints[0].Y(), path.myPoints[0].Z() );
734 theNewNodes.push_back( nPrev );
736 for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
738 n = theMesh->AddNode( path.myPoints[iP].X(), path.myPoints[iP].Y(), path.myPoints[iP].Z() );
739 theNewNodes.push_back( n );
741 const SMDS_MeshElement* elem = theMesh->AddEdge( nPrev, n );
742 theNewEdges.push_back( elem );
744 theGroup->Add( elem );
751 gp_XYZ pMid = 0.5 * ( path.myPoints[0] + path.myPoints.back() );
752 if ( isVectorOK[ iSeg ])
754 // find the most distant point of a path
756 for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
758 double dist = Abs( theSegments[iSeg].myVector * ( path.myPoints[iP] - path.myPoints[0] ));
759 if ( dist > maxDist )
762 theSegments[iSeg].myMidProjPoint = path.myPoints[iP];
765 if ( maxDist < Precision::Confusion() ) // planar case
766 theSegments[iSeg].myMidProjPoint =
767 pMid + theSegments[iSeg].myVector.XYZ().Normalized() * path.myLength * planarCoef;
769 theSegments[iSeg].myVector = gp_Vec( pMid, theSegments[iSeg].myMidProjPoint );