X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=784853234b4157551317d0220806b50bffa94f06;hp=b2e4a5816de1b1c5b4f8c09f657f2bdbbd8a232c;hb=1071ba77da7f945f7793c8c69a5887d82deb0342;hpb=dcc7add229567325a0ded0ffbf47c13d55b02f1f diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index b2e4a5816..784853234 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2014 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 @@ -6,7 +6,7 @@ // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either -// version 2.1 of the License. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -478,7 +478,10 @@ bool SMESH_MesherHelper::toCheckPosOnShape(int shapeID ) const void SMESH_MesherHelper::setPosOnShapeValidity(int shapeID, bool ok ) const { - ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok)); + std::map< int,bool >::iterator sh_ok = + ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok)).first; + if ( !ok ) + sh_ok->second = ok; } //======================================================================= @@ -661,9 +664,10 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, const bool force, double distXYZ[4]) const { - int shapeID = n->getshapeId(); + int shapeID = n->getshapeId(); bool infinit = ( Precision::IsInfinite( uv.X() ) || Precision::IsInfinite( uv.Y() )); - if ( force || toCheckPosOnShape( shapeID ) || infinit ) + bool zero = ( uv.X() == 0. && uv.Y() == 0. ); + if ( force || toCheckPosOnShape( shapeID ) || infinit || zero ) { // check that uv is correct TopLoc_Location loc; @@ -898,7 +902,8 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, { int shapeID = n->getshapeId(); bool infinit = Precision::IsInfinite( u ); - if ( force || toCheckPosOnShape( shapeID ) || infinit ) + bool zero = ( u == 0. ); + if ( force || toCheckPosOnShape( shapeID ) || infinit || zero ) { TopLoc_Location loc; double f,l; Handle(Geom_Curve) curve = BRep_Tool::Curve( E,loc,f,l ); @@ -1195,10 +1200,11 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1, } else // ( force3d || F.IsNull() ) { - P = ( SMESH_TNodeXYZ( n1 ) + - SMESH_TNodeXYZ( n2 ) + - SMESH_TNodeXYZ( n3 ) + - SMESH_TNodeXYZ( n4 ) ) / 4; + P = calcTFI (0.5, 0.5, + SMESH_TNodeXYZ(n1), SMESH_TNodeXYZ(n2), + SMESH_TNodeXYZ(n3), SMESH_TNodeXYZ(n4), + SMESH_TNodeXYZ(n12), SMESH_TNodeXYZ(n23), + SMESH_TNodeXYZ(n34), SMESH_TNodeXYZ(n41)); centralNode = meshDS->AddNode( P.X(), P.Y(), P.Z() ); if ( !F.IsNull() ) // force3d @@ -1423,8 +1429,16 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, return getMediumNodeOnComposedWire(n1,n2,force3d); } E = TopoDS::Edge(meshDS->IndexToShape( edgeID = pos.first )); - u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]); - u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]); + try { + u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]); + u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]); + } + catch ( Standard_Failure& f ) + { + // issue 22502 / a node is on VERTEX not belonging to E + // issue 22568 / both nodes are on non-connected VERTEXes + return getMediumNodeOnComposedWire(n1,n2,force3d); + } } if ( !force3d & uvOK[0] && uvOK[1] ) @@ -1522,67 +1536,108 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_ const SMDS_MeshNode* n2, bool force3d) { - gp_Pnt middle = 0.5 * XYZ(n1) + 0.5 * XYZ(n2); + SMESH_TNodeXYZ p1( n1 ), p2( n2 ); + gp_Pnt middle = 0.5 * p1 + 0.5 * p2; SMDS_MeshNode* n12 = AddNode( middle.X(), middle.Y(), middle.Z() ); // To find position on edge and 3D position for n12, // project to 2 edges and select projection most close to - double u = 0, distMiddleProj = Precision::Infinite(), distXYZ[4]; - int iOkEdge = 0; - TopoDS_Edge edges[2]; + TopoDS_Edge bestEdge; + double u = 0, distMiddleProj = Precision::Infinite(), distXYZ[4], f,l; + + // get shapes under the nodes + TopoDS_Shape shape[2]; + int nbShapes = 0; for ( int is2nd = 0; is2nd < 2; ++is2nd ) { - // get an edge const SMDS_MeshNode* n = is2nd ? n2 : n1; - TopoDS_Shape shape = GetSubShapeByNode( n, GetMeshDS() ); - if ( shape.IsNull() || shape.ShapeType() != TopAbs_EDGE ) - continue; + TopoDS_Shape S = GetSubShapeByNode( n, GetMeshDS() ); + if ( !S.IsNull() ) + shape[ nbShapes++ ] = S; + } + // get EDGEs + vector< TopoDS_Shape > edges; + for ( int iS = 0; iS < nbShapes; ++iS ) + { + switch ( shape[iS].ShapeType() ) { + case TopAbs_EDGE: + { + edges.push_back( shape[iS] ); + break; + } + case TopAbs_VERTEX: + { + TopoDS_Shape edge; + if ( nbShapes == 2 && iS==0 && shape[1-iS].ShapeType() == TopAbs_VERTEX ) + edge = GetCommonAncestor( shape[iS], shape[1-iS], *myMesh, TopAbs_EDGE ); - // project to get U of projection and distance from middle to projection - TopoDS_Edge edge = edges[ is2nd ] = TopoDS::Edge( shape ); - double node2MiddleDist = middle.Distance( XYZ(n) ); - double foundU = GetNodeU( edge, n ); - CheckNodeU( edge, n12, foundU, 2*BRep_Tool::Tolerance(edge), /*force=*/true, distXYZ ); - if ( distXYZ[0] < node2MiddleDist ) + if ( edge.IsNull() ) + { + PShapeIteratorPtr eIt = GetAncestors( shape[iS], *myMesh, TopAbs_EDGE ); + while( const TopoDS_Shape* e = eIt->next() ) + edges.push_back( *e ); + } + break; + } + case TopAbs_FACE: { - distMiddleProj = distXYZ[0]; - u = foundU; - iOkEdge = is2nd; + if ( nbShapes == 1 || shape[1-iS].ShapeType() < TopAbs_EDGE ) + for ( TopExp_Explorer e( shape[iS], TopAbs_EDGE ); e.More(); e.Next() ) + edges.push_back( e.Current() ); + break; + } + default: + continue; } } - if ( Precision::IsInfinite( distMiddleProj )) + // project to get U of projection and distance from middle to projection + for ( size_t iE = 0; iE < edges.size(); ++iE ) { - // both projections failed; set n12 on the edge of n1 with U of a common vertex - TopoDS_Vertex vCommon; - if ( TopExp::CommonVertex( edges[0], edges[1], vCommon )) - u = BRep_Tool::Parameter( vCommon, edges[0] ); - else + const TopoDS_Edge& edge = TopoDS::Edge( edges[ iE ]); + distXYZ[0] = distMiddleProj; + double testU = 0; + CheckNodeU( edge, n12, testU, 2 * BRep_Tool::Tolerance(edge), /*force=*/true, distXYZ ); + if ( distXYZ[0] < distMiddleProj ) { - double f,l, u0 = GetNodeU( edges[0], n1 ); - BRep_Tool::Range( edges[0],f,l ); - u = ( fabs(u0-f) < fabs(u0-l) ) ? f : l; + distMiddleProj = distXYZ[0]; + u = testU; + bestEdge = edge; } - iOkEdge = 0; - distMiddleProj = 0; - } - - // move n12 to position of a successfull projection - double tol = BRep_Tool::Tolerance(edges[ iOkEdge ]); - if ( !force3d && distMiddleProj > 2*tol ) - { - TopLoc_Location loc; double f,l; - Handle(Geom_Curve) curve = BRep_Tool::Curve( edges[iOkEdge],loc,f,l ); - gp_Pnt p = curve->Value( u ); - GetMeshDS()->MoveNode( n12, p.X(), p.Y(), p.Z() ); } - - //if ( mySetElemOnShape ) node is not elem! + // { + // // both projections failed; set n12 on the edge of n1 with U of a common vertex + // TopoDS_Vertex vCommon; + // if ( TopExp::CommonVertex( edges[0], edges[1], vCommon )) + // u = BRep_Tool::Parameter( vCommon, edges[0] ); + // else + // { + // double f,l, u0 = GetNodeU( edges[0], n1 ); + // BRep_Tool::Range( edges[0],f,l ); + // u = ( fabs(u0-f) < fabs(u0-l) ) ? f : l; + // } + // iOkEdge = 0; + // distMiddleProj = 0; + // } + + if ( !bestEdge.IsNull() ) { - int edgeID = GetMeshDS()->ShapeToIndex( edges[iOkEdge] ); - if ( edgeID != n12->getshapeId() ) - GetMeshDS()->UnSetNodeOnShape( n12 ); - GetMeshDS()->SetNodeOnEdge(n12, edgeID, u); + // move n12 to position of a successfull projection + //double tol = BRep_Tool::Tolerance(edges[ iOkEdge ]); + if ( !force3d /*&& distMiddleProj > 2*tol*/ ) + { + TopLoc_Location loc; + Handle(Geom_Curve) curve = BRep_Tool::Curve( bestEdge,loc,f,l ); + gp_Pnt p = curve->Value( u ).Transformed( loc ); + GetMeshDS()->MoveNode( n12, p.X(), p.Y(), p.Z() ); + } + //if ( mySetElemOnShape ) node is not elem! + { + int edgeID = GetMeshDS()->ShapeToIndex( bestEdge ); + if ( edgeID != n12->getshapeId() ) + GetMeshDS()->UnSetNodeOnShape( n12 ); + GetMeshDS()->SetNodeOnEdge(n12, edgeID, u); + } } myTLinkNodeMap.insert( make_pair( SMESH_TLink(n1,n2), n12 )); @@ -2920,7 +2975,7 @@ namespace { // Structures used by FixQuadraticElements() //======================================================================= #define __DMP__(txt) \ - //cout << txt + // cout << txt #define MSG(txt) __DMP__(txt<GetSubShape()); -// const SMDS_MeshNode* inFaceNode = uvHelper->GetNodeUVneedInFaceNode() ? GetNodeInFace() : 0; -// gp_XY uv1 = uvHelper->GetNodeUV( face, _sides[i]->node1(), inFaceNode ); -// gp_XY uv2 = uvHelper->GetNodeUV( face, _sides[i]->node2(), inFaceNode ); -// norm.SetCoord( uv1.Y() - uv2.Y(), uv2.X() - uv1.X(), 0 ); - -// const QLink* otherLink = _sides[(i + 1) % _sides.size()]; -// const SMDS_MeshNode* otherNode = -// otherLink->node1() == _sides[i]->node1() ? otherLink->node2() : otherLink->node1(); -// gp_XY pIn = uvHelper->GetNodeUV( face, otherNode, inFaceNode ); -// vecOut.SetCoord( uv1.X() - pIn.X(), uv1.Y() - pIn.Y(), 0 ); -// } -// else { - norm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2())); - gp_XYZ pIn = ( XYZ( _sides[0]->node1() ) + - XYZ( _sides[0]->node2() ) + - XYZ( _sides[1]->node1() )) / 3.; - vecOut.SetXYZ( _sides[i]->MiddlePnt() - pIn ); - //} + gp_Vec norm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2())); + gp_XYZ pIn = ( _sides[ (i+1)%3 ]->MiddlePnt() + + _sides[ (i+2)%3 ]->MiddlePnt() ) / 2.; + gp_Vec vecOut = ( _sides[i]->MiddlePnt() - pIn ); + if ( norm * vecOut < 0 ) norm.Reverse(); double mag2 = norm.SquareMagnitude(); @@ -3402,10 +3441,26 @@ namespace { // Structures used by FixQuadraticElements() int iL1 = (iL + 1) % 3, iL2 = (iL + 2) % 3; // indices of the two other links of triangle TLinkInSet link1 = theLinks.find( _sides[iL1] ); TLinkInSet link2 = theLinks.find( _sides[iL2] ); - if ( link1 == theLinks.end() || link2 == theLinks.end() ) - return thePrevLen; - const QFace* f1 = link1->NextFace( this ); // adjacent faces - const QFace* f2 = link2->NextFace( this ); + + const QFace *f1 = 0, *f2 = 0; // adjacent faces + bool isBndLink1 = true, isBndLink2 = true; + if ( link1 != theLinks.end() && link2 != theLinks.end() ) + { + f1 = link1->NextFace( this ); + f2 = link2->NextFace( this ); + + isBndLink1 = ( theLink->MediumPos() > (*link1)->MediumPos() ); + isBndLink2 = ( theLink->MediumPos() > (*link2)->MediumPos() ); + if ( theStep == theFirstStep ) // (issue 22541) quad-dominant mesh + { + if ( !isBndLink1 && !f1 ) + f1 = (*link1)->GetContinuesFace( this ); // get a quadrangle face + if ( !isBndLink2 && !f2 ) + f2 = (*link2)->GetContinuesFace( this ); + } + } + else if ( _sides.size() < 4 ) + return thePrevLen; // propagate to adjacent faces till limit step or boundary double len1 = thePrevLen + (theLink->MiddlePnt() - _sides[iL1]->MiddlePnt()).Modulus(); @@ -3414,7 +3469,7 @@ namespace { // Structures used by FixQuadraticElements() gp_Vec linkDir2(0,0,0); try { OCC_CATCH_SIGNALS; - if ( f1 && theLink->MediumPos() <= (*link1)->MediumPos() ) + if ( f1 && !isBndLink1 ) len1 = f1->MoveByBoundary ( *link1, theRefVec, theLinks, theFaceHelper, len1, theStep-1, &linkDir1, theSign); else @@ -3425,7 +3480,7 @@ namespace { // Structures used by FixQuadraticElements() } try { OCC_CATCH_SIGNALS; - if ( f2 && theLink->MediumPos() <= (*link2)->MediumPos() ) + if ( f2 && !isBndLink2 ) len2 = f2->MoveByBoundary ( *link2, theRefVec, theLinks, theFaceHelper, len2, theStep-1, &linkDir2, theSign); else @@ -3448,7 +3503,7 @@ namespace { // Structures used by FixQuadraticElements() MSG(string(theStep,'.')<<" Move "<< theLink->_mediumNode->GetID()<< " by " << refProj * ( 1 - r ) << " following " << - (choose1 ? *link1->_qlink : *link2->_qlink)); + (choose1 ? *link1->_qlink : *link2->_qlink)); // warning: link1 can be invalid if ( theLinkNorm ) *theLinkNorm = linkNorm; } @@ -3538,7 +3593,7 @@ namespace { // Structures used by FixQuadraticElements() } else if ( _faces.size() > 1 ) // not found, set NULL by the first face { - _faces.insert( ++_faces.begin(), 0 ); + _faces.insert( ++_faces.begin(), (QFace*) 0 ); } } //================================================================================