From: eap Date: Fri, 9 Nov 2018 10:28:18 +0000 (+0300) Subject: PAL0023627: [IMACS] ASERIS: project point to the mesh X-Git-Tag: V9_3_0a1~49 X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=commitdiff_plain;h=441a2df90cb97ea7d035771a49e28dc53469e3d6 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. --- diff --git a/doc/salome/gui/SMESH/input/smesh_module.rst b/doc/salome/gui/SMESH/input/smesh_module.rst index 8c2fda6e6..1178488c8 100644 --- a/doc/salome/gui/SMESH/input/smesh_module.rst +++ b/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 diff --git a/idl/SMESH_MeshEditor.idl b/idl/SMESH_MeshEditor.idl index 2ec0484ea..5cb9c85cc 100644 --- a/idl/SMESH_MeshEditor.idl +++ b/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 - // by two nodes, else it is at node1ID1 + // a point is defined as follows: + // ( 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 - // by two nodes, else it is at node2ID1 + // point 2 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 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 - * plane passing through pairs of points specified by each PolySegment structure. + * 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 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. diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx index ceaeb1acb..cea314220 100644 --- a/src/SMESH/SMESH_MeshEditor.cxx +++ b/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] ); - 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 + // 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 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::min() ); diff --git a/src/SMESH/SMESH_MeshEditor.hxx b/src/SMESH/SMESH_MeshEditor.hxx index b6c404a15..8e605ab92 100644 --- a/src/SMESH/SMESH_MeshEditor.hxx +++ b/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 - // by two nodes, else it is at myNode1 - const SMDS_MeshNode* myNode1[2]; - const SMDS_MeshNode* myNode2[2]; - - gp_Vec myVector; // vector on the plane; to use a default plane set vector = (0,0,0) + // 2 points, each defined as follows: + // ( myNode1 && myNode2 ) ==> point is in the middle of an edge defined by two nodes + // ( myNode1 && !myNode2 ) ==> point is at myNode1 + // else ==> point is at myXYZ + const SMDS_MeshNode* myNode1[2]; + const SMDS_MeshNode* myNode2[2]; + gp_XYZ myXYZ [2]; + + // 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; diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index 0fc388e75..f42e24c11 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.cxx +++ b/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 ) - { - case POS_LEFT: - { - // point is most close to an edge - gp_Vec edge( xyz[ pos._index ], xyz[ pos._index+1 ]); - 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: + + double minDist2 = Precision::Infinite(); + for ( std::set< PointPos >::iterator posIt = pntPosSet.begin(); posIt != pntPosSet.end(); ++posIt) { - // point is inside the face - double distToFacePlane = Abs( tmpPnt.Y() ); - if ( closestPnt ) + PointPos pos = *posIt; + if ( pos._name != pntPosSet.begin()->_name ) + break; + switch ( pos._name ) + { + case POS_LEFT: // point is most close to an edge + { + gp_Vec edge( xyz[ pos._index ], xyz[ pos._index+1 ]); + 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(); + double dist2 = point.SquareDistance( proj ); + if ( dist2 < minDist2 ) + { + if ( closestPnt ) *closestPnt = proj; + minDist2 = dist2; + } + break; + } + + case POS_RIGHT: // point is inside the face { - if ( distToFacePlane < std::numeric_limits::min() ) { - *closestPnt = point.XYZ(); + double distToFacePlane = Abs( tmpPnt.Y() ); + if ( closestPnt ) + { + if ( distToFacePlane < std::numeric_limits::min() ) { + *closestPnt = point.XYZ(); + } + else { + tmpPnt.SetY( 0 ); + trsf.Inverted().Transforms( tmpPnt ); + *closestPnt = tmpPnt; + } } - 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: - { - // point is most close to a node - gp_Vec distVec( point, xyz[ pos._index ]); - return distVec.Magnitude(); - } - default:; } - return badDistance; + return Sqrt( minDist2 ); } //======================================================================= diff --git a/src/SMESH_I/SMESH_MeshEditor_i.cxx b/src/SMESH_I/SMESH_MeshEditor_i.cxx index 343831685..9ca11239b 100644 --- a/src/SMESH_I/SMESH_MeshEditor_i.cxx +++ b/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 diff --git a/src/SMESH_SWIG/smeshBuilder.py b/src/SMESH_SWIG/smeshBuilder.py index fed6214bb..e1af32d91 100755 --- a/src/SMESH_SWIG/smeshBuilder.py +++ b/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 ` - 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,12 +6675,12 @@ 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: + Parameters: segments: list of :class:`SMESH.PolySegment` defining positions of cutting planes. groupName: optional name of a group where created mesh segments will be added. - + """ editor = self.editor if isPreview: