1 // Copyright (C) 2018 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(): myLength(0.0), myFace(0) {}
57 bool SetCutAtCorner( const SMESH_NodeXYZ& cornerNode,
58 const SMDS_MeshElement* face,
59 const gp_XYZ& plnNorm,
60 const gp_XYZ& plnOrig );
62 void AddPoint( const gp_XYZ& p );
64 bool Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig );
66 bool ReachSamePoint( const Path& other );
68 static void Remove( std::vector< Path > & paths, size_t& i );
71 //================================================================================
73 * \brief Return true if this Path meats another
75 //================================================================================
77 bool Path::ReachSamePoint( const Path& other )
79 return ( mySrcPntInd != other.mySrcPntInd &&
80 myFace == other.myFace );
83 //================================================================================
85 * \brief Remove a path from a vector
87 //================================================================================
89 void Path::Remove( std::vector< Path > & paths, size_t& i )
91 if ( paths.size() > 1 )
93 size_t j = paths.size() - 1; // last item to be removed
96 paths[ i ].myPoints.swap( paths[ j ].myPoints );
97 paths[ i ].myLength = paths[ j ].myLength;
98 paths[ i ].mySrcPntInd = paths[ j ].mySrcPntInd;
99 paths[ i ].myFace = paths[ j ].myFace;
100 paths[ i ].myNode1 = paths[ j ].myNode1;
101 paths[ i ].myNode2 = paths[ j ].myNode2;
102 paths[ i ].myNodeInd1 = paths[ j ].myNodeInd1;
103 paths[ i ].myNodeInd2 = paths[ j ].myNodeInd2;
104 paths[ i ].myDot1 = paths[ j ].myDot1;
105 paths[ i ].myDot2 = paths[ j ].myDot2;
113 //================================================================================
115 * \brief Store a point that is at a node of a face if the face is intersected by plane.
116 * Return false if the node is a sole intersection point of the face and the plane
118 //================================================================================
120 bool Path::SetCutAtCorner( const SMESH_NodeXYZ& cornerNode,
121 const SMDS_MeshElement* face,
122 const gp_XYZ& plnNorm,
123 const gp_XYZ& plnOrig )
125 if ( face == myFace )
127 myNodeInd1 = face->GetNodeIndex( cornerNode._node );
128 myNodeInd2 = ( myNodeInd1 + 1 ) % face->NbCornerNodes();
129 int ind3 = ( myNodeInd1 + 2 ) % face->NbCornerNodes();
130 myNode1.Set( face->GetNode( ind3 ));
131 myNode2.Set( face->GetNode( myNodeInd2 ));
133 myDot1 = plnNorm * ( myNode1 - plnOrig );
134 myDot2 = plnNorm * ( myNode2 - plnOrig );
136 bool ok = ( myDot1 * myDot2 < 0 );
137 if ( !ok && myDot1 * myDot2 == 0 )
139 ok = ( myDot1 != myDot2 );
141 ok = ( myFace->GetNodeIndex(( myDot1 == 0 ? myNode1 : myNode2 )._node ) < 0 );
147 AddPoint( cornerNode );
152 //================================================================================
154 * \brief Store a point and update myLength
156 //================================================================================
158 void Path::AddPoint( const gp_XYZ& p )
160 if ( !myPoints.empty() )
161 myLength += ( p - myPoints.back() ).Modulus();
164 myPoints.push_back( p );
167 //================================================================================
169 * \brief Try to find the next point
170 * \param [in] plnNorm - cutting plane normal
171 * \param [in] plnOrig - cutting plane origin
173 //================================================================================
175 bool Path::Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig )
177 int nodeInd3 = ( myNodeInd1 + 1 ) % myFace->NbCornerNodes();
178 if ( myNodeInd2 == nodeInd3 )
179 nodeInd3 = ( myNodeInd1 + 2 ) % myFace->NbCornerNodes();
181 SMESH_NodeXYZ node3 = myFace->GetNode( nodeInd3 );
182 double dot3 = plnNorm * ( node3 - plnOrig );
184 if ( dot3 * myDot1 < 0. )
187 myNodeInd2 = nodeInd3;
190 else if ( dot3 * myDot2 < 0. )
193 myNodeInd1 = nodeInd3;
196 else if ( dot3 == 0. )
198 SMDS_ElemIteratorPtr fIt = node3._node->GetInverseElementIterator(SMDSAbs_Face);
199 while ( fIt->more() )
200 if ( SetCutAtCorner( node3, fIt->next(), plnNorm, plnOrig ))
204 else if ( myDot2 == 0. )
206 SMESH_NodeXYZ node2 = myNode2; // copy as myNode2 changes in SetCutAtCorner()
207 SMDS_ElemIteratorPtr fIt = node2._node->GetInverseElementIterator(SMDSAbs_Face);
208 while ( fIt->more() )
209 if ( SetCutAtCorner( node2, fIt->next(), plnNorm, plnOrig ))
214 double r = Abs( myDot1 / ( myDot2 - myDot1 ));
215 AddPoint( myNode1 * ( 1 - r ) + myNode2 * r );
218 myAvoidSet.insert( myFace );
219 myFace = SMESH_MeshAlgos::FindFaceInSet( myNode1._node, myNode2._node,
220 myElemSet, myAvoidSet,
221 &myNodeInd1, &myNodeInd2 );
225 //================================================================================
227 * \brief Compute a path between two points of PolySegment
229 struct PolyPathCompute
231 SMESH_MeshAlgos::TListOfPolySegments& mySegments; //!< inout PolySegment's
232 std::vector< Path >& myPaths; //!< path of each of segments to compute
234 mutable std::vector< std::string > myErrors;
236 PolyPathCompute( SMESH_MeshAlgos::TListOfPolySegments& theSegments,
237 std::vector< Path >& thePaths,
239 mySegments( theSegments ),
242 myErrors( theSegments.size() )
247 #define SMESH_CAUGHT myErrors[i] =
248 void operator() ( const int i ) const
251 const_cast< PolyPathCompute* >( this )->Compute( i );
252 SMESH_CATCH( SMESH::returnError );
256 //================================================================================
258 * \brief Compute a path of a given segment
260 //================================================================================
262 void Compute( const int iSeg )
264 SMESH_MeshAlgos::PolySegment& polySeg = mySegments[ iSeg ];
267 gp_XYZ plnNorm = ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ) ^ polySeg.myVector.XYZ();
268 gp_XYZ plnOrig = polySeg.myXYZ[1];
270 // Find paths connecting the 2 end points of polySeg
272 std::vector< Path > paths; paths.reserve(10);
274 // 1) initialize paths; two paths starts at each end point
276 for ( int iP = 0; iP < 2; ++iP ) // loop on the polySeg end points
279 path.mySrcPntInd = iP;
280 size_t nbPaths = paths.size();
282 if ( polySeg.myFace[ iP ]) // the end point lies on polySeg.myFace[ iP ]
284 // check coincidence of polySeg.myXYZ[ iP ] with nodes
285 const double tol = 1e-20;
286 SMESH_NodeXYZ nodes[4];
287 for ( int i = 0; i < 3 && !polySeg.myNode1[ iP ]; ++i )
289 nodes[ i ] = polySeg.myFace[ iP ]->GetNode( i );
290 if (( nodes[ i ] - polySeg.myXYZ[ iP ]).SquareModulus() < tol*tol )
291 polySeg.myNode1[ iP ] = nodes[ i ].Node();
293 nodes[ 3 ] = nodes[ 0 ];
295 // check coincidence of polySeg.myXYZ[ iP ] with edges
296 for ( int i = 0; i < 3 && !polySeg.myNode1[ iP ]; ++i )
298 SMDS_LinearEdge edge( nodes[i].Node(), nodes[i+1].Node() );
299 if ( SMESH_MeshAlgos::GetDistance( &edge, polySeg.myXYZ[ iP ]) < tol )
301 polySeg.myNode1[ iP ] = nodes[ i ].Node();
302 polySeg.myNode2[ iP ] = nodes[ i + 1 ].Node();
306 if ( !polySeg.myNode1[ iP ] ) // polySeg.myXYZ[ iP ] is within polySeg.myFace[ iP ]
309 for ( int i = 0; i < 3; ++i )
310 dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
313 int iCut = 0; // index of a cut edge
314 if ( dot[ 1 ] * dot[ 2 ] < 0. ) iCut = 1;
315 else if ( dot[ 2 ] * dot[ 3 ] < 0. ) iCut = 2;
317 // initialize path so as if it entered the face via iCut-th edge
318 path.myFace = polySeg.myFace[ iP ];
319 path.myNodeInd1 = iCut;
320 path.myNodeInd2 = iCut + 1;
321 path.myNode1.Set( nodes[ iCut ].Node() );
322 path.myNode2.Set( nodes[ iCut + 1 ].Node() );
323 path.myDot1 = dot[ iCut ];
324 path.myDot2 = dot[ iCut + 1 ];
325 path.myPoints.clear();
326 path.AddPoint( polySeg.myXYZ[ iP ]);
327 paths.push_back( path );
329 path.Extend( plnNorm, plnOrig ); // to get another edge cut
330 path.myFace = polySeg.myFace[ iP ];
331 if ( path.myDot1 == 0. ) // cut at a node
333 path.myNodeInd1 = ( iCut + 2 ) % 3;
334 path.myNodeInd2 = ( iCut + 3 ) % 3;
335 path.myNode2.Set( path.myFace->GetNode( path.myNodeInd2 ));
336 path.myDot2 = dot[ path.myNodeInd2 ];
340 path.myNodeInd1 = path.myFace->GetNodeIndex( path.myNode1.Node() );
341 path.myNodeInd2 = path.myFace->GetNodeIndex( path.myNode2.Node() );
343 path.myPoints.clear();
344 path.AddPoint( polySeg.myXYZ[ iP ]);
345 paths.push_back( path );
349 if ( polySeg.myNode2[ iP ] && polySeg.myNode2[ iP ] != polySeg.myNode1[ iP ] )
351 // the end point is on an edge
352 while (( path.myFace = SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
353 polySeg.myNode2[ iP ],
359 path.myNode1.Set( polySeg.myNode1[ iP ]);
360 path.myNode2.Set( polySeg.myNode2[ iP ]);
361 path.myDot1 = plnNorm * ( path.myNode1 - plnOrig );
362 path.myDot2 = plnNorm * ( path.myNode2 - plnOrig );
363 path.myPoints.clear();
364 path.AddPoint( polySeg.myXYZ[ iP ]);
365 path.myAvoidSet.insert( path.myFace );
366 paths.push_back( path );
368 if ( nbPaths == paths.size() )
369 throw SALOME_Exception ( SMESH_Comment("No face edge found by point ") << iP+1
370 << " in a PolySegment " << iSeg );
372 else if ( polySeg.myNode1[ iP ] ) // the end point is at a node
374 std::set<const SMDS_MeshNode* > nodes;
375 SMDS_ElemIteratorPtr fIt = polySeg.myNode1[ iP ]->GetInverseElementIterator(SMDSAbs_Face);
376 while ( fIt->more() )
378 path.myPoints.clear();
379 if ( path.SetCutAtCorner( polySeg.myNode1[ iP ], fIt->next(), plnNorm, plnOrig ))
381 if (( path.myDot1 * path.myDot2 != 0 ) ||
382 ( nodes.insert( path.myDot1 == 0 ? path.myNode1._node : path.myNode2._node ).second ))
383 paths.push_back( path );
388 // look for a one-segment path
389 for ( size_t i = 0; i < nbPaths; ++i )
390 for ( size_t j = nbPaths; j < paths.size(); ++j )
391 if ( paths[i].myFace == paths[j].myFace )
393 myPaths[ iSeg ].myPoints.push_back( paths[i].myPoints[0] );
394 myPaths[ iSeg ].myPoints.push_back( paths[j].myPoints[0] );
399 // 2) extend paths and compose the shortest one connecting the two points
401 myPaths[ iSeg ].myLength = 1e100;
403 while ( paths.size() >= 2 )
405 for ( size_t i = 0; i < paths.size(); ++i )
407 Path& path = paths[ i ];
408 if ( !path.Extend( plnNorm, plnOrig ) || // path reached a mesh boundary
409 path.myLength > myPaths[ iSeg ].myLength ) // path is longer than others
411 Path::Remove( paths, i );
415 // join paths that reach same point
416 for ( size_t j = 0; j < paths.size(); ++j )
418 if ( i != j && paths[i].ReachSamePoint( paths[j] ))
420 double distLast = ( paths[i].myPoints.back() - paths[j].myPoints.back() ).Modulus();
421 double fullLength = ( paths[i].myLength + paths[j].myLength + distLast );
422 if ( fullLength < myPaths[ iSeg ].myLength )
424 myPaths[ iSeg ].myLength = fullLength;
425 std::vector< gp_XYZ > & allPoints = myPaths[ iSeg ].myPoints;
426 allPoints.swap( paths[i].myPoints );
427 allPoints.insert( allPoints.end(),
428 paths[j].myPoints.rbegin(),
429 paths[j].myPoints.rend() );
431 Path::Remove( paths, i );
432 Path::Remove( paths, j );
436 if ( !paths.empty() && (int) paths[0].myPoints.size() > myMesh->NbFaces() )
437 throw SALOME_Exception(LOCALIZED( "Infinite loop in MakePolyLine()"));
440 if ( myPaths[ iSeg ].myPoints.empty() )
441 throw SALOME_Exception( SMESH_Comment("Can't find a full path for PolySegment #") << iSeg );
444 double d00 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.front() ).SquareModulus();
445 double d01 = ( polySeg.myXYZ[0] - myPaths[ iSeg ].myPoints.back() ).SquareModulus();
447 std::reverse( myPaths[ iSeg ].myPoints.begin(), myPaths[ iSeg ].myPoints.end() );
449 } // PolyPathCompute::Compute()
451 }; // struct PolyPathCompute
455 //=======================================================================
456 //function : MakePolyLine
457 //purpose : Create a polyline consisting of 1D mesh elements each lying on a 2D element of
459 //=======================================================================
461 void SMESH_MeshAlgos::MakePolyLine( SMDS_Mesh* theMesh,
462 TListOfPolySegments& theSegments,
463 std::vector<const SMDS_MeshElement*>& theNewEdges,
464 std::vector< const SMDS_MeshNode* >& theNewNodes,
465 SMDS_MeshGroup* theGroup,
466 SMESH_ElementSearcher* theSearcher)
468 std::vector< Path > segPaths( theSegments.size() ); // path of each of segments
470 SMESH_ElementSearcher* searcher = theSearcher;
471 SMESHUtils::Deleter<SMESH_ElementSearcher> delSearcher;
474 searcher = SMESH_MeshAlgos::GetElementSearcher( *theMesh );
475 delSearcher._obj = searcher;
478 // get cutting planes
480 std::vector< bool > isVectorOK( theSegments.size(), true );
481 const double planarCoef = 0.333; // plane height in planar case
483 for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
485 PolySegment& polySeg = theSegments[ iSeg ];
487 gp_XYZ p1 = SMESH_NodeXYZ( polySeg.myNode1[0] );
488 gp_XYZ p2 = SMESH_NodeXYZ( polySeg.myNode1[1] );
489 if ( polySeg.myNode2[0] ) p1 = 0.5 * ( p1 + SMESH_NodeXYZ( polySeg.myNode2[0] ));
490 if ( polySeg.myNode2[1] ) p2 = 0.5 * ( p2 + SMESH_NodeXYZ( polySeg.myNode2[1] ));
492 polySeg.myFace[0] = polySeg.myFace[1] = 0;
493 if ( !polySeg.myNode1[0] && !polySeg.myNode2[0] )
495 p1 = searcher->Project( polySeg.myXYZ[0], SMDSAbs_Face, &polySeg.myFace[0] );
497 if ( !polySeg.myNode1[1] && !polySeg.myNode2[1] )
499 p2 = searcher->Project( polySeg.myXYZ[1], SMDSAbs_Face, &polySeg.myFace[1] );
501 polySeg.myXYZ[0] = p1;
502 polySeg.myXYZ[1] = p2;
504 gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
506 isVectorOK[ iSeg ] = ( plnNorm.Modulus() > std::numeric_limits<double>::min() );
507 if ( !isVectorOK[ iSeg ])
509 gp_XYZ pMid = 0.5 * ( p1 + p2 );
510 const SMDS_MeshElement* face;
511 polySeg.myMidProjPoint = searcher->Project( pMid, SMDSAbs_Face, &face );
512 polySeg.myVector = polySeg.myMidProjPoint.XYZ() - pMid;
515 SMESH_MeshAlgos::FaceNormal( face, faceNorm );
517 if ( polySeg.myVector.Magnitude() < Precision::Confusion() ||
518 polySeg.myVector * faceNorm < Precision::Confusion() )
520 polySeg.myVector = faceNorm;
521 polySeg.myMidProjPoint = pMid + faceNorm * ( p1 - p2 ).Modulus() * planarCoef;
526 polySeg.myVector = plnNorm ^ ( p1 - p2 );
530 // assure that inverse elements are constructed, avoid their concurrent building in threads
531 theMesh->nodesIterator()->next()->NbInverseElements();
535 PolyPathCompute algo( theSegments, segPaths, theMesh );
536 OSD_Parallel::For( 0, theSegments.size(), algo, theSegments.size() == 1 );
538 for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
539 if ( !algo.myErrors[ iSeg ].empty() )
540 throw SALOME_Exception( algo.myErrors[ iSeg ].c_str() );
544 const SMDS_MeshNode *n, *nPrev = 0;
546 for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
548 const Path& path = segPaths[iSeg];
549 if ( path.myPoints.size() < 2 )
552 double tol = path.myLength / path.myPoints.size() / 1000.;
553 if ( !nPrev || ( SMESH_NodeXYZ( nPrev ) - path.myPoints[0] ).SquareModulus() > tol*tol )
555 nPrev = theMesh->AddNode( path.myPoints[0].X(), path.myPoints[0].Y(), path.myPoints[0].Z() );
556 theNewNodes.push_back( nPrev );
558 for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
560 n = theMesh->AddNode( path.myPoints[iP].X(), path.myPoints[iP].Y(), path.myPoints[iP].Z() );
561 theNewNodes.push_back( n );
563 const SMDS_MeshElement* elem = theMesh->AddEdge( nPrev, n );
564 theNewEdges.push_back( elem );
566 theGroup->Add( elem );
573 gp_XYZ pMid = 0.5 * ( path.myPoints[0] + path.myPoints.back() );
574 if ( isVectorOK[ iSeg ])
576 // find the most distant point of a path
578 for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
580 double dist = Abs( theSegments[iSeg].myVector * ( path.myPoints[iP] - path.myPoints[0] ));
581 if ( dist > maxDist )
584 theSegments[iSeg].myMidProjPoint = path.myPoints[iP];
587 if ( maxDist < Precision::Confusion() ) // planar case
588 theSegments[iSeg].myMidProjPoint =
589 pMid + theSegments[iSeg].myVector.XYZ().Normalized() * path.myLength * planarCoef;
591 theSegments[iSeg].myVector = gp_Vec( pMid, theSegments[iSeg].myMidProjPoint );