PAL0023627: [IMACS] ASERIS: project point to the mesh

Algorithm of poly-line construction modified to enable specifying
   arbitrary points serving as end points of the poly-line segment.

doc/salome/gui/SMESH/input/smesh_module.rst

@@ -618,7 +618,7 @@ SMESH_MeshEditor
 
 .. py:class:: SMESH_MeshEditor.Sew_Error
 
-      Enumeration of errors SMESH_MeshEditor.Sewing... methods
+      Enumeration of errors of SMESH_MeshEditor.Sewing... methods
 
       .. py:attribute::
          SEW_OK

idl/SMESH_MeshEditor.idl

@@ -60,17 +60,23 @@ module SMESH
   // structure used in MakePolyLine() to define a cutting plane
   struct PolySegment
   {
-    // point 1: if node1ID2 > 0, then the point is in the middle of a face edge defined
+    // a point is defined as follows:
-    //          by two nodes, else it is at node1ID1
+    // ( node*ID1 > 0 && node*ID2 > 0 ) ==> point is in the middle of an edge defined by two nodes
+    // ( node*ID1 > 0 && node*ID2 <=0 ) ==> point is at node*ID1
+    // else                             ==> point is at xyz*
+
+    // point 1
     long node1ID1;
     long node1ID2;
+    PointStruct xyz1;
 
-    // point 2: if node2ID2 > 0, then the point is in the middle of a face edge defined
+    // point 2
-    //          by two nodes, else it is at node2ID1
     long node2ID1;
     long node2ID2;
+    PointStruct xyz2;
 
-    DirStruct vector; // vector on the plane; to use a default plane set vector = (0,0,0)
+    // vector on the plane; to use a default plane set vector = (0,0,0)
+    DirStruct vector;
   };
   typedef sequence<PolySegment> ListOfPolySegments;
 
@@ -1255,8 +1261,8 @@ module SMESH
 
     /*!
      * \brief Create a polyline consisting of 1D mesh elements each lying on a 2D element of
-     *        the initial mesh. Positions of new nodes are found by cutting the mesh by the
+     *        the initial triangle mesh. Positions of new nodes are found by cutting the mesh by
-     *        plane passing through pairs of points specified by each PolySegment structure.
+     *        the plane passing through pairs of points specified by each PolySegment structure.
      *        If there are several paths connecting a pair of points, the shortest path is
      *        selected by the module. Position of the cutting plane is defined by the two
      *        points and an optional vector lying on the plane specified by a PolySegment.

src/SMESH/SMESH_MeshEditor.cxx

@@ -13050,14 +13050,15 @@ namespace // utils for MakePolyLine
     std::vector< gp_XYZ >   myPoints;
     double                  myLength;
 
-    int                     mySrcPntInd; //!< start point index
     const SMDS_MeshElement* myFace;
-    SMESH_NodeXYZ           myNode1;
+    SMESH_NodeXYZ           myNode1; // nodes of the edge the path entered myFace
     SMESH_NodeXYZ           myNode2;
     int                     myNodeInd1;
     int                     myNodeInd2;
     double                  myDot1;
     double                  myDot2;
+
+    int                     mySrcPntInd; //!< start point index
     TIDSortedElemSet        myElemSet, myAvoidSet;
 
     Path(): myLength(0.0), myFace(0) {}
@@ -13250,6 +13251,7 @@ namespace // utils for MakePolyLine
       myErrors( theSegments.size() )
     {
     }
+
 #undef SMESH_CAUGHT
 #define SMESH_CAUGHT myErrors[i] =
     void operator() ( const int i ) const
@@ -13259,6 +13261,7 @@ namespace // utils for MakePolyLine
       SMESH_CATCH( SMESH::returnError );
     }
 #undef SMESH_CAUGHT
+
     //================================================================================
     /*!
      * \brief Compute a path of a given segment
@@ -13269,21 +13272,15 @@ namespace // utils for MakePolyLine
     {
       SMESH_MeshEditor::PolySegment& polySeg = mySegments[ iSeg ];
 
-      // get a cutting plane
+      // the cutting plane
+      gp_XYZ plnNorm = ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ) ^ polySeg.myVector.XYZ();
+      gp_XYZ plnOrig = polySeg.myXYZ[1];
 
-      gp_XYZ p1 = SMESH_NodeXYZ( polySeg.myNode1[0] );
+      // Find paths connecting the 2 end points of polySeg
-      gp_XYZ p2 = SMESH_NodeXYZ( polySeg.myNode1[1] );
-      if ( polySeg.myNode2[0] ) p1 = 0.5 * ( p1 + SMESH_NodeXYZ( polySeg.myNode2[0] ));
-      if ( polySeg.myNode2[1] ) p2 = 0.5 * ( p2 + SMESH_NodeXYZ( polySeg.myNode2[1] ));
-
-      gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
-      gp_XYZ plnOrig = p2;
-
-      // find paths connecting the 2 end points of polySeg
 
       std::vector< Path > paths; paths.reserve(10);
 
-      // initialize paths
+      // 1) initialize paths; two paths starts at each end point
 
       for ( int iP = 0; iP < 2; ++iP ) // loop on the polySeg end points
       {
@@ -13291,8 +13288,76 @@ namespace // utils for MakePolyLine
         path.mySrcPntInd = iP;
         size_t nbPaths = paths.size();
 
+        if ( polySeg.myFace[ iP ]) // the end point lies on polySeg.myFace[ iP ]
+        {
+          // check coincidence of polySeg.myXYZ[ iP ] with nodes
+          const double tol = 1e-20;
+          SMESH_NodeXYZ nodes[4];
+          for ( int i = 0; i < 3 &&  !polySeg.myNode1[ iP ]; ++i )
+          {
+            nodes[ i ] = polySeg.myFace[ iP ]->GetNode( i );
+            if (( nodes[ i ] - polySeg.myXYZ[ iP ]).SquareModulus() < tol*tol )
+              polySeg.myNode1[ iP ] = nodes[ i ].Node();
+          }
+          nodes[ 3 ] = nodes[ 0 ];
+
+          // check coincidence of polySeg.myXYZ[ iP ] with edges
+          for ( int i = 0; i < 3 &&  !polySeg.myNode1[ iP ]; ++i )
+          {
+            SMDS_LinearEdge edge( nodes[i].Node(), nodes[i+1].Node() );
+            if ( SMESH_MeshAlgos::GetDistance( &edge, polySeg.myXYZ[ iP ]) < tol )
+            {
+              polySeg.myNode1[ iP ] = nodes[ i     ].Node();
+              polySeg.myNode2[ iP ] = nodes[ i + 1 ].Node();
+            }
+          }
+
+          if ( !polySeg.myNode1[ iP ] ) // polySeg.myXYZ[ iP ] is within polySeg.myFace[ iP ]
+          {
+            double dot[ 4 ];
+            for ( int i = 0; i < 3; ++i )
+              dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
+            dot[ 3 ] = dot[ 0 ];
+
+            int iCut = 0; // index of a cut edge
+            if      ( dot[ 1 ] * dot[ 2 ] < 0. ) iCut = 1;
+            else if ( dot[ 2 ] * dot[ 3 ] < 0. ) iCut = 2;
+
+            // initialize path so as if it entered the face via iCut-th edge
+            path.myFace = polySeg.myFace[ iP ];
+            path.myNodeInd1 = iCut;
+            path.myNodeInd2 = iCut + 1;
+            path.myNode1.Set( nodes[ iCut     ].Node() );
+            path.myNode2.Set( nodes[ iCut + 1 ].Node() );
+            path.myDot1 = dot[ iCut ];
+            path.myDot2 = dot[ iCut + 1 ];
+            path.myPoints.clear();
+            path.AddPoint( polySeg.myXYZ[ iP ]);
+            paths.push_back( path );
+
+            path.Extend( plnNorm, plnOrig ); // to get another edge cut
+            path.myFace = polySeg.myFace[ iP ];
+            if ( path.myDot1 == 0. ) // cut at a node
+            {
+              path.myNodeInd1 = ( iCut + 2 ) % 3;
+              path.myNodeInd2 = ( iCut + 3 ) % 3;
+              path.myNode2.Set( path.myFace->GetNode( path.myNodeInd2 ));
+              path.myDot2 = dot[ path.myNodeInd2 ];
+            }
+            else
+            {
+              path.myNodeInd1 = path.myFace->GetNodeIndex( path.myNode1.Node() );
+              path.myNodeInd2 = path.myFace->GetNodeIndex( path.myNode2.Node() );
+            }
+            path.myPoints.clear();
+            path.AddPoint( polySeg.myXYZ[ iP ]);
+            paths.push_back( path );
+          }
+        }
+
         if ( polySeg.myNode2[ iP ] && polySeg.myNode2[ iP ] != polySeg.myNode1[ iP ] )
         {
+          // the end point is on an edge
           while (( path.myFace = SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
                                                                  polySeg.myNode2[ iP ],
                                                                  path.myElemSet,
@@ -13305,7 +13370,7 @@ namespace // utils for MakePolyLine
             path.myDot1 = plnNorm * ( path.myNode1 - plnOrig );
             path.myDot2 = plnNorm * ( path.myNode2 - plnOrig );
             path.myPoints.clear();
-            path.AddPoint( 0.5 * ( path.myNode1 + path.myNode2 ));
+            path.AddPoint( polySeg.myXYZ[ iP ]);
             path.myAvoidSet.insert( path.myFace );
             paths.push_back( path );
           }
@@ -13313,7 +13378,7 @@ namespace // utils for MakePolyLine
             throw SALOME_Exception ( SMESH_Comment("No face edge found by point ") << iP+1
                                      << " in a PolySegment " << iSeg );
         }
-        else // an end point is at node
+        else if ( polySeg.myNode1[ iP ] ) // the end point is at a node
         {
           std::set<const SMDS_MeshNode* > nodes;
           SMDS_ElemIteratorPtr fIt = polySeg.myNode1[ iP ]->GetInverseElementIterator(SMDSAbs_Face);
@@ -13340,7 +13405,7 @@ namespace // utils for MakePolyLine
             }
       }
 
-      // extend paths
+      // 2) extend paths and compose the shortest one connecting the two points
 
       myPaths[ iSeg ].myLength = 1e100;
 
@@ -13349,7 +13414,7 @@ namespace // utils for MakePolyLine
         for ( size_t i = 0; i < paths.size(); ++i )
         {
           Path& path = paths[ i ];
-          if ( !path.Extend( plnNorm, plnOrig ) ||         // path reached a mesh boundary
+          if ( !path.Extend( plnNorm, plnOrig ) ||        // path reached a mesh boundary
                path.myLength > myPaths[ iSeg ].myLength ) // path is longer than others
           {
             Path::Remove( paths, i );
@@ -13424,6 +13489,18 @@ void SMESH_MeshEditor::MakePolyLine( TListOfPolySegments&   theSegments,
     if ( polySeg.myNode2[0] ) p1 = 0.5 * ( p1 + SMESH_NodeXYZ( polySeg.myNode2[0] ));
     if ( polySeg.myNode2[1] ) p2 = 0.5 * ( p2 + SMESH_NodeXYZ( polySeg.myNode2[1] ));
 
+    polySeg.myFace[0] = polySeg.myFace[1] = 0;
+    if ( !polySeg.myNode1[0] && !polySeg.myNode2[0] )
+    {
+      p1 = searcher->Project( polySeg.myXYZ[0], SMDSAbs_Face, &polySeg.myFace[0] );
+    }
+    if ( !polySeg.myNode1[1] && !polySeg.myNode2[1] )
+    {
+      p2 = searcher->Project( polySeg.myXYZ[1], SMDSAbs_Face, &polySeg.myFace[1] );
+    }
+    polySeg.myXYZ[0] = p1;
+    polySeg.myXYZ[1] = p2;
+
     gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
 
     isVectorOK[ iSeg ] = ( plnNorm.Modulus() > std::numeric_limits<double>::min() );

src/SMESH/SMESH_MeshEditor.hxx

@@ -721,12 +721,19 @@ public:
   // structure used in MakePolyLine() to define a cutting plane
   struct PolySegment
   {
-    // 2 points: if myNode2 != 0, then the point is the middle of a face edge defined
+    // 2 points, each defined as follows:
-    //           by two nodes, else it is at myNode1
+    // ( myNode1 &&  myNode2 ) ==> point is in the middle of an edge defined by two nodes
-    const SMDS_MeshNode* myNode1[2];
+    // ( myNode1 && !myNode2 ) ==> point is at myNode1
-    const SMDS_MeshNode* myNode2[2];
+    // else                    ==> point is at myXYZ
+    const SMDS_MeshNode*    myNode1[2];
+    const SMDS_MeshNode*    myNode2[2];
+    gp_XYZ                  myXYZ  [2];
 
-    gp_Vec myVector; // vector on the plane; to use a default plane set vector = (0,0,0)
+    // face on which myXYZ projects (found by MakePolyLine())
+    const SMDS_MeshElement* myFace [2];
+
+    // vector on the plane; to use a default plane set vector = (0,0,0)
+    gp_Vec myVector;
 
     // point to return coordinates of a middle of the two points, projected to mesh
     gp_Pnt myMidProjPoint;

src/SMESHUtils/SMESH_MeshAlgos.cxx

@@ -1476,9 +1476,9 @@ namespace
 
   //================================================================================
   /*!
-   * \brief Return of a point relative to a segment
+   * \brief Return position of a point relative to a segment
    *  \param point2D      - the point to analyze position of
-   *  \param xyVec        - end points of segments
+   *  \param segEnds      - end points of segments
    *  \param index0       - 0-based index of the first point of segment
    *  \param posToFindOut - flags of positions to detect
    *  \retval PointPos - point position
@@ -1607,46 +1607,63 @@ double SMESH_MeshAlgos::GetDistance( const SMDS_MeshFace* face,
   }
 
   // compute distance
-  PointPos pos = *pntPosSet.begin();
+
-  switch ( pos._name )
+  double minDist2 = Precision::Infinite();
+  for ( std::set< PointPos >::iterator posIt = pntPosSet.begin(); posIt != pntPosSet.end(); ++posIt)
   {
-  case POS_LEFT:
+    PointPos pos = *posIt;
-  {
+    if ( pos._name != pntPosSet.begin()->_name )
-    // point is most close to an edge
+      break;
-    gp_Vec edge( xyz[ pos._index ], xyz[ pos._index+1 ]);
+    switch ( pos._name )
-    gp_Vec n1p ( xyz[ pos._index ], point  );
-    double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
-    // projection of the point on the edge
-    gp_XYZ proj = xyz[ pos._index ] + u * edge.XYZ();
-    if ( closestPnt ) *closestPnt = proj;
-    return point.Distance( proj );
-  }
-  case POS_RIGHT:
-  {
-    // point is inside the face
-    double distToFacePlane = Abs( tmpPnt.Y() );
-    if ( closestPnt )
     {
-      if ( distToFacePlane < std::numeric_limits<double>::min() ) {
+    case POS_LEFT: // point is most close to an edge
-        *closestPnt = point.XYZ();
+    {
-      }
+      gp_Vec edge( xyz[ pos._index ], xyz[ pos._index+1 ]);
-      else {
+      gp_Vec n1p ( xyz[ pos._index ], point  );
-        tmpPnt.SetY( 0 );
+      double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
-        trsf.Inverted().Transforms( tmpPnt );
+      // projection of the point on the edge
-        *closestPnt = tmpPnt;
+      gp_XYZ proj = xyz[ pos._index ] + u * edge.XYZ();
+      double dist2 = point.SquareDistance( proj );
+      if ( dist2 < minDist2 )
+      {
+        if ( closestPnt ) *closestPnt = proj;
+        minDist2 = dist2;
       }
+      break;
+    }
+
+    case POS_RIGHT: // point is inside the face
+    {
+      double distToFacePlane = Abs( tmpPnt.Y() );
+      if ( closestPnt )
+      {
+        if ( distToFacePlane < std::numeric_limits<double>::min() ) {
+          *closestPnt = point.XYZ();
+        }
+        else {
+          tmpPnt.SetY( 0 );
+          trsf.Inverted().Transforms( tmpPnt );
+          *closestPnt = tmpPnt;
+        }
+      }
+      return distToFacePlane;
+    }
+
+    case POS_VERTEX: // point is most close to a node
+    {
+      double dist2 = point.SquareDistance( xyz[ pos._index ]);
+      if ( dist2 < minDist2 )
+      {
+        if ( closestPnt ) *closestPnt = xyz[ pos._index ];
+        minDist2 = dist2;
+      }
+      break;
+    }
+    default:;
+      return badDistance;
     }
-    return distToFacePlane;
   }
-  case POS_VERTEX:
+  return Sqrt( minDist2 );
-  {
-    // point is most close to a node
-    gp_Vec distVec( point, xyz[ pos._index ]);
-    return distVec.Magnitude();
-  }
-  default:;
-  }
-  return badDistance;
 }
 
 //=======================================================================

src/SMESH_I/SMESH_MeshEditor_i.cxx

@@ -7224,15 +7224,15 @@ void SMESH_MeshEditor_i::MakePolyLine(SMESH::ListOfPolySegments& theSegments,
     segOut.myNode2[0] = meshDS->FindNode( segIn.node1ID2 );
     segOut.myNode1[1] = meshDS->FindNode( segIn.node2ID1 );
     segOut.myNode2[1] = meshDS->FindNode( segIn.node2ID2 );
+    segOut.myXYZ[0].SetCoord( segIn.xyz1.x,
+                              segIn.xyz1.y,
+                              segIn.xyz1.z);
+    segOut.myXYZ[1].SetCoord( segIn.xyz2.x,
+                              segIn.xyz2.y,
+                              segIn.xyz2.z);
     segOut.myVector.SetCoord( segIn.vector.PS.x,
                               segIn.vector.PS.y,
                               segIn.vector.PS.z );
-    if ( !segOut.myNode1[0] )
-      THROW_SALOME_CORBA_EXCEPTION( SMESH_Comment( "Invalid node ID: ") << segIn.node1ID1,
-                                    SALOME::BAD_PARAM );
-    if ( !segOut.myNode1[1] )
-      THROW_SALOME_CORBA_EXCEPTION( SMESH_Comment( "Invalid node ID: ") << segIn.node2ID1,
-                                    SALOME::BAD_PARAM );
   }
 
   // get a static ElementSearcher

src/SMESH_SWIG/smeshBuilder.py

@@ -4536,7 +4536,7 @@ class Mesh(metaclass = MeshMeta):
 
         Parameters:
                 IDsOfElements: the faces to be splitted
-                Diag13:        is used to choose a diagonal for splitting.
+                Diag13 (boolean):        is used to choose a diagonal for splitting.
 
         Returns:
             True in case of success, False otherwise.
@@ -4552,7 +4552,7 @@ class Mesh(metaclass = MeshMeta):
         Parameters:
                 theObject: the object from which the list of elements is taken,
                         this is :class:`mesh, sub-mesh, group or filter <SMESH.SMESH_IDSource>`
-                Diag13:    is used to choose a diagonal for splitting.
+                Diag13 (boolean):    is used to choose a diagonal for splitting.
 
         Returns:
             True in case of success, False otherwise.
@@ -6665,7 +6665,7 @@ class Mesh(metaclass = MeshMeta):
     def MakePolyLine(self, segments, groupName='', isPreview=False ):
         """    
         Create a polyline consisting of 1D mesh elements each lying on a 2D element of
-        the initial mesh. Positions of new nodes are found by cutting the mesh by the
+        the initial triangle mesh. Positions of new nodes are found by cutting the mesh by the
         plane passing through pairs of points specified by each :class:`SMESH.PolySegment` structure.
         If there are several paths connecting a pair of points, the shortest path is
         selected by the module. Position of the cutting plane is defined by the two
@@ -6675,7 +6675,7 @@ class Mesh(metaclass = MeshMeta):
         The vector goes from the middle point to the projection point. In case of planar
         mesh, the vector is normal to the mesh.
 
-        *segments* [i].vector returns the used vector which goes from the middle point to its projection.
+        In preview mode, *segments* [i].vector returns the used vector which goes from the middle point to its projection.
 
         Parameters:
             segments: list of :class:`SMESH.PolySegment` defining positions of cutting planes.