X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=b0fd1351a89e21a14adeadee55e6f3d1ff56ee4f;hp=13b06a3965960dc8f883e8f0f3f48ab8f360bca5;hb=d4a710ce52f6e76786a7b3845e2f7975dc9a00b1;hpb=9357f5c87098aff2b95b754d69f66c76d2df9c24 diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index 13b06a396..b0fd1351a 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -1,23 +1,23 @@ -// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2012 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 +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS // -// 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. +// 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. // -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. // -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // -// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // // File: SMESH_MesherHelper.cxx @@ -30,6 +30,7 @@ #include "SMDS_EdgePosition.hxx" #include "SMDS_VolumeTool.hxx" #include "SMESH_subMesh.hxx" +#include "SMESH_ProxyMesh.hxx" #include #include @@ -39,6 +40,7 @@ #include #include #include +#include #include #include #include @@ -89,9 +91,16 @@ SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh) SMESH_MesherHelper::~SMESH_MesherHelper() { - TID2Projector::iterator i_proj = myFace2Projector.begin(); - for ( ; i_proj != myFace2Projector.end(); ++i_proj ) - delete i_proj->second; + { + TID2ProjectorOnSurf::iterator i_proj = myFace2Projector.begin(); + for ( ; i_proj != myFace2Projector.end(); ++i_proj ) + delete i_proj->second; + } + { + TID2ProjectorOnCurve::iterator i_proj = myEdge2Projector.begin(); + for ( ; i_proj != myEdge2Projector.end(); ++i_proj ) + delete i_proj->second; + } } //======================================================================= @@ -105,7 +114,7 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) { SMESHDS_Mesh* meshDS = GetMeshDS(); // we can create quadratic elements only if all elements - // created on subshapes of given shape are quadratic + // created on sub-shapes of given shape are quadratic // also we have to fill myTLinkNodeMap myCreateQuadratic = true; mySeamShapeIds.clear(); @@ -115,35 +124,50 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) int nbOldLinks = myTLinkNodeMap.size(); - TopExp_Explorer exp( aSh, subType ); - for (; exp.More() && myCreateQuadratic; exp.Next()) { - if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) { - if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) { - while(it->more()) { - const SMDS_MeshElement* e = it->next(); - if ( e->GetType() != elemType || !e->IsQuadratic() ) { - myCreateQuadratic = false; - break; - } - else { - // fill TLinkNodeMap - switch ( e->NbNodes() ) { - case 3: - AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break; - case 6: - AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3)); - AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4)); - AddTLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break; - case 8: - AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4)); - AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5)); - AddTLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6)); - AddTLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7)); - break; - default: + if ( !myMesh->HasShapeToMesh() ) + { + if (( myCreateQuadratic = myMesh->NbFaces( ORDER_QUADRATIC ))) + { + SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(); + while ( fIt->more() ) + AddTLinks( static_cast< const SMDS_MeshFace* >( fIt->next() )); + } + } + else + { + TopExp_Explorer exp( aSh, subType ); + TopTools_MapOfShape checkedSubShapes; + for (; exp.More() && myCreateQuadratic; exp.Next()) { + if ( !checkedSubShapes.Add( exp.Current() )) + continue; // needed if aSh is compound of solids + if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) { + if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) { + while(it->more()) { + const SMDS_MeshElement* e = it->next(); + if ( e->GetType() != elemType || !e->IsQuadratic() ) { myCreateQuadratic = false; break; } + else { + // fill TLinkNodeMap + switch ( e->NbNodes() ) { + case 3: + AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break; + case 6: + AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3)); + AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4)); + AddTLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break; + case 8: + AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4)); + AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5)); + AddTLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6)); + AddTLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7)); + break; + default: + myCreateQuadratic = false; + break; + } + } } } } @@ -170,7 +194,7 @@ void SMESH_MesherHelper::SetSubShape(const int aShID) { if ( aShID == myShapeID ) return; - if ( aShID > 1 ) + if ( aShID > 0 ) SetSubShape( GetMeshDS()->IndexToShape( aShID )); else SetSubShape( TopoDS_Shape() ); @@ -202,9 +226,16 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() ) { const TopoDS_Face& face = TopoDS::Face( eF.Current() ); - BRepAdaptor_Surface surface( face ); - if ( surface.IsUPeriodic() || surface.IsVPeriodic() ) + TopLoc_Location loc; + Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc ); + + if ( surface->IsUPeriodic() || surface->IsVPeriodic() || + surface->IsUClosed() || surface->IsVClosed() ) { + //while ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface ))) + //surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface(); + GeomAdaptor_Surface surf( surface ); + for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next()) { // look for a seam edge @@ -216,13 +247,13 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) )) { myParIndex |= U_periodic; - myPar1[0] = surface.FirstUParameter(); - myPar2[0] = surface.LastUParameter(); + myPar1[0] = surf.FirstUParameter(); + myPar2[0] = surf.LastUParameter(); } else { myParIndex |= V_periodic; - myPar1[1] = surface.FirstVParameter(); - myPar2[1] = surface.LastVParameter(); + myPar1[1] = surf.FirstVParameter(); + myPar2[1] = surf.LastVParameter(); } // store seam shape indices, negative if shape encounters twice int edgeID = meshDS->ShapeToIndex( edge ); @@ -248,7 +279,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) //function : GetNodeUVneedInFaceNode //purpose : Check if inFaceNode argument is necessary for call GetNodeUV(F,..) // Return true if the face is periodic. -// If F is Null, answer about subshape set through IsQuadraticSubMesh() or +// If F is Null, answer about sub-shape set through IsQuadraticSubMesh() or // * SetSubShape() //======================================================================= @@ -284,9 +315,9 @@ bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode* node, //======================================================================= TopoDS_Shape SMESH_MesherHelper::GetSubShapeByNode(const SMDS_MeshNode* node, - SMESHDS_Mesh* meshDS) + const SMESHDS_Mesh* meshDS) { - int shapeID = node->GetPosition()->GetShapeId(); + int shapeID = node->getshapeId(); if ( 0 < shapeID && shapeID <= meshDS->MaxShapeIndex() ) return meshDS->IndexToShape( shapeID ); else @@ -308,6 +339,100 @@ void SMESH_MesherHelper::AddTLinkNode(const SMDS_MeshNode* n1, myTLinkNodeMap.insert( make_pair(link,n12)); } +//================================================================================ +/*! + * \brief Add quadratic links of edge to own data structure + */ +//================================================================================ + +void SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge) +{ + if ( edge->IsQuadratic() ) + AddTLinkNode(edge->GetNode(0), edge->GetNode(1), edge->GetNode(2)); +} + +//================================================================================ +/*! + * \brief Add quadratic links of face to own data structure + */ +//================================================================================ + +void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f) +{ + if ( !f->IsPoly() ) + switch ( f->NbNodes() ) { + case 6: + AddTLinkNode(f->GetNode(0),f->GetNode(1),f->GetNode(3)); + AddTLinkNode(f->GetNode(1),f->GetNode(2),f->GetNode(4)); + AddTLinkNode(f->GetNode(2),f->GetNode(0),f->GetNode(5)); break; + case 8: + AddTLinkNode(f->GetNode(0),f->GetNode(1),f->GetNode(4)); + AddTLinkNode(f->GetNode(1),f->GetNode(2),f->GetNode(5)); + AddTLinkNode(f->GetNode(2),f->GetNode(3),f->GetNode(6)); + AddTLinkNode(f->GetNode(3),f->GetNode(0),f->GetNode(7)); + default:; + } +} + +//================================================================================ +/*! + * \brief Add quadratic links of volume to own data structure + */ +//================================================================================ + +void SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume) +{ + if ( volume->IsQuadratic() ) + { + SMDS_VolumeTool vTool( volume ); + const SMDS_MeshNode** nodes = vTool.GetNodes(); + set addedLinks; + for ( int iF = 1; iF < vTool.NbFaces(); ++iF ) + { + const int nbN = vTool.NbFaceNodes( iF ); + const int* iNodes = vTool.GetFaceNodesIndices( iF ); + for ( int i = 0; i < nbN; ) + { + int iN1 = iNodes[i++]; + int iN12 = iNodes[i++]; + int iN2 = iNodes[i++]; + if ( iN1 > iN2 ) std::swap( iN1, iN2 ); + int linkID = iN1 * vTool.NbNodes() + iN2; + pair< set::iterator, bool > it_isNew = addedLinks.insert( linkID ); + if ( it_isNew.second ) + AddTLinkNode( nodes[iN1], nodes[iN2], nodes[iN12] ); + else + addedLinks.erase( it_isNew.first ); // each link encounters only twice + } + } + } +} + +//================================================================================ +/*! + * \brief Return true if position of nodes on the shape hasn't yet been checked or + * the positions proved to be invalid + */ +//================================================================================ + +bool SMESH_MesherHelper::toCheckPosOnShape(int shapeID ) const +{ + map< int,bool >::const_iterator id_ok = myNodePosShapesValidity.find( shapeID ); + return ( id_ok == myNodePosShapesValidity.end() || !id_ok->second ); +} + +//================================================================================ +/*! + * \brief Set validity of positions of nodes on the shape. + * Once set, validity is not changed + */ +//================================================================================ + +void SMESH_MesherHelper::setPosOnShapeValidity(int shapeID, bool ok ) const +{ + ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok)); +} + //======================================================================= //function : GetUVOnSeam //purpose : Select UV on either of 2 pcurves of a seam edge, closest to the given UV @@ -353,10 +478,10 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, { // node has position on face const SMDS_FacePosition* fpos = - static_cast(n->GetPosition().get()); + static_cast(n->GetPosition()); uv.SetCoord(fpos->GetUParameter(),fpos->GetVParameter()); if ( check ) - uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), BRep_Tool::Tolerance( F )); + uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F )); } else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE) { @@ -364,8 +489,8 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, // corresponding edge from face, get pcurve for this // edge and retrieve value from this pcurve const SMDS_EdgePosition* epos = - static_cast(n->GetPosition().get()); - int edgeID = Pos->GetShapeId(); + static_cast(n->GetPosition()); + int edgeID = n->getshapeId(); TopoDS_Edge E = TopoDS::Edge(GetMeshDS()->IndexToShape(edgeID)); double f, l, u = epos->GetUParameter(); Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l); @@ -373,14 +498,14 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, if ( validU ) uv = C2d->Value( u ); else - uv.SetCoord(0.,0.); + uv.SetCoord( Precision::Infinite(),0.); if ( check || !validU ) - uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), BRep_Tool::Tolerance( E ),/*force=*/ !validU ); + uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ),/*force=*/ !validU ); // for a node on a seam edge select one of UVs on 2 pcurves if ( n2 && IsSeamShape( edgeID ) ) { - uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 )); + uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0, check )); } else { // adjust uv to period @@ -400,7 +525,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, } else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX) { - if ( int vertexID = n->GetPosition()->GetShapeId() ) { + if ( int vertexID = n->getshapeId() ) { const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID)); try { uv = BRep_Tool::Parameters( V, F ); @@ -452,6 +577,10 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 )); } } + else + { + uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F )); + } if ( check ) *check = uvOK; @@ -468,19 +597,28 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, const SMDS_MeshNode* n, gp_XY& uv, const double tol, - const bool force) const + const bool force, + double distXYZ[4]) const { - if ( force || !myOkNodePosShapes.count( n->GetPosition()->GetShapeId() )) + int shapeID = n->getshapeId(); + bool infinit = ( Precision::IsInfinite( uv.X() ) || Precision::IsInfinite( uv.Y() )); + if ( force || toCheckPosOnShape( shapeID ) || infinit ) { // check that uv is correct TopLoc_Location loc; Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc ); - gp_Pnt nodePnt = XYZ( n ); + gp_Pnt nodePnt = XYZ( n ), surfPnt(0,0,0); + double dist = 0; if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() ); - if ( Precision::IsInfinite( uv.X() ) || - Precision::IsInfinite( uv.Y() ) || - nodePnt.Distance( surface->Value( uv.X(), uv.Y() )) > tol ) + if ( infinit || + (dist = nodePnt.Distance( surfPnt = surface->Value( uv.X(), uv.Y() ))) > tol ) { + setPosOnShapeValidity( shapeID, false ); + if ( !infinit && distXYZ ) { + surfPnt.Transform( loc ); + distXYZ[0] = dist; + distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z(); + } // uv incorrect, project the node to surface GeomAPI_ProjectPointOnSurf& projector = GetProjector( F, loc, tol ); projector.Perform( nodePnt ); @@ -491,16 +629,27 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, } Quantity_Parameter U,V; projector.LowerDistanceParameters(U,V); - if ( nodePnt.Distance( surface->Value( U, V )) > tol ) + uv.SetCoord( U,V ); + surfPnt = surface->Value( U, V ); + dist = nodePnt.Distance( surfPnt ); + if ( distXYZ ) { + surfPnt.Transform( loc ); + distXYZ[0] = dist; + distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z(); + } + if ( dist > tol ) { MESSAGE( "SMESH_MesherHelper::CheckNodeUV(), invalid projection" ); return false; } - uv.SetCoord( U,V ); + // store the fixed UV on the face + if ( myShape.IsSame(F) && shapeID == myShapeID ) + const_cast(n)->SetPosition + ( SMDS_PositionPtr( new SMDS_FacePosition( U, V ))); } else if ( uv.Modulus() > numeric_limits::min() ) { - ((SMESH_MesherHelper*) this)->myOkNodePosShapes.insert( n->GetPosition()->GetShapeId() ); + setPosOnShapeValidity( shapeID, true ); } } return true; @@ -517,8 +666,8 @@ GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face& { Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc ); int faceID = GetMeshDS()->ShapeToIndex( F ); - TID2Projector& i2proj = const_cast< TID2Projector&>( myFace2Projector ); - TID2Projector::iterator i_proj = i2proj.find( faceID ); + TID2ProjectorOnSurf& i2proj = const_cast< TID2ProjectorOnSurf&>( myFace2Projector ); + TID2ProjectorOnSurf::iterator i_proj = i2proj.find( faceID ); if ( i_proj == i2proj.end() ) { if ( tol == 0 ) tol = BRep_Tool::Tolerance( F ); @@ -587,7 +736,14 @@ gp_XY SMESH_MesherHelper::GetMiddleUV(const Handle(Geom_Surface)& surface, const gp_XY& p1, const gp_XY& p2) { - return applyIn2D( surface, p1, p2, & AverageUV ); + // NOTE: + // the proper place of getting basic surface seems to be in applyIn2D() + // but we put it here to decrease a risk of regressions just before releasing a version + Handle(Geom_Surface) surf = surface; + while ( !surf.IsNull() && surf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface ))) + surf = Handle(Geom_RectangularTrimmedSurface)::DownCast( surf )->BasisSurface(); + + return applyIn2D( surf, p1, p2, & AverageUV ); } //======================================================================= @@ -597,23 +753,43 @@ gp_XY SMESH_MesherHelper::GetMiddleUV(const Handle(Geom_Surface)& surface, double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge& E, const SMDS_MeshNode* n, + const SMDS_MeshNode* inEdgeNode, bool* check) { double param = 0; - const SMDS_PositionPtr Pos = n->GetPosition(); - if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE) { - const SMDS_EdgePosition* epos = - static_cast(n->GetPosition().get()); + const SMDS_PositionPtr pos = n->GetPosition(); + if ( pos->GetTypeOfPosition()==SMDS_TOP_EDGE ) + { + const SMDS_EdgePosition* epos = static_cast( pos ); param = epos->GetUParameter(); } - else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX) { - SMESHDS_Mesh * meshDS = GetMeshDS(); - int vertexID = n->GetPosition()->GetShapeId(); - const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID)); - param = BRep_Tool::Parameter( V, E ); + else if( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX ) + { + if ( inEdgeNode && TopExp::FirstVertex( E ).IsSame( TopExp::LastVertex( E ))) // issue 0020128 + { + Standard_Real f,l; + BRep_Tool::Range( E, f,l ); + double uInEdge = GetNodeU( E, inEdgeNode ); + param = ( fabs( uInEdge - f ) < fabs( l - uInEdge )) ? f : l; + } + else + { + SMESHDS_Mesh * meshDS = GetMeshDS(); + int vertexID = n->getshapeId(); + const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID)); + param = BRep_Tool::Parameter( V, E ); + } } if ( check ) - *check = CheckNodeU( E, n, param, BRep_Tool::Tolerance( E )); + { + double tol = BRep_Tool::Tolerance( E ); + double f,l; BRep_Tool::Range( E, f,l ); + bool force = ( param < f-tol || param > l+tol ); + if ( !force && pos->GetTypeOfPosition()==SMDS_TOP_EDGE ) + force = ( GetMeshDS()->ShapeToIndex( E ) != n->getshapeId() ); + + *check = CheckNodeU( E, n, param, 2*tol, force ); + } return param; } @@ -627,11 +803,12 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, const SMDS_MeshNode* n, double& u, const double tol, - const bool force) const + const bool force, + double distXYZ[4]) const { - if ( force || !myOkNodePosShapes.count( n->GetPosition()->GetShapeId() )) + int shapeID = n->getshapeId(); + if ( force || toCheckPosOnShape( shapeID )) { - // check that u is correct TopLoc_Location loc; double f,l; Handle(Geom_Curve) curve = BRep_Tool::Curve( E,loc,f,l ); if ( curve.IsNull() ) // degenerated edge @@ -644,34 +821,157 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, } else { - gp_Pnt nodePnt = SMESH_MeshEditor::TNodeXYZ( n ); + gp_Pnt nodePnt = SMESH_TNodeXYZ( n ); if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() ); - if ( nodePnt.Distance( curve->Value( u )) > tol ) + gp_Pnt curvPnt = curve->Value( u ); + double dist = nodePnt.Distance( curvPnt ); + if ( distXYZ ) { + curvPnt.Transform( loc ); + distXYZ[0] = dist; + distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z(); + } + if ( dist > tol ) { + setPosOnShapeValidity( shapeID, false ); // u incorrect, project the node to the curve - GeomAPI_ProjectPointOnCurve projector( nodePnt, curve, f, l ); - if ( projector.NbPoints() < 1 ) + int edgeID = GetMeshDS()->ShapeToIndex( E ); + TID2ProjectorOnCurve& i2proj = const_cast< TID2ProjectorOnCurve&>( myEdge2Projector ); + TID2ProjectorOnCurve::iterator i_proj = + i2proj.insert( make_pair( edgeID, (GeomAPI_ProjectPointOnCurve*) 0 )).first; + if ( !i_proj->second ) + { + i_proj->second = new GeomAPI_ProjectPointOnCurve(); + i_proj->second->Init( curve, f, l ); + } + GeomAPI_ProjectPointOnCurve* projector = i_proj->second; + projector->Perform( nodePnt ); + if ( projector->NbPoints() < 1 ) { MESSAGE( "SMESH_MesherHelper::CheckNodeU() failed to project" ); return false; } - Quantity_Parameter U = projector.LowerDistanceParameter(); - if ( nodePnt.Distance( curve->Value( U )) > tol ) + Quantity_Parameter U = projector->LowerDistanceParameter(); + u = double( U ); + curvPnt = curve->Value( u ); + dist = nodePnt.Distance( curvPnt ); + if ( distXYZ ) { + curvPnt.Transform( loc ); + distXYZ[0] = dist; + distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z(); + } + if ( dist > tol ) { MESSAGE( "SMESH_MesherHelper::CheckNodeU(), invalid projection" ); + MESSAGE("distance " << dist << " " << tol ); return false; } - u = double( U ); + // store the fixed U on the edge + if ( myShape.IsSame(E) && shapeID == myShapeID ) + const_cast(n)->SetPosition + ( SMDS_PositionPtr( new SMDS_EdgePosition( U ))); } else if ( fabs( u ) > numeric_limits::min() ) { - ((SMESH_MesherHelper*) this)->myOkNodePosShapes.insert( n->GetPosition()->GetShapeId() ); + setPosOnShapeValidity( shapeID, true ); + } + if (( u < f-tol || u > l+tol ) && force ) + { + // node is on vertex but is set on periodic but trimmed edge (issue 0020890) + try + { + // do not use IsPeriodic() as Geom_TrimmedCurve::IsPeriodic () returns false + double period = curve->Period(); + u = ( u < f ) ? u + period : u - period; + } + catch (Standard_Failure& exc) + { + return false; + } } } } return true; } +//======================================================================= +//function : GetMediumPos +//purpose : Return index and type of the shape (EDGE or FACE only) to +// set a medium node on +//======================================================================= + +std::pair SMESH_MesherHelper::GetMediumPos(const SMDS_MeshNode* n1, + const SMDS_MeshNode* n2) +{ + TopAbs_ShapeEnum shapeType = TopAbs_SHAPE; + int shapeID = -1; + TopoDS_Shape shape; + + if (( myShapeID == n1->getshapeId() || myShapeID == n2->getshapeId() ) && myShapeID > 0 ) + { + shapeType = myShape.ShapeType(); + shapeID = myShapeID; + } + else if ( n1->getshapeId() == n2->getshapeId() ) + { + shapeID = n2->getshapeId(); + shape = GetSubShapeByNode( n1, GetMeshDS() ); + } + else + { + const SMDS_TypeOfPosition Pos1 = n1->GetPosition()->GetTypeOfPosition(); + const SMDS_TypeOfPosition Pos2 = n2->GetPosition()->GetTypeOfPosition(); + + if ( Pos1 == SMDS_TOP_3DSPACE || Pos2 == SMDS_TOP_3DSPACE ) + { + } + else if ( Pos1 == SMDS_TOP_FACE || Pos2 == SMDS_TOP_FACE ) + { + if ( Pos1 != SMDS_TOP_FACE || Pos2 != SMDS_TOP_FACE ) + { + if ( Pos1 != SMDS_TOP_FACE ) std::swap( n1,n2 ); + TopoDS_Shape F = GetSubShapeByNode( n1, GetMeshDS() ); + TopoDS_Shape S = GetSubShapeByNode( n2, GetMeshDS() ); + if ( IsSubShape( S, F )) + { + shapeType = TopAbs_FACE; + shapeID = n1->getshapeId(); + } + } + } + else if ( Pos1 == SMDS_TOP_EDGE && Pos2 == SMDS_TOP_EDGE ) + { + TopoDS_Shape E1 = GetSubShapeByNode( n1, GetMeshDS() ); + TopoDS_Shape E2 = GetSubShapeByNode( n2, GetMeshDS() ); + shape = GetCommonAncestor( E1, E2, *myMesh, TopAbs_FACE ); + } + else if ( Pos1 == SMDS_TOP_VERTEX && Pos2 == SMDS_TOP_VERTEX ) + { + TopoDS_Shape V1 = GetSubShapeByNode( n1, GetMeshDS() ); + TopoDS_Shape V2 = GetSubShapeByNode( n2, GetMeshDS() ); + shape = GetCommonAncestor( V1, V2, *myMesh, TopAbs_EDGE ); + if ( shape.IsNull() ) shape = GetCommonAncestor( V1, V2, *myMesh, TopAbs_FACE ); + } + else // VERTEX and EDGE + { + if ( Pos1 != SMDS_TOP_VERTEX ) std::swap( n1,n2 ); + TopoDS_Shape V = GetSubShapeByNode( n1, GetMeshDS() ); + TopoDS_Shape E = GetSubShapeByNode( n2, GetMeshDS() ); + if ( IsSubShape( V, E )) + shape = E; + else + shape = GetCommonAncestor( V, E, *myMesh, TopAbs_FACE ); + } + } + + if ( !shape.IsNull() ) + { + if ( shapeID < 1 ) + shapeID = GetMeshDS()->ShapeToIndex( shape ); + shapeType = shape.ShapeType(); + } + return make_pair( shapeID, shapeType ); +} + //======================================================================= //function : GetMediumNode //purpose : Return existing or create new medium nodes between given ones @@ -694,68 +994,57 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, SMDS_MeshNode* n12; SMESHDS_Mesh* meshDS = GetMeshDS(); + if ( IsSeamShape( n1->getshapeId() )) + // to get a correct UV of a node on seam, the second node must have checked UV + std::swap( n1, n2 ); + // get type of shape for the new medium node int faceID = -1, edgeID = -1; - const SMDS_PositionPtr Pos1 = n1->GetPosition(); - const SMDS_PositionPtr Pos2 = n2->GetPosition(); - - if( myShape.IsNull() ) - { - if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) { - faceID = Pos1->GetShapeId(); - } - else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) { - faceID = Pos2->GetShapeId(); - } - - if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) { - edgeID = Pos1->GetShapeId(); - } - if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) { - edgeID = Pos2->GetShapeId(); - } - } - // get positions of the given nodes on shapes TopoDS_Edge E; double u [2]; TopoDS_Face F; gp_XY uv[2]; bool uvOK[2] = { false, false }; - TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType(); - if ( faceID>0 || shapeType == TopAbs_FACE) + + pair pos = GetMediumPos( n1, n2 ); + + // get positions of the given nodes on shapes + if ( pos.second == TopAbs_FACE ) { - if( myShape.IsNull() ) - F = TopoDS::Face(meshDS->IndexToShape(faceID)); - else { - F = TopoDS::Face(myShape); - faceID = myShapeID; - } + F = TopoDS::Face(meshDS->IndexToShape( faceID = pos.first )); uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]); uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]); } - else if (edgeID>0 || shapeType == TopAbs_EDGE) + else if ( pos.second == TopAbs_EDGE ) { - if( myShape.IsNull() ) - E = TopoDS::Edge(meshDS->IndexToShape(edgeID)); - else { - E = TopoDS::Edge(myShape); - edgeID = myShapeID; + const SMDS_PositionPtr Pos1 = n1->GetPosition(); + const SMDS_PositionPtr Pos2 = n2->GetPosition(); + if ( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE && + Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE && + n1->getshapeId() != n2->getshapeId() ) + { + // issue 0021006 + return getMediumNodeOnComposedWire(n1,n2,force3d); } - u[0] = GetNodeU(E,n1, force3d ? 0 : &uvOK[0]); - u[1] = GetNodeU(E,n2, force3d ? 0 : &uvOK[1]); + 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]); } - if(!force3d) + + if ( !force3d & uvOK[0] && uvOK[1] ) { // we try to create medium node using UV parameters of // nodes, else - medium between corresponding 3d points if( ! F.IsNull() ) { - if ( uvOK[0] && uvOK[1] ) + //if ( uvOK[0] && uvOK[1] ) { - if ( IsDegenShape( Pos1->GetShapeId() )) + 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( Pos2->GetShapeId() )) + } + 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); @@ -792,31 +1081,106 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, } } } + // 3d variant double x = ( n1->X() + n2->X() )/2.; double y = ( n1->Y() + n2->Y() )/2.; double z = ( n1->Z() + n2->Z() )/2.; n12 = meshDS->AddNode(x,y,z); + if ( !F.IsNull() ) { gp_XY UV = ( uv[0] + uv[1] ) / 2.; - CheckNodeUV( F, n12, UV, BRep_Tool::Tolerance( F ), /*force=*/true); + CheckNodeUV( F, n12, UV, 2*BRep_Tool::Tolerance( F ), /*force=*/true); meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y() ); } else if ( !E.IsNull() ) { double U = ( u[0] + u[1] ) / 2.; - CheckNodeU( E, n12, U, BRep_Tool::Tolerance( E ), /*force=*/true); + CheckNodeU( E, n12, U, 2*BRep_Tool::Tolerance( E ), /*force=*/true); meshDS->SetNodeOnEdge(n12, edgeID, U); } - else if ( myShapeID > 1 ) + else if ( myShapeID > 0 ) { meshDS->SetNodeInVolume(n12, myShapeID); } + myTLinkNodeMap.insert( make_pair( link, n12 )); return n12; } +//================================================================================ +/*! + * \brief Makes a medium node if nodes reside different edges + */ +//================================================================================ + +const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_MeshNode* n1, + const SMDS_MeshNode* n2, + bool force3d) +{ + gp_Pnt middle = 0.5 * XYZ(n1) + 0.5 * XYZ(n2); + 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]; + 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; + + // 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 ) + { + distMiddleProj = distXYZ[0]; + u = foundU; + iOkEdge = is2nd; + } + } + if ( Precision::IsInfinite( distMiddleProj )) + { + // 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; + } + + // 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() ); + } + + GetMeshDS()->SetNodeOnEdge(n12, edges[iOkEdge], u); + + myTLinkNodeMap.insert( make_pair( SMESH_TLink(n1,n2), n12 )); + + return n12; +} + //======================================================================= //function : AddNode //purpose : Creates a node @@ -972,6 +1336,45 @@ SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1, return elem; } +//======================================================================= +//function : AddPolygonalFace +//purpose : Creates polygon, with additional nodes in quadratic mesh +//======================================================================= + +SMDS_MeshFace* SMESH_MesherHelper::AddPolygonalFace (const vector& nodes, + const int id, + const bool force3d) +{ + SMESHDS_Mesh * meshDS = GetMeshDS(); + SMDS_MeshFace* elem = 0; + + if(!myCreateQuadratic) { + if(id) + elem = meshDS->AddPolygonalFaceWithID(nodes, id); + else + elem = meshDS->AddPolygonalFace(nodes); + } + else { + vector newNodes; + for ( int i = 0; i < nodes.size(); ++i ) + { + const SMDS_MeshNode* n1 = nodes[i]; + const SMDS_MeshNode* n2 = nodes[(i+1)%nodes.size()]; + const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d); + newNodes.push_back( n1 ); + newNodes.push_back( n12 ); + } + if(id) + elem = meshDS->AddPolygonalFaceWithID(newNodes, id); + else + elem = meshDS->AddPolygonalFace(newNodes); + } + if ( mySetElemOnShape && myShapeID > 0 ) + meshDS->SetMeshElementOnShape( elem, myShapeID ); + + return elem; +} + //======================================================================= //function : AddVolume //purpose : Creates quadratic or linear prism @@ -1162,6 +1565,62 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, return elem; } +//======================================================================= +//function : AddPolyhedralVolume +//purpose : Creates polyhedron. In quadratic mesh, adds medium nodes +//======================================================================= + +SMDS_MeshVolume* +SMESH_MesherHelper::AddPolyhedralVolume (const std::vector& nodes, + const std::vector& quantities, + const int id, + const bool force3d) +{ + SMESHDS_Mesh * meshDS = GetMeshDS(); + SMDS_MeshVolume* elem = 0; + if(!myCreateQuadratic) + { + if(id) + elem = meshDS->AddPolyhedralVolumeWithID(nodes, quantities, id); + else + elem = meshDS->AddPolyhedralVolume(nodes, quantities); + } + else + { + vector newNodes; + vector newQuantities; + for ( int iFace=0, iN=0; iFace < quantities.size(); ++iFace) + { + int nbNodesInFace = quantities[iFace]; + newQuantities.push_back(0); + for ( int i = 0; i < nbNodesInFace; ++i ) + { + const SMDS_MeshNode* n1 = nodes[ iN + i ]; + newNodes.push_back( n1 ); + newQuantities.back()++; + + const SMDS_MeshNode* n2 = nodes[ iN + ( i+1==nbNodesInFace ? 0 : i+1 )]; +// if ( n1->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE && +// n2->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE ) + { + const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d); + newNodes.push_back( n12 ); + newQuantities.back()++; + } + } + iN += nbNodesInFace; + } + if(id) + elem = meshDS->AddPolyhedralVolumeWithID( newNodes, newQuantities, id ); + else + elem = meshDS->AddPolyhedralVolume( newNodes, newQuantities ); + } + if ( mySetElemOnShape && myShapeID > 0 ) + meshDS->SetMeshElementOnShape( elem, myShapeID ); + + return elem; +} + //======================================================================= //function : LoadNodeColumns //purpose : Load nodes bound to face into a map of node columns @@ -1170,40 +1629,107 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap, const TopoDS_Face& theFace, const TopoDS_Edge& theBaseEdge, - SMESHDS_Mesh* theMesh) + SMESHDS_Mesh* theMesh, + SMESH_ProxyMesh* theProxyMesh) { - SMESHDS_SubMesh* faceSubMesh = theMesh->MeshElements( theFace ); + return LoadNodeColumns(theParam2ColumnMap, + theFace, + std::list(1,theBaseEdge), + theMesh, + theProxyMesh); +} + +//======================================================================= +//function : LoadNodeColumns +//purpose : Load nodes bound to face into a map of node columns +//======================================================================= + +bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap, + const TopoDS_Face& theFace, + const std::list& theBaseSide, + SMESHDS_Mesh* theMesh, + SMESH_ProxyMesh* theProxyMesh) +{ + // get a right submesh of theFace + + const SMESHDS_SubMesh* faceSubMesh = 0; + if ( theProxyMesh ) + { + faceSubMesh = theProxyMesh->GetSubMesh( theFace ); + if ( !faceSubMesh || + faceSubMesh->NbElements() == 0 || + theProxyMesh->IsTemporary( faceSubMesh->GetElements()->next() )) + { + // can use a proxy sub-mesh with not temporary elements only + faceSubMesh = 0; + theProxyMesh = 0; + } + } + if ( !faceSubMesh ) + faceSubMesh = theMesh->MeshElements( theFace ); if ( !faceSubMesh || faceSubMesh->NbElements() == 0 ) return false; - // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them + // get data of edges for normalization of params - map< double, const SMDS_MeshNode*> sortedBaseNodes; - if ( !SMESH_Algo::GetSortedNodesOnEdge( theMesh, theBaseEdge,/*noMedium=*/true, sortedBaseNodes) - || sortedBaseNodes.size() < 2 ) - return false; + vector< double > length; + double fullLen = 0; + list::const_iterator edge; + { + for ( edge = theBaseSide.begin(); edge != theBaseSide.end(); ++edge ) + { + double len = std::max( 1e-10, SMESH_Algo::EdgeLength( *edge )); + fullLen += len; + length.push_back( len ); + } + } - int nbRows = faceSubMesh->NbElements() / ( sortedBaseNodes.size()-1 ) + 1; - map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin(); - double f = u_n->first, range = sortedBaseNodes.rbegin()->first - f; - for ( ; u_n != sortedBaseNodes.end(); u_n++ ) + // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them + edge = theBaseSide.begin(); + for ( int iE = 0; edge != theBaseSide.end(); ++edge, ++iE ) { - double par = ( u_n->first - f ) / range; - vector& nCol = theParam2ColumnMap[ par ]; - nCol.resize( nbRows ); - nCol[0] = u_n->second; + map< double, const SMDS_MeshNode*> sortedBaseNodes; + SMESH_Algo::GetSortedNodesOnEdge( theMesh, *edge,/*noMedium=*/true, sortedBaseNodes); + if ( sortedBaseNodes.empty() ) continue; + + double f, l; + BRep_Tool::Range( *edge, f, l ); + if ( edge->Orientation() == TopAbs_REVERSED ) std::swap( f, l ); + const double coeff = 1. / ( l - f ) * length[iE] / fullLen; + const double prevPar = theParam2ColumnMap.empty() ? 0 : theParam2ColumnMap.rbegin()->first; + map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin(); + for ( ; u_n != sortedBaseNodes.end(); u_n++ ) + { + double par = prevPar + coeff * ( u_n->first - f ); + TParam2ColumnMap::iterator u2nn = + theParam2ColumnMap.insert( theParam2ColumnMap.end(), make_pair( par, TNodeColumn())); + u2nn->second.push_back( u_n->second ); + } + } + TParam2ColumnMap::iterator par_nVec_2, par_nVec_1 = theParam2ColumnMap.begin(); + if ( theProxyMesh ) + { + for ( ; par_nVec_1 != theParam2ColumnMap.end(); ++par_nVec_1 ) + { + const SMDS_MeshNode* & n = par_nVec_1->second[0]; + n = theProxyMesh->GetProxyNode( n ); + } } + int nbRows = 1 + faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 ); + // fill theParam2ColumnMap column by column by passing from nodes on // theBaseEdge up via mesh faces on theFace - TParam2ColumnMap::iterator par_nVec_2 = theParam2ColumnMap.begin(); - TParam2ColumnMap::iterator par_nVec_1 = par_nVec_2++; + par_nVec_2 = theParam2ColumnMap.begin(); + par_nVec_1 = par_nVec_2++; TIDSortedElemSet emptySet, avoidSet; for ( ; par_nVec_2 != theParam2ColumnMap.end(); ++par_nVec_1, ++par_nVec_2 ) { vector& nCol1 = par_nVec_1->second; vector& nCol2 = par_nVec_2->second; + nCol1.resize( nbRows ); + nCol2.resize( nbRows ); int i1, i2, iRow = 0; const SMDS_MeshNode *n1 = nCol1[0], *n2 = nCol2[0]; @@ -1226,10 +1752,11 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap, } avoidSet.insert( face ); } - if ( iRow + 1 < nbRows ) // compact if necessary - nCol1.resize( iRow + 1 ), nCol2.resize( iRow + 1 ); + // set a real height + nCol1.resize( iRow + 1 ); + nCol2.resize( iRow + 1 ); } - return true; + return theParam2ColumnMap.size() > 1 && theParam2ColumnMap.begin()->second.size() > 1; } //======================================================================= @@ -1306,7 +1833,64 @@ bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, SMESH_Mesh* aMes return aMesh->GetMeshDS()->ShapeToIndex( shape ) || // PAL16202 - shape.ShapeType() == TopAbs_COMPOUND && aMesh->GetMeshDS()->IsGroupOfSubShapes( shape ); + (shape.ShapeType() == TopAbs_COMPOUND && aMesh->GetMeshDS()->IsGroupOfSubShapes( shape )); +} + +//================================================================================ +/*! + * \brief Return maximal tolerance of shape + */ +//================================================================================ + +double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape ) +{ + double tol = Precision::Confusion(); + TopExp_Explorer exp; + for ( exp.Init( shape, TopAbs_FACE ); exp.More(); exp.Next() ) + tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Face( exp.Current()))); + for ( exp.Init( shape, TopAbs_EDGE ); exp.More(); exp.Next() ) + tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Edge( exp.Current()))); + for ( exp.Init( shape, TopAbs_VERTEX ); exp.More(); exp.Next() ) + tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Vertex( exp.Current()))); + + return tol; +} + +//================================================================================ +/*! + * \brief Check if the first and last vertices of an edge are the same + * \param anEdge - the edge to check + * \retval bool - true if same + */ +//================================================================================ + +bool SMESH_MesherHelper::IsClosedEdge( const TopoDS_Edge& anEdge ) +{ + if ( anEdge.Orientation() >= TopAbs_INTERNAL ) + return IsClosedEdge( TopoDS::Edge( anEdge.Oriented( TopAbs_FORWARD ))); + return TopExp::FirstVertex( anEdge ).IsSame( TopExp::LastVertex( anEdge )); +} + +//================================================================================ +/*! + * \brief Wrapper over TopExp::FirstVertex() and TopExp::LastVertex() fixing them + * in the case of INTERNAL edge + */ +//================================================================================ + +TopoDS_Vertex SMESH_MesherHelper::IthVertex( const bool is2nd, + TopoDS_Edge anEdge, + const bool CumOri ) +{ + if ( anEdge.Orientation() >= TopAbs_INTERNAL ) + anEdge.Orientation( TopAbs_FORWARD ); + + const TopAbs_Orientation tgtOri = is2nd ? TopAbs_REVERSED : TopAbs_FORWARD; + TopoDS_Iterator vIt( anEdge, CumOri ); + while ( vIt.More() && vIt.Value().Orientation() != tgtOri ) + vIt.Next(); + + return ( vIt.More() ? TopoDS::Vertex(vIt.Value()) : TopoDS_Vertex() ); } //======================================================================= @@ -1354,16 +1938,98 @@ double SMESH_MesherHelper::GetOtherParam(const double param) const return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i]; } +namespace { + + //======================================================================= + /*! + * \brief Iterator on ancestors of the given type + */ + //======================================================================= + + struct TAncestorsIterator : public SMDS_Iterator + { + TopTools_ListIteratorOfListOfShape _ancIter; + TopAbs_ShapeEnum _type; + TopTools_MapOfShape _encountered; + TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type) + : _ancIter( ancestors ), _type( type ) + { + if ( _ancIter.More() ) { + if ( _ancIter.Value().ShapeType() != _type ) next(); + else _encountered.Add( _ancIter.Value() ); + } + } + virtual bool more() + { + return _ancIter.More(); + } + virtual const TopoDS_Shape* next() + { + const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0; + if ( _ancIter.More() ) + for ( _ancIter.Next(); _ancIter.More(); _ancIter.Next()) + if ( _ancIter.Value().ShapeType() == _type && _encountered.Add( _ancIter.Value() )) + break; + return s; + } + }; + +} // namespace + +//======================================================================= +/*! + * \brief Return iterator on ancestors of the given type + */ +//======================================================================= + +PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape, + const SMESH_Mesh& mesh, + TopAbs_ShapeEnum ancestorType) +{ + return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType)); +} + +//======================================================================= +//function : GetCommonAncestor +//purpose : Find a common ancestors of two shapes of the given type +//======================================================================= + +TopoDS_Shape SMESH_MesherHelper::GetCommonAncestor(const TopoDS_Shape& shape1, + const TopoDS_Shape& shape2, + const SMESH_Mesh& mesh, + TopAbs_ShapeEnum ancestorType) +{ + TopoDS_Shape commonAnc; + if ( !shape1.IsNull() && !shape2.IsNull() ) + { + PShapeIteratorPtr ancIt = GetAncestors( shape1, mesh, ancestorType ); + while ( const TopoDS_Shape* anc = ancIt->next() ) + if ( IsSubShape( shape2, *anc )) + { + commonAnc = *anc; + break; + } + } + return commonAnc; +} + +//#include + //======================================================================= namespace { // Structures used by FixQuadraticElements() //======================================================================= #define __DMP__(txt) \ -//cout << txt + //cout << txt #define MSG(txt) __DMP__(txt< < 1/15 * + return middleNodeMove2 < 1/15./15. * linkLen2; + } struct QFace; // --------------------------------------- @@ -1399,17 +2065,19 @@ namespace { // Structures used by FixQuadraticElements() void Move(const gp_Vec& move, bool sum=false) const { _nodeMove += move; _nbMoves += sum ? (_nbMoves==0) : 1; } gp_XYZ Move() const { return _nodeMove.XYZ() / _nbMoves; } - bool IsMoved() const { return (_nbMoves > 0 && !IsStraight()); } - bool IsStraight() const { return _nodeMove.SquareMagnitude() <= straightTol2; } - + bool IsMoved() const { return (_nbMoves > 0 /*&& !IsStraight()*/); } + bool IsStraight() const + { return isStraightLink( (XYZ(node1())-XYZ(node2())).SquareModulus(), + _nodeMove.SquareMagnitude()); + } bool operator<(const QLink& other) const { return (node1()->GetID() == other.node1()->GetID() ? node2()->GetID() < other.node2()->GetID() : node1()->GetID() < other.node1()->GetID()); } - struct PtrComparator { - bool operator() (const QLink* l1, const QLink* l2 ) const { return *l1 < *l2; } - }; +// struct PtrComparator { +// bool operator() (const QLink* l1, const QLink* l2 ) const { return *l1 < *l2; } +// }; }; // --------------------------------------------------------- /*! @@ -1423,7 +2091,7 @@ namespace { // Structures used by FixQuadraticElements() TChainLink(const QLink* qlink=0):_qlink(qlink) { _qfaces[0] = _qfaces[1] = 0; } - void SetFace(const QFace* face) { int iF = _qfaces[0] ? 1 : 0; _qfaces[iF]=face; } + void SetFace(const QFace* face) const { int iF = _qfaces[0] ? 1 : 0; _qfaces[iF]=face; } bool IsBoundary() const { return !_qfaces[1]; } @@ -1443,6 +2111,8 @@ namespace { // Structures used by FixQuadraticElements() const QLink* operator->() const { return _qlink; } gp_Vec Normal() const; + + bool IsStraight() const; }; // -------------------------------------------------------------------- typedef list< TChainLink > TChain; @@ -1456,7 +2126,7 @@ namespace { // Structures used by FixQuadraticElements() /*! * \brief Face shared by two volumes and bound by QLinks */ - struct QFace: public TIDSortedElemSet + struct QFace: public TIDSortedNodeSet { mutable const SMDS_MeshElement* _volumes[2]; mutable vector< const QLink* > _sides; @@ -1495,6 +2165,8 @@ namespace { // Structures used by FixQuadraticElements() bool Contains( const SMDS_MeshNode* node ) const { return count(node); } + bool IsSpoiled(const QLink* bentLink ) const; + TLinkInSet GetBoundaryLink( const TLinkSet& links, const TChainLink& avoidLink, TLinkInSet * notBoundaryLink = 0, @@ -1539,7 +2211,7 @@ namespace { // Structures used by FixQuadraticElements() ostream& operator << (ostream& out, const QFace& f) { out <<"QFace nodes: "/*<< &f << " "*/; - for ( TIDSortedElemSet::const_iterator n = f.begin(); n != f.end(); ++n ) + for ( TIDSortedNodeSet::const_iterator n = f.begin(); n != f.end(); ++n ) out << (*n)->GetID() << " "; out << " \tvolumes: " << (f._volumes[0] ? f._volumes[0]->GetID() : 0) << " " @@ -1582,7 +2254,7 @@ namespace { // Structures used by FixQuadraticElements() } //================================================================================ /*! - * \brief Make up chain of links + * \brief Make up a chain of links * \param iSide - link to add first * \param chain - chain to fill in * \param pos - postion of medium nodes the links should have @@ -1601,22 +2273,34 @@ namespace { // Structures used by FixQuadraticElements() if ( _sides.size() != 4 ) { // triangle - visit all my continous faces MSGBEG( *this ); + TLinkSet links; list< const QFace* > faces( 1, this ); - for (list< const QFace* >::iterator fIt = faces.begin(); fIt != faces.end(); ++fIt ) { - const QFace* face = *fIt; + while ( !faces.empty() ) { + const QFace* face = faces.front(); for ( int i = 0; i < face->_sides.size(); ++i ) { if ( !face->_sideIsAdded[i] && face->_sides[i] ) { face->_sideIsAdded[i] = true; - TChain::iterator chLink = chain.insert( chain.begin(), TChainLink(face->_sides[i])); + // find a face side in the chain + TLinkInSet chLink = links.insert( TChainLink(face->_sides[i])).first; +// TChain::iterator chLink = chain.begin(); +// for ( ; chLink != chain.end(); ++chLink ) +// if ( chLink->_qlink == face->_sides[i] ) +// break; +// if ( chLink == chain.end() ) +// chLink = chain.insert( chain.begin(), TChainLink(face->_sides[i])); + // add a face to a chained link and put a continues face in the queue chLink->SetFace( face ); - if ( face->_sides[i]->MediumPos() >= pos ) + if ( face->_sides[i]->MediumPos() == pos ) if ( const QFace* contFace = face->_sides[i]->GetContinuesFace( face )) - faces.push_back( contFace ); + if ( contFace->_sides.size() == 3 ) + faces.push_back( contFace ); } } + faces.pop_front(); } if ( error < ERR_TRI ) error = ERR_TRI; + chain.insert( chain.end(), links.begin(),links.end() ); return false; } _sideIsAdded[iSide] = true; // not to add this link to chain again @@ -1632,10 +2316,11 @@ namespace { // Structures used by FixQuadraticElements() // propagate from quadrangle to neighbour faces if ( link->MediumPos() >= pos ) { int nbLinkFaces = link->_faces.size(); - if ( nbLinkFaces == 4 || nbLinkFaces < 4 && link->OnBoundary()) { + if ( nbLinkFaces == 4 || (/*nbLinkFaces < 4 && */link->OnBoundary())) { // hexahedral mesh or boundary quadrangles - goto a continous face if ( const QFace* f = link->GetContinuesFace( this )) - return f->GetLinkChain( *chLink, chain, pos, error ); + if ( f->_sides.size() == 4 ) + return f->GetLinkChain( *chLink, chain, pos, error ); } else { TChainLink chLink(link); // side face of prismatic mesh - visit all faces of iSide @@ -1682,7 +2367,7 @@ namespace { // Structures used by FixQuadraticElements() TLinkInSet link = links.find( _sides[iL] ); if ( link == linksEnd ) continue; if ( (*link)->MediumPos() > SMDS_TOP_FACE ) - continue; // We work on faces here, don't go into a volume + continue; // We work on faces here, don't go inside a solid // check link if ( link->IsBoundary() ) { @@ -1825,10 +2510,11 @@ namespace { // Structures used by FixQuadraticElements() // propagate to adjacent faces till limit step or boundary double len1 = thePrevLen + (theLink->MiddlePnt() - _sides[iL1]->MiddlePnt()).Modulus(); double len2 = thePrevLen + (theLink->MiddlePnt() - _sides[iL2]->MiddlePnt()).Modulus(); - gp_Vec linkDir1, linkDir2; + gp_Vec linkDir1(0,0,0); // initialize to avoid valgrind error ("Conditional jump...") + gp_Vec linkDir2(0,0,0); try { OCC_CATCH_SIGNALS; - if ( f1 ) + if ( f1 && theLink->MediumPos() <= (*link1)->MediumPos() ) len1 = f1->MoveByBoundary ( *link1, theRefVec, theLinks, theFaceHelper, len1, theStep-1, &linkDir1, theSign); else @@ -1839,7 +2525,7 @@ namespace { // Structures used by FixQuadraticElements() } try { OCC_CATCH_SIGNALS; - if ( f2 ) + if ( f2 && theLink->MediumPos() <= (*link2)->MediumPos() ) len2 = f2->MoveByBoundary ( *link2, theRefVec, theLinks, theFaceHelper, len2, theStep-1, &linkDir2, theSign); else @@ -1869,6 +2555,37 @@ namespace { // Structures used by FixQuadraticElements() return fullLen; } + //================================================================================ + /*! + * \brief Checks if the face is distorted due to bentLink + */ + //================================================================================ + + bool QFace::IsSpoiled(const QLink* bentLink ) const + { + // code is valid for convex faces only + gp_XYZ gc(0,0,0); + for ( TIDSortedNodeSet::const_iterator n = begin(); n!=end(); ++n) + gc += XYZ( *n ) / size(); + for (unsigned i = 0; i < _sides.size(); ++i ) + { + if ( _sides[i] == bentLink ) continue; + gp_Vec linkNorm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2())); + gp_Vec vecOut( gc, _sides[i]->MiddlePnt() ); + if ( linkNorm * vecOut < 0 ) + linkNorm.Reverse(); + double mag2 = linkNorm.SquareMagnitude(); + if ( mag2 > numeric_limits::min() ) + linkNorm /= sqrt( mag2 ); + gp_Vec vecBent ( _sides[i]->MiddlePnt(), bentLink->MediumPnt()); + gp_Vec vecStraight( _sides[i]->MiddlePnt(), bentLink->MiddlePnt()); + if ( vecBent * linkNorm > -0.1*vecStraight.Magnitude() ) + return true; + } + return false; + + } + //================================================================================ /*! * \brief Find pairs of continues faces @@ -1887,7 +2604,9 @@ namespace { // Structures used by FixQuadraticElements() if ( _faces.empty() ) return; - int iFaceCont = -1; + int iFaceCont = -1, nbBoundary = 0, iBoundary[2]={-1,-1}; + if ( _faces[0]->IsBoundary() ) + iBoundary[ nbBoundary++ ] = 0; for ( int iF = 1; iFaceCont < 0 && iF < _faces.size(); ++iF ) { // look for a face bounding none of volumes bound by _faces[0] @@ -1898,8 +2617,21 @@ namespace { // Structures used by FixQuadraticElements() _faces[iF]->_volumes[iV] == _faces[0]->_volumes[1]); if ( !sameVol ) iFaceCont = iF; + if ( _faces[iF]->IsBoundary() ) + iBoundary[ nbBoundary++ ] = iF; + } + // Set continues faces: arrange _faces to have + // _faces[0] continues to _faces[1] + // _faces[2] continues to _faces[3] + if ( nbBoundary == 2 ) // bnd faces are continues + { + if (( iBoundary[0] < 2 ) != ( iBoundary[1] < 2 )) + { + int iNear0 = iBoundary[0] < 2 ? 1-iBoundary[0] : 5-iBoundary[0]; + std::swap( _faces[ iBoundary[1] ], _faces[iNear0] ); + } } - if ( iFaceCont > 0 ) // continues faces found, set one by the other + else if ( iFaceCont > 0 ) // continues faces found { if ( iFaceCont != 1 ) std::swap( _faces[1], _faces[iFaceCont] ); @@ -1949,6 +2681,27 @@ namespace { // Structures used by FixQuadraticElements() if (_qfaces[1]) norm += _qfaces[1]->_normal; return norm; } + //================================================================================ + /*! + * \brief Test link curvature taking into account size of faces + */ + //================================================================================ + + bool TChainLink::IsStraight() const + { + bool isStraight = _qlink->IsStraight(); + if ( isStraight && _qfaces[0] && !_qfaces[1] ) + { + int i = _qfaces[0]->LinkIndex( _qlink ); + int iOpp = ( i + 2 ) % _qfaces[0]->_sides.size(); + gp_XYZ mid1 = _qlink->MiddlePnt(); + gp_XYZ mid2 = _qfaces[0]->_sides[ iOpp ]->MiddlePnt(); + double faceSize2 = (mid1-mid2).SquareModulus(); + isStraight = _qlink->_nodeMove.SquareMagnitude() < 1/10./10. * faceSize2; + } + return isStraight; + } + //================================================================================ /*! * \brief Move medium nodes of vertical links of pentahedrons adjacent by side faces @@ -1967,7 +2720,7 @@ namespace { // Structures used by FixQuadraticElements() bndLinks1.insert( lnk->_qlink ); else interLinks.insert( lnk->_qlink ); - isCurved = isCurved || !(*lnk)->IsStraight(); + isCurved = isCurved || !lnk->IsStraight(); } if ( !isCurved ) return; // no need to move @@ -2014,41 +2767,20 @@ namespace { // Structures used by FixQuadraticElements() TLinkSet linkSet( allLinks.begin(), allLinks.end()); TLinkInSet linkIt = linkSet.begin(), linksEnd = linkSet.end(); - // move in 2d if we are on geom face -// TopoDS_Face face; -// TopLoc_Location loc; -// SMESH_MesherHelper faceHelper( *helper.GetMesh()); -// while ( linkIt->IsBoundary()) ++linkIt; -// if ( linkIt == linksEnd ) return; -// if ( (*linkIt)->MediumPos() == SMDS_TOP_FACE ) { -// bool checkPos = true; -// TopoDS_Shape f = helper.GetSubShapeByNode( (*linkIt)->_mediumNode, helper.GetMeshDS() ); -// if ( !f.IsNull() && f.ShapeType() == TopAbs_FACE ) { -// face = TopoDS::Face( f ); -// helper.GetNodeUV( face, (*linkIt)->_mediumNode, 0, &checkPos); -// if (checkPos) -// face.Nullify(); -// else -// faceHelper.SetSubShape( face ); -// } -// } for ( linkIt = linkSet.begin(); linkIt != linksEnd; ++linkIt) { - if ( linkIt->IsBoundary() && !(*linkIt)->IsStraight() && linkIt->_qfaces[0]) + if ( linkIt->IsBoundary() && !linkIt->IsStraight() && linkIt->_qfaces[0]) { -// if ( !face.IsNull() ) { -// const SMDS_MeshNode* inFaceNode = -// faceHelper.GetNodeUVneedInFaceNode() ? linkIt->_qfaces[0]->GetNodeInFace() : 0; -// gp_XY uvm = helper.GetNodeUV( face, (*linkIt)->_mediumNode, inFaceNode ); -// gp_XY uv1 = helper.GetNodeUV( face, (*linkIt)->node1(), inFaceNode); -// gp_XY uv2 = helper.GetNodeUV( face, (*linkIt)->node2(), inFaceNode); -// gp_XY uvMove = uvm - helper.GetMiddleUV( BRep_Tool::Surface(face,loc), uv1, uv2); -// gp_Vec move( uvMove.X(), uvMove.Y(), 0 ); -// linkIt->_qfaces[0]->MoveByBoundary( *linkIt, move, linkSet, &faceHelper ); -// } -// else { - linkIt->_qfaces[0]->MoveByBoundary( *linkIt, (*linkIt)->_nodeMove, linkSet ); - //} + // move iff a boundary link is bent towards inside of a face (issue 0021084) + const QFace* face = linkIt->_qfaces[0]; + gp_XYZ pIn = ( face->_sides[0]->MiddlePnt() + + face->_sides[1]->MiddlePnt() + + face->_sides[2]->MiddlePnt() ) / 3.; + gp_XYZ insideDir( pIn - (*linkIt)->MiddlePnt()); + bool linkBentInside = ((*linkIt)->_nodeMove.Dot( insideDir ) > 0 ); + //if ( face->IsSpoiled( linkIt->_qlink )) + if ( linkBentInside ) + face->MoveByBoundary( *linkIt, (*linkIt)->_nodeMove, linkSet ); } } } @@ -2061,7 +2793,7 @@ namespace { // Structures used by FixQuadraticElements() enum TSplitTriaResult { _OK, _NO_CORNERS, _FEW_ROWS, _MANY_ROWS, _NO_SIDELINK, _BAD_MIDQUAD, _NOT_RECT, - _NO_MIDQUAD, _NO_UPTRIA, _BAD_SET_SIZE, _BAD_CORNER, _BAD_START, _NO_BOTLINK }; + _NO_MIDQUAD, _NO_UPTRIA, _BAD_SET_SIZE, _BAD_CORNER, _BAD_START, _NO_BOTLINK, _TWISTED_CHAIN }; TSplitTriaResult splitTrianglesIntoChains( TChain & allLinks, vector< TChain> & resultChains, @@ -2127,6 +2859,8 @@ namespace { // Structures used by FixQuadraticElements() const QFace* botTria = botLink->_qfaces[0]; // bottom triangle bound by botLink if ( !botTria ) { // the column ends + if ( botLink == startLink ) + return _TWISTED_CHAIN; // issue 0020951 linkSet.erase( botLink ); if ( iRow != rowChains.size() ) return _FEW_ROWS; // different nb of rows in columns @@ -2186,8 +2920,10 @@ namespace { // Structures used by FixQuadraticElements() // next bottom link ends at the new corner linkSet.erase( botLink ); botLink = upTria->GetLinkByNode( linkSet, (isCase2 ? *sideLink : *midQuadLink), corner ); - if ( botLink == linksEnd || botLink == (isCase2 ? midQuadLink : sideLink)) + if ( botLink == linksEnd || botLink == midQuadLink || botLink == sideLink) return _NO_BOTLINK; + if ( midQuadLink == startLink || sideLink == startLink ) + return _TWISTED_CHAIN; // issue 0020951 linkSet.erase( midQuadLink ); linkSet.erase( sideLink ); @@ -2231,7 +2967,7 @@ namespace { // Structures used by FixQuadraticElements() return _OK; } -} +} //namespace //======================================================================= /*! @@ -2244,37 +2980,55 @@ namespace { // Structures used by FixQuadraticElements() void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) { - // apply algorithm to solids or geom faces + // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion + if ( getenv("NO_FixQuadraticElements") ) + return; + + // 0. Apply algorithm to solids or geom faces // ---------------------------------------------- if ( myShape.IsNull() ) { if ( !myMesh->HasShapeToMesh() ) return; SetSubShape( myMesh->GetShapeToMesh() ); +#ifdef _DEBUG_ + int nbSolids = 0; + TopTools_IndexedMapOfShape solids; + TopExp::MapShapes(myShape,TopAbs_SOLID,solids); + nbSolids = solids.Extent(); +#endif TopTools_MapOfShape faces; // faces not in solid or in not meshed solid for ( TopExp_Explorer f(myShape,TopAbs_FACE,TopAbs_SOLID); f.More(); f.Next() ) { - faces.Add( f.Current() ); + faces.Add( f.Current() ); // not in solid } for ( TopExp_Explorer s(myShape,TopAbs_SOLID); s.More(); s.Next() ) { if ( myMesh->GetSubMesh( s.Current() )->IsEmpty() ) { // get faces of solid for ( TopExp_Explorer f( s.Current(), TopAbs_FACE); f.More(); f.Next() ) - faces.Add( f.Current() ); + faces.Add( f.Current() ); // in not meshed solid } else { // fix nodes in the solid and its faces +#ifdef _DEBUG_ + MSG("FIX SOLID " << nbSolids-- << " #" << GetMeshDS()->ShapeToIndex(s.Current())); +#endif SMESH_MesherHelper h(*myMesh); h.SetSubShape( s.Current() ); h.FixQuadraticElements(false); } } // fix nodes on geom faces +#ifdef _DEBUG_ + int nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; +#endif for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) { + MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key())); SMESH_MesherHelper h(*myMesh); h.SetSubShape( fIt.Key() ); h.FixQuadraticElements(true); } + //perf_print_all_meters(1); return; } - // Find out type of elements and get iterator on them + // 1. Find out type of elements and get iterator on them // --------------------------------------------------- SMDS_ElemIteratorPtr elemIt; @@ -2293,7 +3047,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) if ( !elemIt || !elemIt->more() || elemType < SMDSAbs_Face ) return; - // Fill in auxiliary data structures + // 2. Fill in auxiliary data structures // ---------------------------------- set< QLink > links; @@ -2302,21 +3056,25 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) set< QFace >::iterator pFace; bool isCurved = false; - bool hasRectFaces = false; - set nbElemNodeSet; + //bool hasRectFaces = false; + //set nbElemNodeSet; + SMDS_VolumeTool volTool; + + TIDSortedNodeSet apexOfPyramid; + const int apexIndex = 4; if ( elemType == SMDSAbs_Volume ) { - SMDS_VolumeTool volTool; while ( elemIt->more() ) // loop on volumes { const SMDS_MeshElement* vol = elemIt->next(); if ( !vol->IsQuadratic() || !volTool.Set( vol )) - return; //continue; + return; + double volMinSize2 = -1.; for ( int iF = 0; iF < volTool.NbFaces(); ++iF ) // loop on faces of volume { int nbN = volTool.NbFaceNodes( iF ); - nbElemNodeSet.insert( nbN ); + //nbElemNodeSet.insert( nbN ); const SMDS_MeshNode** faceNodes = volTool.GetFaceNodes( iF ); vector< const QLink* > faceLinks( nbN/2 ); for ( int iN = 0; iN < nbN; iN += 2 ) // loop on links of a face @@ -2325,19 +3083,26 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) QLink link( faceNodes[iN], faceNodes[iN+2], faceNodes[iN+1] ); pLink = links.insert( link ).first; faceLinks[ iN/2 ] = & *pLink; - if ( !isCurved ) - isCurved = !link.IsStraight(); - if ( link.MediumPos() == SMDS_TOP_3DSPACE && !link.IsStraight() ) - return; // already fixed + + if ( link.MediumPos() == SMDS_TOP_3DSPACE ) + { + if ( !link.IsStraight() ) + return; // already fixed + } + else if ( !isCurved ) + { + if ( volMinSize2 < 0 ) volMinSize2 = volTool.MinLinearSize2(); + isCurved = !isStraightLink( volMinSize2, link._nodeMove.SquareMagnitude() ); + } } // store QFace pFace = faces.insert( QFace( faceLinks )).first; if ( pFace->NbVolumes() == 0 ) pFace->AddSelfToLinks(); pFace->SetVolume( vol ); - hasRectFaces = hasRectFaces || - ( volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_HEXA || - volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_PENTA ); +// hasRectFaces = hasRectFaces || +// ( volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_HEXA || +// volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_PENTA ); #ifdef _DEBUG_ if ( nbN == 6 ) pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],faceNodes[4]); @@ -2346,6 +3111,9 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) faceNodes[4],faceNodes[6] ); #endif } + // collect pyramid apexes for further correction + if ( vol->NbCornerNodes() == 5 ) + apexOfPyramid.insert( vol->GetNode( apexIndex )); } set< QLink >::iterator pLink = links.begin(); for ( ; pLink != links.end(); ++pLink ) @@ -2358,7 +3126,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) const SMDS_MeshElement* face = elemIt->next(); if ( !face->IsQuadratic() ) continue; - nbElemNodeSet.insert( face->NbNodes() ); + //nbElemNodeSet.insert( face->NbNodes() ); int nbN = face->NbNodes()/2; vector< const QLink* > faceLinks( nbN ); for ( int iN = 0; iN < nbN; ++iN ) // loop on links of a face @@ -2373,27 +3141,28 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) // store QFace pFace = faces.insert( QFace( faceLinks )).first; pFace->AddSelfToLinks(); - hasRectFaces = ( hasRectFaces || nbN == 4 ); + //hasRectFaces = ( hasRectFaces || nbN == 4 ); } } if ( !isCurved ) return; // no curved edges of faces - // Compute displacement of medium nodes - // ------------------------------------- + // 3. Compute displacement of medium nodes + // --------------------------------------- - // two loops on faces: the first is to treat boundary links, the second is for internal ones + // two loops on QFaces: the first is to treat boundary links, the second is for internal ones TopLoc_Location loc; // not treat boundary of volumic submesh int isInside = ( elemType == SMDSAbs_Volume && volumeOnly ) ? 1 : 0; for ( ; isInside < 2; ++isInside ) { MSG( "--------------- LOOP (inside=" << isInside << ") ------------------"); SMDS_TypeOfPosition pos = isInside ? SMDS_TOP_3DSPACE : SMDS_TOP_FACE; + SMDS_TypeOfPosition bndPos = isInside ? SMDS_TOP_FACE : SMDS_TOP_EDGE; for ( pFace = faces.begin(); pFace != faces.end(); ++pFace ) { if ( bool(isInside) == pFace->IsBoundary() ) continue; - for ( int dir = 0; dir < 2; ++dir ) // 2 directions of propagation from quadrangle + for ( int dir = 0; dir < 2; ++dir ) // 2 directions of propagation from the quadrangle { MSG( "CHAIN"); // make chain of links connected via continues faces @@ -2426,8 +3195,13 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) { TChain& chain = chains[iC]; if ( chain.empty() ) continue; - if ( chain.front()->IsStraight() && chain.back()->IsStraight() ) { - MSG("3D straight"); + if ( chain.front().IsStraight() && chain.back().IsStraight() ) { + MSG("3D straight - ignore"); + continue; + } + if ( chain.front()->MediumPos() > bndPos || + chain.back() ->MediumPos() > bndPos ) { + MSG("Internal chain - ignore"); continue; } // mesure chain length and compute link position along the chain @@ -2456,14 +3230,17 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) TopoDS_Face face; bool checkUV = true; - if ( !isInside ) { - // compute node displacement of end links in parametric space of face - const SMDS_MeshNode* nodeOnFace = (*(++chain.begin()))->_mediumNode; + if ( !isInside ) + { + // compute node displacement of end links of chain in parametric space of face + TChainLink& linkOnFace = *(++chain.begin()); + const SMDS_MeshNode* nodeOnFace = linkOnFace->_mediumNode; TopoDS_Shape f = GetSubShapeByNode( nodeOnFace, GetMeshDS() ); if ( !f.IsNull() && f.ShapeType() == TopAbs_FACE ) { face = TopoDS::Face( f ); Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc); + bool isStraight[2]; for ( int is1 = 0; is1 < 2; ++is1 ) // move0 or move1 { TChainLink& link = is1 ? chain.back() : chain.front(); @@ -2474,12 +3251,26 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) // uvMove = uvm - uv12 gp_XY uvMove = applyIn2D(surf, uvm, uv12, gp_XY_Subtracted, /*inPeriod=*/false); ( is1 ? move1 : move0 ).SetCoord( uvMove.X(), uvMove.Y(), 0 ); + if ( !is1 ) // correct nodeOnFace for move1 (issue 0020919) + nodeOnFace = (*(++chain.rbegin()))->_mediumNode; + isStraight[is1] = isStraightLink( (uv2-uv1).SquareModulus(), + 10 * uvMove.SquareModulus()); } - if ( move0.SquareMagnitude() < straightTol2 && - move1.SquareMagnitude() < straightTol2 ) { - MSG("2D straight"); + if ( isStraight[0] && isStraight[1] ) { + MSG("2D straight - ignore"); continue; // straight - no need to move nodes of internal links } + + // check if a chain is already fixed + gp_XY uvm = GetNodeUV( face, linkOnFace->_mediumNode, 0, &checkUV); + gp_XY uv1 = GetNodeUV( face, linkOnFace->node1(), nodeOnFace, &checkUV); + gp_XY uv2 = GetNodeUV( face, linkOnFace->node2(), nodeOnFace, &checkUV); + gp_XY uv12 = GetMiddleUV( surf, uv1, uv2); + if (( uvm - uv12 ).SquareModulus() > 1e-10 ) + { + MSG("Already fixed - ignore"); + continue; + } } } gp_Trsf trsf; @@ -2529,7 +3320,8 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) { gp_XY uv0 = GetNodeUV( face, (*link0)->_mediumNode, 0, &checkUV); gp_XY uv2 = GetNodeUV( face, (*link2)->_mediumNode, 0, &checkUV); - MSG( "uv0: "< vols( 100 ); +// vector volSize( 100 ); +// int nbVols; +// bool ok; for ( pLink = links.begin(); pLink != links.end(); ++pLink ) { if ( pLink->IsMoved() ) { - //gp_Pnt p = pLink->MediumPnt() + pLink->Move(); gp_Pnt p = pLink->MiddlePnt() + pLink->Move(); GetMeshDS()->MoveNode( pLink->_mediumNode, p.X(), p.Y(), p.Z()); + // +// gp_Pnt pNew = pLink->MiddlePnt() + pLink->Move(); +// if ( pLink->MediumPos() != SMDS_TOP_3DSPACE ) +// { +// // avoid making distorted volumes near boundary +// SMDS_ElemIteratorPtr volIt = +// (*pLink)._mediumNode->GetInverseElementIterator( SMDSAbs_Volume ); +// for ( nbVols = 0; volIt->more() && volTool.Set( volIt->next() ); ++nbVols ) +// { +// vols [ nbVols ] = volTool.Element(); +// volSize[ nbVols ] = volTool.GetSize(); +// } +// gp_Pnt pOld = pLink->MediumPnt(); +// const_cast( pLink->_mediumNode )->setXYZ( pNew.X(), pNew.Y(), pNew.Z() ); +// ok = true; +// while ( nbVols-- && ok ) +// { +// volTool.Set( vols[ nbVols ]); +// ok = ( volSize[ nbVols ] * volTool.GetSize() > 1e-20 ); +// } +// if ( !ok ) +// { +// const_cast( pLink->_mediumNode )->setXYZ( pOld.X(), pOld.Y(), pOld.Z() ); +// MSG( "Do NOT move \t" << pLink->_mediumNode->GetID() +// << " because of distortion of volume " << vols[ nbVols+1 ]->GetID()); +// continue; +// } +// } +// GetMeshDS()->MoveNode( pLink->_mediumNode, pNew.X(), pNew.Y(), pNew.Z() ); } } -} -//======================================================================= -/*! - * \brief Iterator on ancestors of the given type - */ -//======================================================================= + //return; -struct TAncestorsIterator : public SMDS_Iterator -{ - TopTools_ListIteratorOfListOfShape _ancIter; - TopAbs_ShapeEnum _type; - TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type) - : _ancIter( ancestors ), _type( type ) - { - if ( _ancIter.More() && _ancIter.Value().ShapeType() != _type ) next(); - } - virtual bool more() - { - return _ancIter.More(); - } - virtual const TopoDS_Shape* next() + // issue 0020982 + // Move the apex of pyramid together with the most curved link + + TIDSortedNodeSet::iterator apexIt = apexOfPyramid.begin(); + for ( ; apexIt != apexOfPyramid.end(); ++apexIt ) { - const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0; - if ( _ancIter.More() ) - for ( _ancIter.Next(); _ancIter.More(); _ancIter.Next()) - if ( _ancIter.Value().ShapeType() == _type ) - break; - return s; - } -}; + SMESH_TNodeXYZ apex = *apexIt; -//======================================================================= -/*! - * \brief Return iterator on ancestors of the given type - */ -//======================================================================= + gp_Vec maxMove( 0,0,0 ); + double maxMoveSize2 = 0; -PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape, - const SMESH_Mesh& mesh, - TopAbs_ShapeEnum ancestorType) -{ - return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType)); + // shift of node index to get medium nodes between the base nodes + const int base2MediumShift = 5; + + // find maximal movement of medium node + SMDS_ElemIteratorPtr volIt = apex._node->GetInverseElementIterator( SMDSAbs_Volume ); + vector< const SMDS_MeshElement* > pyramids; + while ( volIt->more() ) + { + const SMDS_MeshElement* pyram = volIt->next(); + if ( pyram->GetEntityType() != SMDSEntity_Quad_Pyramid ) continue; + pyramids.push_back( pyram ); + + for ( int iBase = 0; iBase < apexIndex; ++iBase ) + { + SMESH_TNodeXYZ medium = pyram->GetNode( iBase + base2MediumShift ); + if ( medium._node->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE ) + { + SMESH_TNodeXYZ n1 = pyram->GetNode( iBase ); + SMESH_TNodeXYZ n2 = pyram->GetNode( ( iBase+1 ) % 4 ); + gp_Pnt middle = 0.5 * ( n1 + n2 ); + gp_Vec move( middle, medium ); + double moveSize2 = move.SquareMagnitude(); + if ( moveSize2 > maxMoveSize2 ) + maxMove = move, maxMoveSize2 = moveSize2; + } + } + } + + // move the apex + if ( maxMoveSize2 > 1e-20 ) + { + apex += maxMove.XYZ(); + GetMeshDS()->MoveNode( apex._node, apex.X(), apex.Y(), apex.Z()); + + // move medium nodes neighboring the apex to the middle + const int base2MediumShift_2 = 9; + for ( unsigned i = 0; i < pyramids.size(); ++i ) + for ( int iBase = 0; iBase < apexIndex; ++iBase ) + { + SMESH_TNodeXYZ base = pyramids[i]->GetNode( iBase ); + const SMDS_MeshNode* medium = pyramids[i]->GetNode( iBase + base2MediumShift_2 ); + gp_XYZ middle = 0.5 * ( apex + base ); + GetMeshDS()->MoveNode( medium, middle.X(), middle.Y(), middle.Z()); + } + } + } } +