X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=914eebff3837e353de40c316cdce4ecd0c76acdb;hp=17702dd30d6a0e1f81589de66ac1692a52734594;hb=refs%2Ftags%2FV9_7_0b1;hpb=6d32f944a0a115b6419184c50b57bf7c4eef5786 diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index 17702dd30..914eebff3 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2019 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -51,6 +51,7 @@ #include #include #include +#include #include #include #include @@ -203,12 +204,12 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis //============================================================================= /*! - * + * Compute the mesh on the given shape */ //============================================================================= -bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape) +bool StdMeshers_Quadrangle_2D::Compute( SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape ) { const TopoDS_Face& F = TopoDS::Face(aShape); aMesh.GetSubMesh( F ); @@ -916,7 +917,7 @@ bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh, std::vector aNbNodes(4); bool IsQuadratic = false; if (!checkNbEdgesForEvaluate(aMesh, aFace, aResMap, aNbNodes, IsQuadratic)) { - std::vector aResVec(SMDSEntity_Last); + std::vector aResVec(SMDSEntity_Last); for (int i=SMDSEntity_Node; i aVec(SMDSEntity_Last,0); + std::vector aVec(SMDSEntity_Last,0); if (IsQuadratic) { aVec[SMDSEntity_Quad_Triangle] = nbFaces3; aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4; @@ -1048,6 +1049,74 @@ namespace return true; } + //================================================================================ + /*! + * \brief Return angle between mesh segments of given EDGEs meeting at theVertexNode + */ + //================================================================================ + + double getAngleByNodes( const int theE1Index, + const int theE2Index, + const SMDS_MeshNode* theVertexNode, + const StdMeshers_FaceSide& theFaceSide, + const gp_Vec& theFaceNormal) + { + int eID1 = theFaceSide.EdgeID( theE1Index ); + int eID2 = theFaceSide.EdgeID( theE2Index ); + + const SMDS_MeshNode *n1 = 0, *n2 = 0; + bool is1st; + SMDS_ElemIteratorPtr segIt = theVertexNode->GetInverseElementIterator( SMDSAbs_Edge ); + while ( segIt->more() ) + { + const SMDS_MeshElement* seg = segIt->next(); + int shapeID = seg->GetShapeID(); + if ( shapeID == eID1 ) + is1st = true; + else if ( shapeID == eID2 ) + is1st = false; + else + continue; + ( is1st ? n1 : n2 ) = seg->GetNodeWrap( 1 + seg->GetNodeIndex( theVertexNode )); + } + + if ( !n1 || !n2 ) + { + std::vector nodes; + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + const SMDS_MeshNode* & n = is2nd ? n2 : n1; + if ( n ) continue; + nodes.clear(); + if ( is2nd ) theFaceSide.GetEdgeNodes( theE2Index, nodes ); + else theFaceSide.GetEdgeNodes( theE1Index, nodes ); + if ( nodes.size() >= 2 ) + { + if ( nodes[0] == theVertexNode ) + n = nodes[1]; + else + n = nodes[ nodes.size() - 2 ]; + } + } + } + double angle = -2 * M_PI; + if ( n1 && n2 ) + { + SMESH_NodeXYZ p1 = n1, p2 = theVertexNode, p3 = n2; + gp_Vec v1( p1, p2 ), v2( p2, p3 ); + try + { + angle = v1.AngleWithRef( v2, theFaceNormal ); + } + catch(...) + { + } + if ( std::isnan( angle )) + angle = -2 * M_PI; + } + return angle; + } + //-------------------------------------------------------------------------------- /*! * \brief EDGE of a FACE @@ -1057,6 +1126,7 @@ namespace TopoDS_Edge myEdge; TopoDS_Vertex my1stVertex; int myIndex; + bool myIsCorner; // is fixed corner double myAngle; // angle at my1stVertex int myNbSegments; // discretization Edge* myPrev; // preceding EDGE @@ -1085,6 +1155,7 @@ namespace // quality criteria to minimize int myOppDiff; + int myIsFixedCorner; double myQuartDiff; double mySumAngle; @@ -1112,6 +1183,11 @@ namespace myOppDiff = ( Abs( myNbSeg[0] - myNbSeg[2] ) + Abs( myNbSeg[1] - myNbSeg[3] )); + myIsFixedCorner = - totNbSeg * ( myCornerE[0]->myIsCorner + + myCornerE[1]->myIsCorner + + myCornerE[2]->myIsCorner + + myCornerE[3]->myIsCorner ); + double nbSideIdeal = totNbSeg / 4.; myQuartDiff = -( Min( Min( myNbSeg[0], myNbSeg[1] ), Min( myNbSeg[2], myNbSeg[3] )) / nbSideIdeal ); @@ -1125,7 +1201,7 @@ namespace }; // first criterion - equality of nbSeg of opposite sides - int crit1() const { return myOppDiff; } + int crit1() const { return myOppDiff + myIsFixedCorner; } // second criterion - equality of nbSeg of adjacent sides and sharpness of angles double crit2() const { return myQuartDiff + mySumAngle; } @@ -1143,17 +1219,19 @@ namespace //================================================================================ /*! * \brief Unite EDGEs to get a required number of sides - * \param [in] theNbCorners - the required number of sides + * \param [in] theNbCorners - the required number of sides, 3 or 4 * \param [in] theConsiderMesh - to considered only meshed VERTEXes * \param [in] theFaceSide - the FACE EDGEs + * \param [in] theFixedVertices - VERTEXes to be used as corners * \param [out] theVertices - the found corner vertices + * \param [out] theHaveConcaveVertices - return if there are concave vertices */ //================================================================================ void uniteEdges( const int theNbCorners, const bool theConsiderMesh, const StdMeshers_FaceSide& theFaceSide, - const TopoDS_Shape& theBaseVertex, + const TopTools_MapOfShape& theFixedVertices, std::vector& theVertices, bool& theHaveConcaveVertices) { @@ -1183,23 +1261,28 @@ namespace // sort edges by angle std::multimap< double, Edge* > edgeByAngle; - int i, iBase = -1, nbConvexAngles = 0, nbSharpAngles = 0; + int i, nbConvexAngles = 0, nbSharpAngles = 0; + const SMDS_MeshNode* vertNode = 0; + gp_Vec faceNormal; const double angTol = 5. / 180 * M_PI; const double sharpAngle = 0.5 * M_PI - angTol; Edge* e = edge0; for ( i = 0; i < nbEdges; ++i, e = e->myNext ) { e->my1stVertex = SMESH_MesherHelper::IthVertex( 0, e->myEdge ); - if ( e->my1stVertex.IsSame( theBaseVertex )) - iBase = e->myIndex; + e->myIsCorner = theFixedVertices.Contains( e->my1stVertex ); e->myAngle = -2 * M_PI; - if ( !theConsiderMesh || theFaceSide.VertexNode( e->myIndex )) + if ( !theConsiderMesh || ( vertNode = theFaceSide.VertexNode( e->myIndex ))) { e->myAngle = SMESH_MesherHelper::GetAngle( e->myPrev->myEdge, e->myEdge, - theFaceSide.Face(), e->my1stVertex ); + theFaceSide.Face(), e->my1stVertex, + &faceNormal ); if ( e->myAngle > 2 * M_PI ) // GetAngle() failed e->myAngle *= -1.; + else if ( vertNode && ( 0. <= e->myAngle ) && ( e->myAngle <= angTol )) + e->myAngle = getAngleByNodes( e->myPrev->myIndex, e->myIndex, + vertNode, theFaceSide, faceNormal ); } edgeByAngle.insert( std::make_pair( e->myAngle, e )); nbConvexAngles += ( e->myAngle > angTol ); @@ -1227,8 +1310,10 @@ namespace // return corners with maximal angles std::set< int > cornerIndices; - if ( iBase != -1 ) - cornerIndices.insert( iBase ); + if ( !theFixedVertices.IsEmpty() ) + for ( i = 0, e = edge0; i < nbEdges; ++i, e = e->myNext ) + if ( e->myIsCorner ) + cornerIndices.insert( e->myIndex ); std::multimap< double, Edge* >::reverse_iterator a2e = edgeByAngle.rbegin(); for (; (int) cornerIndices.size() < theNbCorners; ++a2e ) @@ -1350,6 +1435,39 @@ namespace return; } + //================================================================================ + /*! + * \brief Remove a seam and degenerated edge from a wire if the shape is + * a quadrangle with a seam inside. + */ + //================================================================================ + + bool removeInternalSeam( std::list& theWire, + SMESH_MesherHelper& theHelper) + { + if ( !theHelper.HasRealSeam() || + theHelper.NbDegeneratedEdges() != 2 ) // 1 EDGE + 1 VERTEX + return false; + + typedef std::list::iterator TEdgeIter; + std::vector< TEdgeIter > edgesToRemove; + edgesToRemove.reserve( 5 ); + for ( TEdgeIter eIt = theWire.begin(); eIt != theWire.end(); ++eIt ) + { + int eID = theHelper.ShapeToIndex( *eIt ); + if ( theHelper.IsRealSeam( eID ) || theHelper.IsDegenShape( eID )) + edgesToRemove.push_back( eIt ); + } + + if ( theWire.size() - edgesToRemove.size() < 4 ) + return false; // cone e.g. + + for ( size_t i = 0; i < edgesToRemove.size(); ++i ) + theWire.erase( edgesToRemove[ i ]); + + return true; + } + } // namespace //================================================================================ @@ -1380,6 +1498,9 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, if ( myHelper ) helper.CopySubShapeInfo( *myHelper ); + if ( removeInternalSeam( theWire, helper )) + theNbDegenEdges = 1; + StdMeshers_FaceSide faceSide( theFace, theWire, &theMesh, /*isFwd=*/true, /*skipMedium=*/true, &helper ); @@ -1410,7 +1531,7 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, if ( theConsiderMesh ) { - const int nbSegments = Max( faceSide.NbPoints()-1, faceSide.NbSegments() ); + const smIdType nbSegments = std::max( faceSide.NbPoints()-1, faceSide.NbSegments() ); if ( nbSegments < nbCorners ) return error(COMPERR_BAD_INPUT_MESH, TComm("Too few boundary nodes: ") << nbSegments); } @@ -1444,7 +1565,20 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, myCheckOri = false; if ( theVertices.size() > 3 ) { - uniteEdges( nbCorners, theConsiderMesh, faceSide, triaVertex, theVertices, myCheckOri ); + TopTools_MapOfShape fixedVertices; + if ( !triaVertex.IsNull() ) + fixedVertices.Add( triaVertex ); + if ( myParams ) + { + const std::vector< int >& vIDs = myParams->GetCorners(); + for ( size_t i = 0; i < vIDs.size(); ++i ) + { + const TopoDS_Shape& vertex = helper.GetMeshDS()->IndexToShape( vIDs[ i ]); + if ( !vertex.IsNull() ) + fixedVertices.Add( vertex ); + } + } + uniteEdges( nbCorners, theConsiderMesh, faceSide, fixedVertices, theVertices, myCheckOri ); } if ( nbCorners == 3 && !triaVertex.IsSame( theVertices[0] )) @@ -1472,7 +1606,8 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, //============================================================================= /*! - * + * Return FaceQuadStruct where sides ordered CCW, top and left sides + * reversed to be co-directed with bottom and right sides */ //============================================================================= @@ -1633,7 +1768,7 @@ bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh& aMes if (anIt==aResMap.end()) { return false; } - std::vector aVec = (*anIt).second; + std::vector aVec = (*anIt).second; IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]); if (nbEdgesInWire.front() == 3) { // exactly 3 edges if (myTriaVertexID>0) { @@ -1655,7 +1790,7 @@ bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh& aMes SMESH_subMesh * sm = aMesh.GetSubMesh(E1); MapShapeNbElemsItr anIt = aResMap.find(sm); if (anIt==aResMap.end()) return false; - std::vector aVec = (*anIt).second; + std::vector aVec = (*anIt).second; if (IsQuadratic) aNbNodes[0] = (aVec[SMDSEntity_Node]-1)/2 + 2; else @@ -1689,7 +1824,7 @@ bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh& aMes if (anIt==aResMap.end()) { return false; } - std::vector aVec = (*anIt).second; + std::vector aVec = (*anIt).second; if (IsQuadratic) aNbNodes[nbSides] = (aVec[SMDSEntity_Node]-1)/2 + 2; else @@ -1726,7 +1861,7 @@ bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh& aMes if (anIt==aResMap.end()) { return false; } - std::vector aVec = (*anIt).second; + std::vector aVec = (*anIt).second; if (IsQuadratic) aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1; else @@ -1767,7 +1902,7 @@ bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh& aMes if (anIt==aResMap.end()) { return false; } - std::vector aVec = (*anIt).second; + std::vector aVec = (*anIt).second; if (IsQuadratic) aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1; else @@ -3044,7 +3179,7 @@ bool StdMeshers_Quadrangle_2D::evaluateQuadPref(SMESH_Mesh & aMesh, nbFaces += (drl+addv)*(nb-1) + (nt-1); } // end new version implementation - std::vector aVec(SMDSEntity_Last); + std::vector aVec(SMDSEntity_Last); for (int i=SMDSEntity_Node; iUVPt( nbhoriz-1, nbvertic-1 ).node ); SMESH_TNodeXYZ a3( quad->UVPt( 0, nbvertic-1 ).node ); + // compute TFI for (int i = 1; i < nbhoriz-1; i++) { SMESH_TNodeXYZ p0( quad->UVPt( i, 0 ).node ); @@ -4347,13 +4483,39 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad) UVPtStruct& uvp = quad->UVPt( i, j ); - gp_Pnt p = myHelper->calcTFI(uvp.x,uvp.y, a0,a1,a2,a3, p0,p1,p2,p3); + gp_Pnt pnew = myHelper->calcTFI(uvp.x,uvp.y, a0,a1,a2,a3, p0,p1,p2,p3); + meshDS->MoveNode( uvp.node, pnew.X(), pnew.Y(), pnew.Z() ); + } + } + // project to surface + double cellSize; + for (int i = 1; i < nbhoriz-1; i++) + { + for (int j = 1; j < nbvertic-1; j++) + { + UVPtStruct& uvp = quad->UVPt( i, j ); + SMESH_NodeXYZ p = uvp.node; + + cellSize = Max( p.SquareDistance( quad->UVPt( i+1, j ).node ), + p.SquareDistance( quad->UVPt( i-1, j ).node )); + cellSize = Max( p.SquareDistance( quad->UVPt( i, j+1 ).node ), cellSize ); + cellSize = Max( p.SquareDistance( quad->UVPt( i, j-1 ).node ), cellSize ); + gp_Pnt2d uv = surface->NextValueOfUV( uvp.UV(), p, 10*tol ); gp_Pnt pnew = surface->Value( uv ); - - meshDS->MoveNode( uvp.node, pnew.X(), pnew.Y(), pnew.Z() ); - uvp.u = uv.X(); - uvp.v = uv.Y(); + bool ok = ( pnew.SquareDistance( p ) < 2 * cellSize ); + if ( !ok ) + { + uv = surface->ValueOfUV( p, 10*tol ); + pnew = surface->Value( uv ); + ok = ( pnew.SquareDistance( p ) < 2 * cellSize ); + } + if ( ok ) + { + meshDS->MoveNode( uvp.node, pnew.X(), pnew.Y(), pnew.Z() ); + uvp.u = uv.X(); + uvp.v = uv.Y(); + } } } } @@ -4614,6 +4776,7 @@ bool StdMeshers_Quadrangle_2D::check() const SMDS_MeshNode* nInFace = 0; if ( myHelper->HasSeam() ) + { for ( int i = 0; i < nbN && !nInFace; ++i ) if ( !myHelper->IsSeamShape( nn[i]->getshapeId() )) { @@ -4622,6 +4785,33 @@ bool StdMeshers_Quadrangle_2D::check() if ( myHelper->IsOnSeam( uv )) nInFace = NULL; } + } + if ( myHelper->GetPeriodicIndex() && !nInFace ) + { + for ( int i = 0; i < nbN && !nInFace; ++i ) + if ( fSubMesh->Contains( nn[i] )) + nInFace = nn[i]; + if ( !nInFace ) + for ( int i = 0; i < nbN && !nInFace; ++i ) + { + SMDS_ElemIteratorPtr fIt = nn[i]->GetInverseElementIterator( SMDSAbs_Face ); + while ( fIt->more() && !nInFace ) + { + const SMDS_MeshElement* face = fIt->next(); + if ( !fSubMesh->Contains( face )) + continue; + for ( int iN = 0, nN = face->NbCornerNodes(); iN < nN; ++iN ) + { + const SMDS_MeshNode* n = face->GetNode( iN ); + if ( fSubMesh->Contains( n )) + { + nInFace = n; + break; + } + } + } + } + } toCheckUV = true; for ( int i = 0; i < nbN; ++i ) @@ -4652,20 +4842,35 @@ bool StdMeshers_Quadrangle_2D::check() default:; } - // if ( isBad && myHelper->HasRealSeam() ) - // { - // // detect a case where a face intersects the seam - // for ( int iPar = 1; iPar < 3; ++iPar ) - // if ( iPar & myHelper->GetPeriodicIndex() ) - // { - // double min = uv[0].Coord( iPar ), max = uv[0].Coord( iPar ); - // for ( int i = 1; i < nbN; ++i ) - // { - // min = Min( min, uv[i].Coord( iPar )); - // max = Max( max, uv[i].Coord( iPar )); - // } - // } - // } + if ( isBad && myHelper->HasRealSeam() ) + { + // fix uv for a case where a face intersects the seam + for ( int iPar = 1; iPar < 3; ++iPar ) + if ( iPar & myHelper->GetPeriodicIndex() ) + { + double max = uv[0].Coord( iPar ); + for ( int i = 1; i < nbN; ++i ) + max = Max( max, uv[i].Coord( iPar )); + + for ( int i = 0; i < nbN; ++i ) + { + double par = uv[i].Coord( iPar ); + double shift = ShapeAnalysis::AdjustByPeriod( par, max, myHelper->GetPeriod( iPar )); + uv[i].SetCoord( iPar, par + shift ); + } + } + double sign1 = getArea( uv[0], uv[1], uv[2] ); + double sign2 = getArea( uv[0], uv[2], uv[3] ); + if ( sign1 * sign2 < 0 ) + { + sign2 = getArea( uv[1], uv[2], uv[3] ); + sign1 = getArea( uv[1], uv[3], uv[0] ); + if ( sign1 * sign2 < 0 ) + continue; // this should not happen + } + isBad = ( sign1 * okSign < 0 ); + } + if ( isBad ) badFaces.push_back ( f ); }