From ce56cca5177c94ff857383602c3c308565e9ece3 Mon Sep 17 00:00:00 2001 From: eap Date: Thu, 18 Apr 2019 21:34:39 +0300 Subject: [PATCH] #16843 EDF 19340 - wrong quadratic mesh --- src/SMESH/SMESH_Algo.cxx | 30 +++++ src/SMESH/SMESH_Algo.hxx | 24 ++++ src/SMESH/SMESH_MesherHelper.cxx | 77 ++++++------- src/SMESH_SWIG/smeshBuilder.py | 4 +- .../StdMeshers_QuadToTriaAdaptor.cxx | 107 +++++++++++++----- 5 files changed, 172 insertions(+), 70 deletions(-) diff --git a/src/SMESH/SMESH_Algo.cxx b/src/SMESH/SMESH_Algo.cxx index 4de1879f1..12371ea64 100644 --- a/src/SMESH/SMESH_Algo.cxx +++ b/src/SMESH/SMESH_Algo.cxx @@ -1268,3 +1268,33 @@ bool SMESH_2D_Algo::FixInternalNodes(const SMESH_ProxyMesh& mesh, } return true; } + +//======================================================================= +//function : IsApplicableToShape +//purpose : Return true if the algorithm can mesh a given shape +//======================================================================= + +bool SMESH_1D_Algo::IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const +{ + return ( !shape.IsNull() && TopExp_Explorer( shape, TopAbs_EDGE ).More() ); +} + +//======================================================================= +//function : IsApplicableToShape +//purpose : Return true if the algorithm can mesh a given shape +//======================================================================= + +bool SMESH_2D_Algo::IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const +{ + return ( !shape.IsNull() && TopExp_Explorer( shape, TopAbs_FACE ).More() ); +} + +//======================================================================= +//function : IsApplicableToShape +//purpose : Return true if the algorithm can mesh a given shape +//======================================================================= + +bool SMESH_3D_Algo::IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const +{ + return ( !shape.IsNull() && TopExp_Explorer( shape, TopAbs_SOLID ).More() ); +} diff --git a/src/SMESH/SMESH_Algo.hxx b/src/SMESH/SMESH_Algo.hxx index 2cc313af4..e35a47120 100644 --- a/src/SMESH/SMESH_Algo.hxx +++ b/src/SMESH/SMESH_Algo.hxx @@ -482,12 +482,28 @@ class SMESH_EXPORT SMESH_1D_Algo: public SMESH_Algo { public: SMESH_1D_Algo(int hypId, SMESH_Gen* gen); + /*! + * \brief Return true if the algorithm can mesh a given shape + * \param [in] aShape - shape to check + * \param [in] toCheckAll - if true, this check returns OK if all shapes are OK, + * else, returns OK if at least one shape is OK + * \retval bool - \c true by default + */ + virtual bool IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const; }; class SMESH_EXPORT SMESH_2D_Algo: public SMESH_Algo { public: SMESH_2D_Algo(int hypId, SMESH_Gen* gen); + /*! + * \brief Return true if the algorithm can mesh a given shape + * \param [in] aShape - shape to check + * \param [in] toCheckAll - if true, this check returns OK if all shapes are OK, + * else, returns OK if at least one shape is OK + * \retval bool - \c true by default + */ + virtual bool IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const; /*! * \brief Method in which an algorithm generating a structured mesh * fixes positions of in-face nodes after there movement @@ -501,6 +517,14 @@ class SMESH_EXPORT SMESH_3D_Algo: public SMESH_Algo { public: SMESH_3D_Algo(int hypId, SMESH_Gen* gen); + /*! + * \brief Return true if the algorithm can mesh a given shape + * \param [in] aShape - shape to check + * \param [in] toCheckAll - if true, this check returns OK if all shapes are OK, + * else, returns OK if at least one shape is OK + * \retval bool - \c true by default + */ + virtual bool IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const; }; #endif diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index eceb92e36..390278faa 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -1673,30 +1673,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, // get positions of the given nodes on shapes if ( pos.second == TopAbs_FACE ) { - F = TopoDS::Face(meshDS->IndexToShape( faceID = pos.first )); + F = TopoDS::Face( meshDS->IndexToShape( faceID = pos.first )); uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]); - if (( !force3d ) && - ( HasDegeneratedEdges() || GetSurface( F )->HasSingularities( 1e-7 ))) - { - // IPAL52850 (degen VERTEX not at singularity) - // project middle point to a surface - SMESH_TNodeXYZ p1( n1 ), p2( n2 ); - gp_Pnt pMid = 0.5 * ( p1 + p2 ); - Handle(ShapeAnalysis_Surface) projector = GetSurface( F ); - gp_Pnt2d uvMid; - if ( uvOK[0] ) - uvMid = projector->NextValueOfUV( uv[0], pMid, BRep_Tool::Tolerance( F )); - else - uvMid = projector->ValueOfUV( pMid, getFaceMaxTol( F )); - if ( projector->Gap() * projector->Gap() < ( p1 - p2 ).SquareModulus() / 4 ) - { - gp_Pnt pProj = projector->Value( uvMid ); - n12 = meshDS->AddNode( pProj.X(), pProj.Y(), pProj.Z() ); - meshDS->SetNodeOnFace( n12, faceID, uvMid.X(), uvMid.Y() ); - myTLinkNodeMap.insert( make_pair ( link, n12 )); - return n12; - } - } uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]); } else if ( pos.second == TopAbs_EDGE ) @@ -1729,26 +1707,43 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, // nodes, else - medium between corresponding 3d points if( ! F.IsNull() ) { - //if ( uvOK[0] && uvOK[1] ) + if ( IsDegenShape( n1->getshapeId() )) { + if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 )); + else uv[0].SetCoord( 2, uv[1].Coord( 2 )); + } + else if ( IsDegenShape( n2->getshapeId() )) { + if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 )); + else uv[1].SetCoord( 2, uv[0].Coord( 2 )); + } + TopLoc_Location loc; + Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc); + gp_XY UV = GetMiddleUV( S, uv[0], uv[1] ); + gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc); + + SMESH_TNodeXYZ p1( n1 ), p2( n2 ); + gp_Pnt pMid = 0.5 * ( p1 + p2 ); + double distMid = pMid.SquareDistance( P ); + double dist12 = ( p1 - p2 ).SquareModulus(); + Handle(ShapeAnalysis_Surface) surfInfo = GetSurface( F ); + if ( distMid > dist12 || + HasDegeneratedEdges() || + surfInfo->HasSingularities( 1e-7 ) ) { - if ( IsDegenShape( n1->getshapeId() )) { - if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 )); - else uv[0].SetCoord( 2, uv[1].Coord( 2 )); - } - else if ( IsDegenShape( n2->getshapeId() )) { - if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 )); - else uv[1].SetCoord( 2, uv[0].Coord( 2 )); - } - TopLoc_Location loc; - Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc); - gp_XY UV = GetMiddleUV( S, uv[0], uv[1] ); - gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc); - n12 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - // if ( mySetElemOnShape ) node is not elem! - meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y()); - myTLinkNodeMap.insert(make_pair(link,n12)); - return n12; + // IPAL52850 (degen VERTEX not at singularity) + // project middle point to a surface + gp_Pnt2d uvMid; + if ( uvOK[0] ) + uvMid = surfInfo->NextValueOfUV( uv[0], pMid, BRep_Tool::Tolerance( F )); + else + uvMid = surfInfo->ValueOfUV( pMid, getFaceMaxTol( F )); + if ( surfInfo->Gap() * surfInfo->Gap() < distMid ) + P = surfInfo->Value( uvMid ); } + n12 = meshDS->AddNode(P.X(), P.Y(), P.Z()); + // if ( mySetElemOnShape ) node is not elem! + meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y()); + myTLinkNodeMap.insert(make_pair(link,n12)); + return n12; } else if ( !E.IsNull() ) { diff --git a/src/SMESH_SWIG/smeshBuilder.py b/src/SMESH_SWIG/smeshBuilder.py index cce91ec02..26332db27 100644 --- a/src/SMESH_SWIG/smeshBuilder.py +++ b/src/SMESH_SWIG/smeshBuilder.py @@ -4490,8 +4490,8 @@ class Mesh(metaclass = MeshMeta): Reorient faces contained in *the2DObject*. Parameters: - the2DObject: is a :class:`mesh, sub-mesh, group or filter ` or list of IDs of 2D elements - theDirection: is a desired direction of normal of *theFace*. + the2DObject: a :class:`mesh, sub-mesh, group or filter ` or list of IDs of 2D elements + theDirection: a desired direction of normal of *theFace*. It can be either a GEOM vector or a list of coordinates [x,y,z]. theFaceOrPoint: defines a face of *the2DObject* whose normal will be compared with theDirection. It can be either ID of face or a point diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx index 993956225..f0d50dc0f 100644 --- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx +++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx @@ -258,6 +258,7 @@ namespace meshDS->RemoveFreeNode( nodesToRemove[i], sm, /*fromGroups=*/false); } } + return; } //================================================================================ @@ -319,7 +320,7 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr set & nodesToMove) { // cout << endl << "Merge " << PrmI->GetID() << " " << PrmJ->GetID() << " " - // << PrmI->GetNode(4) << PrmJ->GetNode(4) << endl; + // << PrmI->GetNode(4)->GetID() << " " << PrmJ->GetNode(4)->GetID() << endl; const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume ); SMESH_TNodeXYZ Pj( Nrem ); @@ -338,6 +339,9 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator; TStdElemIterator itEnd; + typedef std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > TNNMap; + TNNMap mediumReplaceMap; + // find and remove coincided faces of merged pyramids vector< const SMDS_MeshElement* > inverseElems // copy inverse elements to avoid iteration on changing container @@ -355,6 +359,11 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr } if ( FJEqual ) { + if ( FJEqual->NbNodes() == 6 ) // find medium nodes to replace + { + mediumReplaceMap.insert( std::make_pair( FJEqual->GetNode(3), FI->GetNode(5) )); + mediumReplaceMap.insert( std::make_pair( FJEqual->GetNode(5), FI->GetNode(3) )); + } removeTmpElement( FI ); removeTmpElement( FJEqual ); myRemovedTrias.insert( FI ); @@ -370,12 +379,27 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr const SMDS_MeshElement* elem = inverseElems[i]; nodes.assign( elem->begin_nodes(), elem->end_nodes() ); nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode; + if ( !mediumReplaceMap.empty() ) + for ( size_t iN = elem->NbCornerNodes(); iN < nodes.size(); ++iN ) + { + TNNMap::iterator n2n = mediumReplaceMap.find( nodes[iN] ); + if ( n2n != mediumReplaceMap.end() ) + nodes[iN] = n2n->second; + } GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size()); } ASSERT( Nrem->NbInverseElements() == 0 ); GetMeshDS()->RemoveFreeNode( Nrem, GetMeshDS()->MeshElements( Nrem->getshapeId()), /*fromGroups=*/false); + if ( !mediumReplaceMap.empty() ) + for ( TNNMap::iterator n2n = mediumReplaceMap.begin(); n2n != mediumReplaceMap.end(); ++n2n ) + { + const SMDS_MeshNode* remNode = n2n->first; + if ( !remNode->IsNull() && remNode->NbInverseElements() == 0 ) + GetMeshDS()->RemoveFreeNode( remNode, 0, /*fromGroups=*/false); + } + return; } //================================================================================ @@ -390,7 +414,7 @@ void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement* PrmI { TIDSortedElemSet adjacentPyrams; bool mergedPyrams = false; - for ( int k=0; k<4; k++ ) // loop on 4 base nodes of PrmI + for ( int k = 0; k < 4; k++ ) // loop on 4 base nodes of PrmI { const SMDS_MeshNode* n = PrmI->GetNode(k); SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); @@ -414,6 +438,7 @@ void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement* PrmI for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm) MergeAdjacent( *prm, nodesToMove, true ); } + return; } //================================================================================ @@ -474,7 +499,7 @@ static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2, // and a segment [PC,P] //======================================================================= -static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, +static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3) { const double EPSILON = 1e-6; @@ -522,7 +547,27 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, Pint = orig + dir * t; - return ( t > 0. && t < segLen ); + bool hasInt = ( t > 0. && t < segLen ); + + if ( hasInt && det < EPSILON ) // t is inaccurate, additionally check + { + gp_XYZ triNorm = edge1 ^ edge2; + gp_XYZ int0vec = Pint.XYZ() - vert0; + gp_XYZ in = triNorm ^ edge1; // dir inside triangle from edge1 + double dot = int0vec * in; + if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON ) + return false; + in = edge2 ^ triNorm; + dot = int0vec * in; + if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON ) + return false; + gp_XYZ int1vec = Pint.XYZ() - vert1; + in = triNorm ^ ( vert2 - vert1 ); + dot = int1vec * in; + if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON ) + return false; + } + return hasInt; } //======================================================================= @@ -817,7 +862,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); SMESH_MesherHelper helper(aMesh); helper.IsQuadraticSubMesh(aShape); - helper.SetElementsOnShape( true ); if ( myElemSearcher ) delete myElemSearcher; vector< SMDS_ElemIteratorPtr > itVec; @@ -877,10 +921,11 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, // degenerate face // add triangles to result map SMDS_MeshFace* NewFace; + helper.SetElementsOnShape( false ); if(!isRev) - NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] ); + NewFace = helper.AddFace( FNodes[0], FNodes[1], FNodes[2] ); else - NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] ); + NewFace = helper.AddFace( FNodes[0], FNodes[2], FNodes[1] ); storeTmpElement( NewFace ); trias.push_back ( NewFace ); quads.push_back( face ); @@ -916,20 +961,25 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, if ( !LimitHeight( PCbest, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/true, aShape )) return false; } - // create node for PCbest + // create node at PCbest + helper.SetElementsOnShape( true ); SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); + // create a pyramid + SMDS_MeshVolume* aPyram; + if ( isRev ) + aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode ); + else + aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); + myPyramids.push_back(aPyram); + // add triangles to result map - for(i=0; i<4; i++) + helper.SetElementsOnShape( false ); + for ( i = 0; i < 4; i++ ) { - trias.push_back ( meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] )); + trias.push_back ( helper.AddFace( NewNode, FNodes[i], FNodes[i+1] )); storeTmpElement( trias.back() ); } - // create a pyramid - if ( isRev ) swap( FNodes[1], FNodes[3]); - SMDS_MeshVolume* aPyram = - helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); - myPyramids.push_back(aPyram); quads.push_back( face ); hasNewTrias = true; @@ -1018,7 +1068,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) vector myPyramids; SMESH_MesherHelper helper(aMesh); helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh()); - helper.SetElementsOnShape( true ); SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh(); @@ -1116,10 +1165,11 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) IsRev = true; } } + helper.SetElementsOnShape( false ); if(!IsRev) - NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] ); + NewFace = helper.AddFace( FNodes[0], FNodes[1], FNodes[2] ); else - NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] ); + NewFace = helper.AddFace( FNodes[0], FNodes[2], FNodes[1] ); storeTmpElement( NewFace ); prxSubMesh->AddElement( NewFace ); continue; @@ -1238,25 +1288,28 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) return false; // create node for Papex + helper.SetElementsOnShape( true ); SMDS_MeshNode* NewNode = helper.AddNode( Papex.X(), Papex.Y(), Papex.Z() ); + // create a pyramid + SMDS_MeshVolume* aPyram; + if(isRev) + aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); + else + aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode ); + myPyramids.push_back(aPyram); + // add triangles to result map + helper.SetElementsOnShape( false ); for ( i = 0; i < 4; i++) { SMDS_MeshFace* NewFace; if(isRev) - NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ); + NewFace = helper.AddFace( NewNode, FNodes[i], FNodes[i+1] ); else - NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] ); + NewFace = helper.AddFace( NewNode, FNodes[i+1], FNodes[i] ); storeTmpElement( NewFace ); prxSubMesh->AddElement( NewFace ); } - // create a pyramid - SMDS_MeshVolume* aPyram; - if(isRev) - aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); - else - aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode ); - myPyramids.push_back(aPyram); } } // end loop on all faces -- 2.39.2