1 // Copyright (C) 2018-2019 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 ( !myAvoidSet.insert( path.myNode1.Node() ).second ||
158 !myAvoidSet.insert( path.myNode2.Node() ).second )
161 if ( isContinuation )
163 if ( ok ) // non-manifold continuation
165 path.myPoints = myPoints;
166 path.myLength = myLength;
167 path.AddPoint( cornerNode );
168 paths->push_back( path );
172 double len = myLength;
173 this->CopyNoPoints( path );
174 this->myLength = len;
175 this->AddPoint( path.myPoints.back() );
180 paths->push_back( path );
188 //================================================================================
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
193 //================================================================================
195 bool Path::SetCutAtCorner( const SMESH_NodeXYZ& cornerNode,
196 const SMDS_MeshElement* face,
197 const gp_XYZ& plnNorm,
198 const gp_XYZ& plnOrig )
200 if ( face == myFace )
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 ));
208 myDot1 = plnNorm * ( myNode1 - plnOrig );
209 myDot2 = plnNorm * ( myNode2 - plnOrig );
211 bool ok = ( myDot1 * myDot2 < 0 );
212 if ( !ok && myDot1 * myDot2 == 0 )
214 ok = ( myDot1 != myDot2 );
216 ok = ( myFace->GetNodeIndex(( myDot1 == 0 ? myNode1 : myNode2 )._node ) < 0 );
222 AddPoint( cornerNode );
227 //================================================================================
229 * \brief Store a point and update myLength
231 //================================================================================
233 void Path::AddPoint( const gp_XYZ& p )
235 if ( !myPoints.empty() )
236 myLength += ( p - myPoints.back() ).Modulus();
239 myPoints.push_back( p );
242 //================================================================================
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
249 //================================================================================
251 bool Path::Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig, std::vector< Path > * paths )
255 int nodeInd3 = ( myNodeInd1 + 1 ) % myFace->NbCornerNodes();
256 if ( myNodeInd2 == nodeInd3 )
257 nodeInd3 = ( myNodeInd1 + 2 ) % myFace->NbCornerNodes();
259 SMESH_NodeXYZ node3 = myFace->GetNode( nodeInd3 );
260 double dot3 = plnNorm * ( node3 - plnOrig );
262 if ( dot3 * myDot1 < 0. )
265 myNodeInd2 = nodeInd3;
268 else if ( dot3 * myDot2 < 0. )
271 myNodeInd1 = nodeInd3;
274 else if ( dot3 == 0. )
276 ok = SetCutAtCorner( node3, plnNorm, plnOrig, paths );
279 else if ( myDot2 == 0. )
281 ok = SetCutAtCorner( myNode2, plnNorm, plnOrig, paths );
285 double r = Abs( myDot1 / ( myDot2 - myDot1 ));
286 AddPoint( myNode1 * ( 1 - r ) + myNode2 * r );
289 myAvoidSet.insert( myFace );
290 const SMDS_MeshElement* nextFace;
292 while (( nextFace = SMESH_MeshAlgos::FindFaceInSet( myNode1._node, myNode2._node,
293 myElemSet, myAvoidSet,
296 if ( ok ) // non-manifold continuation
298 paths->push_back( *this );
299 paths->back().myFace = nextFace;
300 paths->back().myNodeInd1 = ind1;
301 paths->back().myNodeInd2 = ind2;
312 myAvoidSet.insert( nextFace );
318 //================================================================================
320 * \brief Compute a path between two points of PolySegment
322 struct PolyPathCompute
324 SMESH_MeshAlgos::TListOfPolySegments& mySegments; //!< inout PolySegment's
325 std::vector< Path >& myPaths; //!< path of each of segments to compute
327 mutable std::vector< std::string > myErrors;
329 PolyPathCompute( SMESH_MeshAlgos::TListOfPolySegments& theSegments,
330 std::vector< Path >& thePaths,
332 mySegments( theSegments ),
335 myErrors( theSegments.size() )
340 #define SMESH_CAUGHT myErrors[i] =
341 void operator() ( const int i ) const
344 const_cast< PolyPathCompute* >( this )->Compute( i );
345 SMESH_CATCH( SMESH::returnError );
349 //================================================================================
351 * \brief Compute a path of a given segment
353 //================================================================================
355 void Compute( const int iSeg )
357 SMESH_MeshAlgos::PolySegment& polySeg = mySegments[ iSeg ];
359 if ( ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ).SquareModulus() == 0 )
361 myPaths[ iSeg ].AddPoint( polySeg.myXYZ[0] );
362 myPaths[ iSeg ].AddPoint( polySeg.myXYZ[1] );
367 gp_XYZ plnNorm = ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ) ^ polySeg.myVector.XYZ();
368 gp_XYZ plnOrig = polySeg.myXYZ[1];
370 // Find paths connecting the 2 end points of polySeg
372 std::vector< Path > paths; paths.reserve(10);
374 // 1) initialize paths; two paths starts at each end point
376 for ( int iP = 0; iP < 2; ++iP ) // loop on the polySeg end points
379 size_t nbPaths = paths.size();
381 if ( polySeg.myFace[ iP ]) // the end point lies on polySeg.myFace[ iP ]
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 )
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();
392 nodes[ 3 ] = nodes[ 0 ];
395 for ( int i = 0; i < 3; ++i )
396 dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
399 // check coincidence of polySeg.myXYZ[ iP ] with edges
400 for ( int i = 0; i < 3 && !polySeg.myNode1[ iP ]; ++i )
402 SMDS_LinearEdge edge( nodes[i].Node(), nodes[i+1].Node() );
403 if ( SMESH_MeshAlgos::GetDistance( &edge, polySeg.myXYZ[ iP ]) < tol )
405 polySeg.myNode1[ iP ] = nodes[ i ].Node();
406 polySeg.myNode2[ iP ] = nodes[ i + 1 ].Node();
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
412 path.myAvoidSet.insert( polySeg.myFace[ iP ]);
413 const SMDS_MeshElement* face2 =
414 SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
415 polySeg.myNode2[ iP ],
419 polySeg.myFace[ iP ] = face2;
422 for ( int i = 0; i < 3; ++i )
424 nodes[ i ] = polySeg.myFace[ iP ]->GetNode( i );
425 dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
428 polySeg.myNode1[ iP ] = polySeg.myNode2[ iP ] = 0;
434 if ( !polySeg.myNode1[ iP ] ) // polySeg.myXYZ[ iP ] is within polySeg.myFace[ iP ]
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;
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 );
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
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 ];
463 path.myNodeInd1 = path.myFace->GetNodeIndex( path.myNode1.Node() );
464 path.myNodeInd2 = path.myFace->GetNodeIndex( path.myNode2.Node() );
466 path.myPoints.clear();
467 path.AddPoint( polySeg.myXYZ[ iP ]);
468 paths.push_back( path );
472 if ( polySeg.myNode2[ iP ] && polySeg.myNode2[ iP ] != polySeg.myNode1[ iP ] )
474 // the end point is on an edge
475 while (( path.myFace = SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
476 polySeg.myNode2[ iP ],
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 ]);
492 if ( nbPaths == paths.size() )
493 throw SALOME_Exception ( SMESH_Comment("No face edge found by point ") << iP+1
494 << " in a PolySegment " << iSeg );
496 if ( path.myDot1 == 0. &&
498 paths.size() - nbPaths >= 2 ) // use a face non-parallel to the plane
500 const SMDS_MeshElement* goodFace = 0;
501 for ( size_t j = nbPaths; j < paths.size(); ++j )
504 if ( path.Extend( plnNorm, plnOrig ))
505 goodFace = paths[j].myFace;
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 )
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() );
522 else if ( polySeg.myNode1[ iP ] ) // the end point is at a node
525 path.SetCutAtCorner( polySeg.myNode1[ iP ], plnNorm, plnOrig, &paths );
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 )
534 myPaths[ iSeg ].myPoints.push_back( paths[i].myPoints[0] );
535 myPaths[ iSeg ].myPoints.push_back( paths[j].myPoints[0] );
539 } // loop on the polySeg end points to initialize all possible paths
542 // 2) extend paths and compose the shortest one connecting the two points
544 myPaths[ iSeg ].myLength = 1e100;
546 while ( paths.size() >= 2 )
548 for ( size_t i = 0; i < paths.size(); ++i )
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
554 Path::Remove( paths, i );
558 // join paths that reach same point
559 for ( size_t j = 0; j < paths.size(); ++j )
561 if ( i != j && paths[i].ReachSamePoint( paths[j] ))
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 )
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() );
574 if ( i < j ) std::swap( i, j );
575 Path::Remove( paths, i );
576 Path::Remove( paths, j );
581 if ( !paths.empty() && (int) paths[0].myPoints.size() > myMesh->NbFaces() )
582 throw SALOME_Exception(LOCALIZED( "Infinite loop in MakePolyLine()"));
585 if ( myPaths[ iSeg ].myPoints.empty() )
586 throw SALOME_Exception( SMESH_Comment("Can't find a full path for PolySegment #") << iSeg );
589 double d00 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.front() ).SquareModulus();
590 double d01 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.back() ).SquareModulus();
592 std::reverse( myPaths[ iSeg ].myPoints.begin(), myPaths[ iSeg ].myPoints.end() );
594 } // PolyPathCompute::Compute()
596 }; // struct PolyPathCompute
600 //=======================================================================
601 //function : MakePolyLine
602 //purpose : Create a polyline consisting of 1D mesh elements each lying on a 2D element of
604 //=======================================================================
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)
613 std::vector< Path > segPaths( theSegments.size() ); // path of each of segments
615 SMESH_ElementSearcher* searcher = theSearcher;
616 SMESHUtils::Deleter<SMESH_ElementSearcher> delSearcher;
619 searcher = SMESH_MeshAlgos::GetElementSearcher( *theMesh );
620 delSearcher._obj = searcher;
623 // get cutting planes
625 std::vector< bool > isVectorOK( theSegments.size(), true );
626 const double planarCoef = 0.333; // plane height in planar case
628 for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
630 PolySegment& polySeg = theSegments[ iSeg ];
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] ));
637 polySeg.myFace[0] = polySeg.myFace[1] = 0;
638 if ( !polySeg.myNode1[0] && !polySeg.myNode2[0] )
640 p1 = searcher->Project( polySeg.myXYZ[0], SMDSAbs_Face, &polySeg.myFace[0] );
642 if ( !polySeg.myNode1[1] && !polySeg.myNode2[1] )
644 p2 = searcher->Project( polySeg.myXYZ[1], SMDSAbs_Face, &polySeg.myFace[1] );
646 polySeg.myXYZ[0] = p1;
647 polySeg.myXYZ[1] = p2;
649 gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
651 isVectorOK[ iSeg ] = ( plnNorm.Modulus() > std::numeric_limits<double>::min() );
652 if ( !isVectorOK[ iSeg ] && ( p1 - p2 ).SquareModulus() > 0. )
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;
660 SMESH_MeshAlgos::FaceNormal( face, faceNorm, /*normalized=*/false );
662 const double tol = Precision::Confusion();
663 if ( polySeg.myVector.Magnitude() < tol || polySeg.myVector * faceNorm < tol )
665 polySeg.myVector = faceNorm;
666 polySeg.myMidProjPoint = pMid + faceNorm * ( p1 - p2 ).Modulus() * planarCoef;
668 plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
669 if ( plnNorm.SquareModulus() == 0 ) // p1-p2 perpendicular to mesh
671 double radius = faceNorm.Modulus();
672 std::vector< const SMDS_MeshElement* > foundElems;
673 while ( plnNorm.SquareModulus() == 0 && radius < 1e200 )
676 searcher->GetElementsInSphere( p1, radius, SMDSAbs_Face, foundElems );
677 searcher->GetElementsInSphere( p2, radius, SMDSAbs_Face, foundElems );
679 polySeg.myVector.SetCoord( 0,0,0 );
680 for ( size_t i = 0; i < foundElems.size(); ++i )
682 SMESH_MeshAlgos::FaceNormal( foundElems[i], faceNorm );
683 polySeg.myVector += faceNorm / foundElems.size();
685 plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
692 polySeg.myVector = plnNorm ^ ( p1 - p2 );
696 // assure that inverse elements are constructed, avoid their concurrent building in threads
697 theMesh->nodesIterator()->next()->NbInverseElements();
701 PolyPathCompute algo( theSegments, segPaths, theMesh );
702 OSD_Parallel::For( 0, theSegments.size(), algo, theSegments.size() == 1 );
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() );
710 const SMDS_MeshNode *n, *nPrev = 0;
712 for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
714 const Path& path = segPaths[iSeg];
715 if ( path.myPoints.size() < 2 )
718 double tol = path.myLength / path.myPoints.size() / 1000.;
719 if ( !nPrev || ( SMESH_NodeXYZ( nPrev ) - path.myPoints[0] ).SquareModulus() > tol*tol )
721 nPrev = theMesh->AddNode( path.myPoints[0].X(), path.myPoints[0].Y(), path.myPoints[0].Z() );
722 theNewNodes.push_back( nPrev );
724 for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
726 n = theMesh->AddNode( path.myPoints[iP].X(), path.myPoints[iP].Y(), path.myPoints[iP].Z() );
727 theNewNodes.push_back( n );
729 const SMDS_MeshElement* elem = theMesh->AddEdge( nPrev, n );
730 theNewEdges.push_back( elem );
732 theGroup->Add( elem );
739 gp_XYZ pMid = 0.5 * ( path.myPoints[0] + path.myPoints.back() );
740 if ( isVectorOK[ iSeg ])
742 // find the most distant point of a path
744 for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
746 double dist = Abs( theSegments[iSeg].myVector * ( path.myPoints[iP] - path.myPoints[0] ));
747 if ( dist > maxDist )
750 theSegments[iSeg].myMidProjPoint = path.myPoints[iP];
753 if ( maxDist < Precision::Confusion() ) // planar case
754 theSegments[iSeg].myMidProjPoint =
755 pMid + theSegments[iSeg].myVector.XYZ().Normalized() * path.myLength * planarCoef;
757 theSegments[iSeg].myVector = gp_Vec( pMid, theSegments[iSeg].myMidProjPoint );