X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;ds=sidebyside;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=4f732effc61f0fa0f569a847094c4184d174144a;hb=4cd2499bddcd3da3ec8900fe825bc98669b789b5;hp=bcdd2096a0e5f22528ef9bec48357005c07124e7;hpb=8fa039a796957b302d86d90b22afc0a998573f83;p=modules%2Fsmesh.git diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index bcdd2096a..4f732effc 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 @@ -26,19 +26,24 @@ // #include "SMESH_MesherHelper.hxx" -#include "SMDS_FacePosition.hxx" #include "SMDS_EdgePosition.hxx" +#include "SMDS_FaceOfNodes.hxx" +#include "SMDS_FacePosition.hxx" +#include "SMDS_IteratorOnIterators.hxx" #include "SMDS_VolumeTool.hxx" +#include "SMESH_Block.hxx" +#include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" +#include #include #include -#include #include #include #include #include #include +#include #include #include #include @@ -76,7 +81,8 @@ namespace { //================================================================================ SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh) - : myParIndex(0), myMesh(&theMesh), myShapeID(0), myCreateQuadratic(false) + : myParIndex(0), myMesh(&theMesh), myShapeID(0), myCreateQuadratic(false), + myFixNodeParameters(false) { myPar1[0] = myPar2[0] = myPar1[1] = myPar2[1] = 0; mySetElemOnShape = ( ! myMesh->HasShapeToMesh() ); @@ -112,45 +118,67 @@ 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(); myDegenShapeIds.clear(); TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE ); + if ( aSh.ShapeType()==TopAbs_COMPOUND ) + { + TopoDS_Iterator subIt( aSh ); + if ( subIt.More() ) + subType = ( subIt.Value().ShapeType()==TopAbs_FACE ) ? TopAbs_EDGE : TopAbs_FACE; + } SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge ); + 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,14 +198,14 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) //======================================================================= //function : SetSubShape -//purpose : Set geomerty to make elements on +//purpose : Set geometry to make elements on //======================================================================= 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() ); @@ -185,7 +213,7 @@ void SMESH_MesherHelper::SetSubShape(const int aShID) //======================================================================= //function : SetSubShape -//purpose : Set geomerty to create elements on +//purpose : Set geometry to create elements on //======================================================================= void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) @@ -209,9 +237,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 @@ -223,13 +258,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 ); @@ -247,6 +282,18 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() )); } } + if ( !myDegenShapeIds.empty() && !myParIndex ) { + if ( surface->IsUPeriodic() || surface->IsUClosed() ) { + myParIndex |= U_periodic; + myPar1[0] = surf.FirstUParameter(); + myPar2[0] = surf.LastUParameter(); + } + else if ( surface->IsVPeriodic() || surface->IsVClosed() ) { + myParIndex |= V_periodic; + myPar1[1] = surf.FirstVParameter(); + myPar2[1] = surf.LastVParameter(); + } + } } } } @@ -255,7 +302,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() //======================================================================= @@ -291,9 +338,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->getshapeId(); + int shapeID = node ? node->getshapeId() : 0; if ( 0 < shapeID && shapeID <= meshDS->MaxShapeIndex() ) return meshDS->IndexToShape( shapeID ); else @@ -315,6 +362,75 @@ 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 @@ -340,6 +456,20 @@ void SMESH_MesherHelper::setPosOnShapeValidity(int shapeID, bool ok ) const ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok)); } +//======================================================================= +//function : ToFixNodeParameters +//purpose : Enables fixing node parameters on EDGEs and FACEs in +// GetNodeU(...,check=true), GetNodeUV(...,check=true), CheckNodeUV() and +// CheckNodeU() in case if a node lies on a shape set via SetSubShape(). +// Default is False +//======================================================================= + +void SMESH_MesherHelper::ToFixNodeParameters(bool toFix) +{ + myFixNodeParameters = toFix; +} + + //======================================================================= //function : GetUVOnSeam //purpose : Select UV on either of 2 pcurves of a seam edge, closest to the given UV @@ -405,7 +535,7 @@ 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(), 10*MaxTolerance( F ),/*force=*/ !validU ); @@ -484,6 +614,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; @@ -500,24 +634,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 { int shapeID = n->getshapeId(); - if ( force || toCheckPosOnShape( shapeID )) + bool infinit = ( Precision::IsInfinite( uv.X() ) || Precision::IsInfinite( uv.Y() )); + if ( force || toCheckPosOnShape( shapeID ) || infinit ) { - double toldis = tol; - double tolmin = 1.e-7*myMesh->GetMeshDS()->getMaxDim(); // nodes coordinates are stored in float format - if (toldis < tolmin) toldis = tolmin; // 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() )) > toldis ) + 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 ); @@ -529,13 +667,20 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, Quantity_Parameter U,V; projector.LowerDistanceParameters(U,V); uv.SetCoord( U,V ); - if ( nodePnt.Distance( surface->Value( U, V )) > toldis ) + 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; } // store the fixed UV on the face - if ( myShape.IsSame(F) && shapeID == myShapeID ) + if ( myShape.IsSame(F) && shapeID == myShapeID && myFixNodeParameters ) const_cast(n)->SetPosition ( SMDS_PositionPtr( new SMDS_FacePosition( U, V ))); } @@ -628,7 +773,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 ); } //======================================================================= @@ -689,15 +841,11 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, double& u, const double tol, const bool force, - double* distance) const + double distXYZ[4]) const { int shapeID = n->getshapeId(); if ( force || toCheckPosOnShape( shapeID )) { - double toldis = tol; - double tolmin = 1.e-7*myMesh->GetMeshDS()->getMaxDim(); // nodes coordinates are stored in float format - if (toldis < tolmin) toldis = tolmin; - // 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 @@ -710,11 +858,16 @@ 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() ); - double dist = nodePnt.Distance( curve->Value( u )); - if ( distance ) *distance = dist; - if ( dist > toldis ) + 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 @@ -736,16 +889,21 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, } Quantity_Parameter U = projector->LowerDistanceParameter(); u = double( U ); - dist = nodePnt.Distance( curve->Value( U )); - if ( distance ) *distance = dist; - if ( dist > toldis ) + 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 " << nodePnt.Distance(curve->Value( U )) << " " << toldis); + MESSAGE("distance " << dist << " " << tol ); return false; } // store the fixed U on the edge - if ( myShape.IsSame(E) && shapeID == myShapeID ) + if ( myShape.IsSame(E) && shapeID == myShapeID && myFixNodeParameters ) const_cast(n)->SetPosition ( SMDS_PositionPtr( new SMDS_EdgePosition( U ))); } @@ -772,6 +930,85 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, 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 @@ -800,72 +1037,51 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, // 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(); - TopoDS_Edge E; double u [2]; TopoDS_Face F; gp_XY uv[2]; bool uvOK[2] = { false, false }; - if( myShape.IsNull() ) - { - if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) { - faceID = n1->getshapeId(); - } - else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) { - faceID = n2->getshapeId(); - } + pair pos = GetMediumPos( n1, n2 ); - if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) { - edgeID = n1->getshapeId(); - } - if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) { - edgeID = n2->getshapeId(); - } - } // get positions of the given nodes on shapes - TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType(); - if ( faceID>0 || shapeType == TopAbs_FACE) + 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 ) { + 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); - - if( myShape.IsNull() ) - E = TopoDS::Edge(meshDS->IndexToShape(edgeID)); - else { - E = TopoDS::Edge(myShape); - edgeID = myShapeID; + n1->getshapeId() != n2->getshapeId() ) + { + // issue 0021006 + 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]); } - 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( n1->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( n2->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); @@ -946,7 +1162,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_ // 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(); + double u = 0, distMiddleProj = Precision::Infinite(), distXYZ[4]; int iOkEdge = 0; TopoDS_Edge edges[2]; for ( int is2nd = 0; is2nd < 2; ++is2nd ) @@ -960,11 +1176,11 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_ // 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 ), foundDist = node2MiddleDist; - CheckNodeU( edge, n12, foundU, 2*BRep_Tool::Tolerance(edge), /*force=*/true, &foundDist ); - if ( foundDist < node2MiddleDist ) + double foundU = GetNodeU( edge, n ); + CheckNodeU( edge, n12, foundU, 2*BRep_Tool::Tolerance(edge), /*force=*/true, distXYZ ); + if ( distXYZ[0] < node2MiddleDist ) { - distMiddleProj = foundDist; + distMiddleProj = distXYZ[0]; u = foundU; iOkEdge = is2nd; } @@ -1007,7 +1223,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_ //purpose : Creates a node //======================================================================= -SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID) +SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID, + double u, double v) { SMESHDS_Mesh * meshDS = GetMeshDS(); SMDS_MeshNode* node = 0; @@ -1017,11 +1234,11 @@ SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID) node = meshDS->AddNode( x, y, z ); if ( mySetElemOnShape && myShapeID > 0 ) { switch ( myShape.ShapeType() ) { - case TopAbs_SOLID: meshDS->SetNodeInVolume( node, myShapeID); break; - case TopAbs_SHELL: meshDS->SetNodeInVolume( node, myShapeID); break; - case TopAbs_FACE: meshDS->SetNodeOnFace( node, myShapeID); break; - case TopAbs_EDGE: meshDS->SetNodeOnEdge( node, myShapeID); break; - case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID); break; + case TopAbs_SOLID: meshDS->SetNodeInVolume( node, myShapeID); break; + case TopAbs_SHELL: meshDS->SetNodeInVolume( node, myShapeID); break; + case TopAbs_FACE: meshDS->SetNodeOnFace( node, myShapeID, u, v); break; + case TopAbs_EDGE: meshDS->SetNodeOnEdge( node, myShapeID, u); break; + case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID); break; default: ; } } @@ -1180,7 +1397,7 @@ SMDS_MeshFace* SMESH_MesherHelper::AddPolygonalFace (const vectorAddVolumeWithID(n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,id); + else + elem = meshDS->AddVolume(n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12); + if ( mySetElemOnShape && myShapeID > 0 ) + meshDS->SetMeshElementOnShape( elem, myShapeID ); + return elem; +} + //======================================================================= //function : AddPolyhedralVolume //purpose : Creates polyhedron. In quadratic mesh, adds medium nodes @@ -1442,6 +1690,43 @@ SMESH_MesherHelper::AddPolyhedralVolume (const std::vector return elem; } +namespace +{ + //================================================================================ + /*! + * \brief Check if a node belongs to any face of sub-mesh + */ + //================================================================================ + + bool isNodeInSubMesh( const SMDS_MeshNode* n, const SMESHDS_SubMesh* sm ) + { + SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face ); + while ( fIt->more() ) + if ( sm->Contains( fIt->next() )) + return true; + return false; + } +} + +//======================================================================= +//function : IsSameElemGeometry +//purpose : Returns true if all elements of a sub-mesh are of same shape +//======================================================================= + +bool SMESH_MesherHelper::IsSameElemGeometry(const SMESHDS_SubMesh* smDS, + SMDSAbs_GeometryType shape, + const bool nullSubMeshRes) +{ + if ( !smDS ) return nullSubMeshRes; + + SMDS_ElemIteratorPtr elemIt = smDS->GetElements(); + while ( elemIt->more() ) + if ( elemIt->next()->GetGeomType() != shape ) + return false; + + return true; +} + //======================================================================= //function : LoadNodeColumns //purpose : Load nodes bound to face into a map of node columns @@ -1450,68 +1735,368 @@ SMESH_MesherHelper::AddPolyhedralVolume (const std::vector 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 ); - if ( !faceSubMesh || faceSubMesh->NbElements() == 0 ) - return false; + return LoadNodeColumns(theParam2ColumnMap, + theFace, + std::list(1,theBaseEdge), + theMesh, + theProxyMesh); +} + +//======================================================================= +//function : LoadNodeColumns +//purpose : Load nodes bound to face into a map of node columns +//======================================================================= - // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them +bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap, + const TopoDS_Face& theFace, + const std::list& theBaseSide, + SMESHDS_Mesh* theMesh, + SMESH_ProxyMesh* theProxyMesh) +{ + // get a right sub-mesh of theFace - map< double, const SMDS_MeshNode*> sortedBaseNodes; - if ( !SMESH_Algo::GetSortedNodesOnEdge( theMesh, theBaseEdge,/*noMedium=*/true, sortedBaseNodes) - || sortedBaseNodes.size() < 2 ) + 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; - 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++ ) + if ( theParam2ColumnMap.empty() ) { - double par = ( u_n->first - f ) / range; - vector& nCol = theParam2ColumnMap[ par ]; - nCol.resize( nbRows ); - nCol[0] = u_n->second; + // get data of edges for normalization of params + + 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 ); + } + } + + // 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 ) + { + map< double, const SMDS_MeshNode*> sortedBaseNN; + SMESH_Algo::GetSortedNodesOnEdge( theMesh, *edge,/*noMedium=*/true, sortedBaseNN); + if ( sortedBaseNN.empty() ) continue; + + map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNN.begin(); + if ( theProxyMesh ) // from sortedBaseNN remove nodes not shared by faces of faceSubMesh + { + const SMDS_MeshNode* n1 = sortedBaseNN.begin()->second; + const SMDS_MeshNode* n2 = sortedBaseNN.rbegin()->second; + bool allNodesAreProxy = ( n1 != theProxyMesh->GetProxyNode( n1 ) && + n2 != theProxyMesh->GetProxyNode( n2 )); + if ( allNodesAreProxy ) + for ( u_n = sortedBaseNN.begin(); u_n != sortedBaseNN.end(); u_n++ ) + u_n->second = theProxyMesh->GetProxyNode( u_n->second ); + + if ( u_n = sortedBaseNN.begin(), !isNodeInSubMesh( u_n->second, faceSubMesh )) + { + while ( ++u_n != sortedBaseNN.end() && !isNodeInSubMesh( u_n->second, faceSubMesh )); + sortedBaseNN.erase( sortedBaseNN.begin(), u_n ); + } + else if ( u_n = --sortedBaseNN.end(), !isNodeInSubMesh( u_n->second, faceSubMesh )) + { + while ( u_n != sortedBaseNN.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh )); + sortedBaseNN.erase( ++u_n, sortedBaseNN.end() ); + } + if ( sortedBaseNN.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; + for ( u_n = sortedBaseNN.begin(); u_n != sortedBaseNN.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 ); + } + } + if ( theParam2ColumnMap.empty() ) + return false; } + // nb rows of nodes + int prevNbRows = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here + int expectedNbRows = faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 ); // to be added + // 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++; + 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( prevNbRows + expectedNbRows ); + nCol2.resize( prevNbRows + expectedNbRows ); - int i1, i2, iRow = 0; - const SMDS_MeshNode *n1 = nCol1[0], *n2 = nCol2[0]; + int i1, i2, foundNbRows = 0; + const SMDS_MeshNode *n1 = nCol1[ prevNbRows-1 ]; + const SMDS_MeshNode *n2 = nCol2[ prevNbRows-1 ]; // find face sharing node n1 and n2 and belonging to faceSubMesh while ( const SMDS_MeshElement* face = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2)) { if ( faceSubMesh->Contains( face )) { - int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes(); + int nbNodes = face->NbCornerNodes(); if ( nbNodes != 4 ) return false; + if ( foundNbRows + 1 > expectedNbRows ) + return false; n1 = face->GetNode( (i2+2) % 4 ); // opposite corner of quadrangle face n2 = face->GetNode( (i1+2) % 4 ); - if ( ++iRow >= nbRows ) - return false; - nCol1[ iRow ] = n1; - nCol2[ iRow ] = n2; - avoidSet.clear(); + nCol1[ prevNbRows + foundNbRows] = n1; + nCol2[ prevNbRows + foundNbRows] = n2; + ++foundNbRows; } avoidSet.insert( face ); } - if ( iRow + 1 < nbRows ) // compact if necessary - nCol1.resize( iRow + 1 ), nCol2.resize( iRow + 1 ); + if ( foundNbRows != expectedNbRows ) + return false; + avoidSet.clear(); } + return ( theParam2ColumnMap.size() > 1 && + theParam2ColumnMap.begin()->second.size() == prevNbRows + expectedNbRows ); +} + +namespace +{ + //================================================================================ + /*! + * \brief Return true if a node is at a corner of a 2D structured mesh of FACE + */ + //================================================================================ + + bool isCornerOfStructure( const SMDS_MeshNode* n, + const SMESHDS_SubMesh* faceSM ) + { + int nbFacesInSM = 0; + if ( n ) { + SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face ); + while ( fIt->more() ) + nbFacesInSM += faceSM->Contains( fIt->next() ); + } + return ( nbFacesInSM == 1 ); + } +} + +//======================================================================= +//function : IsStructured +//purpose : Return true if 2D mesh on FACE is structured +//======================================================================= + +bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM ) +{ + SMESHDS_SubMesh* fSM = faceSM->GetSubMeshDS(); + if ( !fSM || fSM->NbElements() == 0 ) + return false; + + list< TopoDS_Edge > edges; + list< int > nbEdgesInWires; + int nbWires = SMESH_Block::GetOrderedEdges( TopoDS::Face( faceSM->GetSubShape() ), + edges, nbEdgesInWires ); + if ( nbWires != 1 ) + return false; + + // algo: find corners of a structure and then analyze nb of faces and + // length of structure sides + + SMESHDS_Mesh* meshDS = faceSM->GetFather()->GetMeshDS(); + + // rotate edges to get the first node being at corner + // (in principle it's not necessary but so far none SALOME algo can make + // such a structured mesh that all corner nodes are not on VERTEXes) + bool isCorner = false; + int nbRemainEdges = nbEdgesInWires.front(); + do { + TopoDS_Vertex V = IthVertex( 0, edges.front() ); + isCorner = isCornerOfStructure( SMESH_Algo::VertexNode( V, meshDS ), fSM); + if ( !isCorner ) { + edges.splice( edges.end(), edges, edges.begin() ); + --nbRemainEdges; + } + } + while ( !isCorner && nbRemainEdges > 0 ); + + if ( !isCorner ) + return false; + + // get all nodes from EDGEs + list< const SMDS_MeshNode* > nodes; + list< TopoDS_Edge >::iterator edge = edges.begin(); + for ( ; edge != edges.end(); ++edge ) + { + map< double, const SMDS_MeshNode* > u2Nodes; + if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, *edge, + /*skipMedium=*/true, u2Nodes )) + return false; + + list< const SMDS_MeshNode* > edgeNodes; + map< double, const SMDS_MeshNode* >::iterator u2n = u2Nodes.begin(); + if ( !nodes.empty() && nodes.back() == u2n->second ) + ++u2n; + map< double, const SMDS_MeshNode* >::iterator u2nEnd = --u2Nodes.end(); + if ( nodes.empty() || nodes.back() != u2nEnd->second ) + ++u2nEnd; + for ( ; u2n != u2nEnd; ++u2n ) + edgeNodes.push_back( u2n->second ); + + if ( edge->Orientation() == TopAbs_REVERSED ) + edgeNodes.reverse(); + nodes.splice( nodes.end(), edgeNodes, edgeNodes.begin(), edgeNodes.end() ); + } + + // get length of structured sides + vector nbEdgesInSide; + int nbEdges = 0; + list< const SMDS_MeshNode* >::iterator n = ++nodes.begin(); + for ( ; n != nodes.end(); ++n ) + { + ++nbEdges; + if ( isCornerOfStructure( *n, fSM )) { + nbEdgesInSide.push_back( nbEdges ); + nbEdges = 0; + } + } + + // checks + if ( nbEdgesInSide.size() != 4 ) + return false; + if ( nbEdgesInSide[0] != nbEdgesInSide[2] ) + return false; + if ( nbEdgesInSide[1] != nbEdgesInSide[3] ) + return false; + if ( nbEdgesInSide[0] * nbEdgesInSide[1] != fSM->NbElements() ) + return false; + return true; } +//================================================================================ +/*! + * \brief Find out elements orientation on a geometrical face + * \param theFace - The face correctly oriented in the shape being meshed + * \retval bool - true if the face normal and the normal of first element + * in the correspoding submesh point in different directions + */ +//================================================================================ + +bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace) +{ + if ( theFace.IsNull() ) + return false; + + // find out orientation of a meshed face + int faceID = GetMeshDS()->ShapeToIndex( theFace ); + TopoDS_Shape aMeshedFace = GetMeshDS()->IndexToShape( faceID ); + bool isReversed = ( theFace.Orientation() != aMeshedFace.Orientation() ); + + const SMESHDS_SubMesh * aSubMeshDSFace = GetMeshDS()->MeshElements( faceID ); + if ( !aSubMeshDSFace ) + return isReversed; + + // find an element with a good normal + gp_Vec Ne; + bool normalOK = false; + gp_XY uv; + SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); + while ( !normalOK && iteratorElem->more() ) // loop on elements on theFace + { + const SMDS_MeshElement* elem = iteratorElem->next(); + if ( elem && elem->NbCornerNodes() > 2 ) + { + SMESH_TNodeXYZ nPnt[3]; + SMDS_ElemIteratorPtr nodesIt = elem->nodesIterator(); + for ( int iN = 0; nodesIt->more() && iN < 3; ++iN) // loop on nodes + nPnt[ iN ] = nodesIt->next(); + + // compute normal + gp_Vec v01( nPnt[0], nPnt[1] ), v02( nPnt[0], nPnt[2] ); + if ( v01.SquareMagnitude() > RealSmall() && + v02.SquareMagnitude() > RealSmall() ) + { + Ne = v01 ^ v02; + if (( normalOK = ( Ne.SquareMagnitude() > RealSmall() ))) + uv = GetNodeUV( theFace, nPnt[0]._node, nPnt[2]._node, &normalOK ); + } + } + } + if ( !normalOK ) + return isReversed; + + // face normal at node position + TopLoc_Location loc; + Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace, loc ); + // if ( surf.IsNull() || surf->Continuity() < GeomAbs_C1 ) + // some surfaces not detected as GeomAbs_C1 are nevertheless correct for meshing + if ( surf.IsNull() || surf->Continuity() < GeomAbs_C0 ) + { + if (!surf.IsNull()) + MESSAGE("surf->Continuity() < GeomAbs_C1 " << (surf->Continuity() < GeomAbs_C1)); + return isReversed; + } + gp_Vec d1u, d1v; gp_Pnt p; + surf->D1( uv.X(), uv.Y(), p, d1u, d1v ); + gp_Vec Nf = (d1u ^ d1v).Transformed( loc ); + + if ( theFace.Orientation() == TopAbs_REVERSED ) + Nf.Reverse(); + + return Ne * Nf < 0.; +} + +//======================================================================= +//function : Count +//purpose : Count nb of sub-shapes +//======================================================================= + +int SMESH_MesherHelper::Count(const TopoDS_Shape& shape, + const TopAbs_ShapeEnum type, + const bool ignoreSame) +{ + if ( ignoreSame ) { + TopTools_IndexedMapOfShape map; + TopExp::MapShapes( shape, type, map ); + return map.Extent(); + } + else { + int nb = 0; + for ( TopExp_Explorer exp( shape, type ); exp.More(); exp.Next() ) + ++nb; + return nb; + } +} + //======================================================================= //function : NbAncestors //purpose : Return number of unique ancestors of the shape @@ -1586,7 +2171,7 @@ 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 )); } //================================================================================ @@ -1609,6 +2194,67 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape ) 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() ); +} + +//================================================================================ +/*! + * \brief Return type of shape contained in a group + * \param group - a shape of type TopAbs_COMPOUND + * \param avoidCompound - not to return TopAbs_COMPOUND + */ +//================================================================================ + +TopAbs_ShapeEnum SMESH_MesherHelper::GetGroupType(const TopoDS_Shape& group, + const bool avoidCompound) +{ + if ( !group.IsNull() ) + { + if ( group.ShapeType() != TopAbs_COMPOUND ) + return group.ShapeType(); + + // iterate on a compound + TopoDS_Iterator it( group ); + if ( it.More() ) + return avoidCompound ? GetGroupType( it.Value() ) : it.Value().ShapeType(); + } + return TopAbs_SHAPE; +} + //======================================================================= //function : IsQuadraticMesh //purpose : Check mesh without geometry for: if all elements on this shape are quadratic, @@ -1623,6 +2269,8 @@ SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh() int NbFacesAndEdges=0; //All faces and edges NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces(); + if ( NbAllEdgsAndFaces == 0 ) + return SMESH_MesherHelper::LINEAR; //Quadratic faces and edges NbQuadFacesAndEdgs = myMesh->NbEdges(ORDER_QUADRATIC) + myMesh->NbFaces(ORDER_QUADRATIC); @@ -1654,6 +2302,81 @@ 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 //======================================================================= @@ -1706,7 +2429,7 @@ 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 IsMoved() const { return (_nbMoves > 0 /*&& !IsStraight()*/); } bool IsStraight() const { return isStraightLink( (XYZ(node1())-XYZ(node2())).SquareModulus(), _nodeMove.SquareMagnitude()); @@ -1752,6 +2475,8 @@ namespace { // Structures used by FixQuadraticElements() const QLink* operator->() const { return _qlink; } gp_Vec Normal() const; + + bool IsStraight() const; }; // -------------------------------------------------------------------- typedef list< TChainLink > TChain; @@ -1929,9 +2654,10 @@ namespace { // Structures used by FixQuadraticElements() // 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(); @@ -1954,10 +2680,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 @@ -2151,7 +2878,7 @@ namespace { // Structures used by FixQuadraticElements() 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 @@ -2162,7 +2889,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 @@ -2241,7 +2968,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] @@ -2252,8 +2981,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] ); @@ -2303,6 +3045,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 @@ -2321,7 +3084,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 @@ -2370,7 +3133,7 @@ namespace { // Structures used by FixQuadraticElements() for ( linkIt = linkSet.begin(); linkIt != linksEnd; ++linkIt) { - if ( linkIt->IsBoundary() && !(*linkIt)->IsStraight() && linkIt->_qfaces[0]) + if ( linkIt->IsBoundary() && !linkIt->IsStraight() && linkIt->_qfaces[0]) { // move iff a boundary link is bent towards inside of a face (issue 0021084) const QFace* face = linkIt->_qfaces[0]; @@ -2568,20 +3331,355 @@ namespace { // Structures used by FixQuadraticElements() return _OK; } + + //================================================================================ + /*! + * \brief Place medium nodes at the link middle for elements whose corner nodes + * are out of geometrical boundary to prevent distorting elements. + * Issue 0020982, note 0013990 + */ + //================================================================================ + + void force3DOutOfBoundary( SMESH_MesherHelper& theHelper, + SMESH_ComputeErrorPtr& theError) + { + SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); + TopoDS_Shape shape = theHelper.GetSubShape().Oriented( TopAbs_FORWARD ); + if ( shape.IsNull() ) return; + + if ( !theError ) theError = SMESH_ComputeError::New(); + + gp_XYZ faceNorm; + + if ( shape.ShapeType() == TopAbs_FACE ) // 2D + { + if ( theHelper.GetMesh()->NbTriangles( ORDER_QUADRATIC ) < 1 ) return; + + SMESHDS_SubMesh* faceSM = meshDS->MeshElements( shape ); + if ( !faceSM ) return; + + const TopoDS_Face& face = TopoDS::Face( shape ); + Handle(Geom_Surface) surface = BRep_Tool::Surface( face ); + + TopExp_Explorer edgeIt( face, TopAbs_EDGE ); + for ( ; edgeIt.More(); edgeIt.Next() ) // loop on EDGEs of a FACE + { + // check if the EDGE needs checking + const TopoDS_Edge& edge = TopoDS::Edge( edgeIt.Current() ); + if ( BRep_Tool::Degenerated( edge ) ) + continue; + if ( theHelper.IsRealSeam( edge ) && + edge.Orientation() == TopAbs_REVERSED ) + continue; + + SMESHDS_SubMesh* edgeSM = meshDS->MeshElements( edge ); + if ( !edgeSM ) continue; + + double f,l; + Handle(Geom2d_Curve) pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l ); + BRepAdaptor_Curve curve3D( edge ); + switch ( curve3D.GetType() ) { + case GeomAbs_Line: continue; + case GeomAbs_Circle: + case GeomAbs_Ellipse: + case GeomAbs_Hyperbola: + case GeomAbs_Parabola: + try + { + gp_Vec D1, D2, Du1, Dv1; gp_Pnt p; + curve3D.D2( 0.5 * ( f + l ), p, D1, D2 ); + gp_Pnt2d uv = pcurve->Value( 0.5 * ( f + l ) ); + surface->D1( uv.X(), uv.Y(), p, Du1, Dv1 ); + gp_Vec fNorm = Du1 ^ Dv1; + if ( fNorm.IsParallel( D2, M_PI * 25./180. )) + continue; // face is normal to the curve3D + + gp_Vec curvNorm = fNorm ^ D1; + if ( edge.Orientation() == TopAbs_REVERSED ) curvNorm.Reverse(); + if ( curvNorm * D2 > 0 ) + continue; // convex edge + } + catch ( Standard_Failure ) + { + continue; + } + } + // get nodes shared by faces that may be distorted + SMDS_NodeIteratorPtr nodeIt; + if ( edgeSM->NbNodes() > 0 ) { + nodeIt = edgeSM->GetNodes(); + } + else { + SMESHDS_SubMesh* vertexSM = meshDS->MeshElements( theHelper.IthVertex( 0, edge )); + if ( !vertexSM ) + vertexSM = meshDS->MeshElements( theHelper.IthVertex( 1, edge )); + if ( !vertexSM ) continue; + nodeIt = vertexSM->GetNodes(); + } + + // find suspicious faces + TIDSortedElemSet checkedFaces; + vector< const SMDS_MeshNode* > nOnEdge( 2 ); + const SMDS_MeshNode* nOnFace; + while ( nodeIt->more() ) + { + const SMDS_MeshNode* n = nodeIt->next(); + SMDS_ElemIteratorPtr faceIt = n->GetInverseElementIterator( SMDSAbs_Face ); + while ( faceIt->more() ) + { + const SMDS_MeshElement* f = faceIt->next(); + if ( !faceSM->Contains( f ) || + f->NbNodes() != 6 || // check quadratic triangles only + !checkedFaces.insert( f ).second ) + continue; + + // get nodes on EDGE and on FACE of a suspicious face + nOnEdge.clear(); nOnFace = 0; + SMDS_MeshElement::iterator triNode = f->begin_nodes(); + for ( int nbN = 0; nbN < 3; ++triNode, ++nbN ) + { + n = *triNode; + if ( n->GetPosition()->GetDim() == 2 ) + nOnFace = n; + else + nOnEdge.push_back( n ); + } + + // check if nOnFace is inside the FACE + if ( nOnFace && nOnEdge.size() == 2 ) + { + theHelper.AddTLinks( static_cast< const SMDS_MeshFace* > ( f )); + if ( !SMESH_Algo::FaceNormal( f, faceNorm, /*normalized=*/false )) + continue; + gp_XYZ edgeDir = SMESH_TNodeXYZ( nOnEdge[0] ) - SMESH_TNodeXYZ( nOnEdge[1] ); + gp_XYZ edgeNorm = faceNorm ^ edgeDir; + n = theHelper.GetMediumNode( nOnEdge[0], nOnEdge[1], true ); + gp_XYZ pN0 = SMESH_TNodeXYZ( nOnEdge[0] ); + gp_XYZ pMedium = SMESH_TNodeXYZ( n ); // on-edge node location + gp_XYZ pFaceN = SMESH_TNodeXYZ( nOnFace ); // on-face node location + double hMedium = edgeNorm * gp_Vec( pN0, pMedium ).XYZ(); + double hFace = edgeNorm * gp_Vec( pN0, pFaceN ).XYZ(); + if ( Abs( hMedium ) > Abs( hFace * 0.6 )) + { + // nOnFace is out of FACE, move a medium on-edge node to the middle + gp_XYZ pMid3D = 0.5 * ( pN0 + SMESH_TNodeXYZ( nOnEdge[1] )); + meshDS->MoveNode( n, pMid3D.X(), pMid3D.Y(), pMid3D.Z() ); + MSG( "move OUT of face " << n ); + theError->myBadElements.push_back( f ); + } + } + } + } + } + if ( !theError->myBadElements.empty() ) + theError->myName = EDITERR_NO_MEDIUM_ON_GEOM; + return; + + } // 2D ============================================================================== + + if ( shape.ShapeType() == TopAbs_SOLID ) // 3D + { + if ( theHelper.GetMesh()->NbTetras ( ORDER_QUADRATIC ) < 1 && + theHelper.GetMesh()->NbPyramids( ORDER_QUADRATIC ) < 1 ) return; + + SMESHDS_SubMesh* solidSM = meshDS->MeshElements( shape ); + if ( !solidSM ) return; + + // check if the SOLID is bound by concave FACEs + vector< TopoDS_Face > concaveFaces; + TopExp_Explorer faceIt( shape, TopAbs_FACE ); + for ( ; faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID + { + const TopoDS_Face& face = TopoDS::Face( faceIt.Current() ); + if ( !meshDS->MeshElements( face )) continue; + + BRepAdaptor_Surface surface( face ); + switch ( surface.GetType() ) { + case GeomAbs_Plane: continue; + case GeomAbs_Cylinder: + case GeomAbs_Cone: + case GeomAbs_Sphere: + try + { + double u = 0.5 * ( surface.FirstUParameter() + surface.LastUParameter() ); + double v = 0.5 * ( surface.FirstVParameter() + surface.LastVParameter() ); + gp_Vec Du1, Dv1, Du2, Dv2, Duv2; gp_Pnt p; + surface.D2( u,v, p, Du1, Dv1, Du2, Dv2, Duv2 ); + gp_Vec fNorm = Du1 ^ Dv1; + if ( face.Orientation() == TopAbs_REVERSED ) fNorm.Reverse(); + bool concaveU = ( fNorm * Du2 > 1e-100 ); + bool concaveV = ( fNorm * Dv2 > 1e-100 ); + if ( concaveU || concaveV ) + concaveFaces.push_back( face ); + } + catch ( Standard_Failure ) + { + concaveFaces.push_back( face ); + } + } + } + if ( concaveFaces.empty() ) + return; + + // fix 2D mesh on the SOLID + for ( faceIt.ReInit(); faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID + { + SMESH_MesherHelper faceHelper( *theHelper.GetMesh() ); + faceHelper.SetSubShape( faceIt.Current() ); + force3DOutOfBoundary( faceHelper, theError ); + } + + // get an iterator over faces on concaveFaces + vector< SMDS_ElemIteratorPtr > faceIterVec( concaveFaces.size() ); + for ( size_t i = 0; i < concaveFaces.size(); ++i ) + faceIterVec[i] = meshDS->MeshElements( concaveFaces[i] )->GetElements(); + typedef SMDS_IteratorOnIterators + < const SMDS_MeshElement*, vector< SMDS_ElemIteratorPtr > > TIterOnIter; + SMDS_ElemIteratorPtr faceIter( new TIterOnIter( faceIterVec )); + + // a seacher to check if a volume is close to a concave face + std::auto_ptr< SMESH_ElementSearcher > faceSearcher + ( SMESH_MeshEditor( theHelper.GetMesh() ).GetElementSearcher( faceIter )); + + // classifier + //BRepClass3d_SolidClassifier solidClassifier( shape ); + + TIDSortedElemSet checkedVols, movedNodes; + for ( faceIt.ReInit(); faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID + { + const TopoDS_Shape& face = faceIt.Current(); + SMESHDS_SubMesh* faceSM = meshDS->MeshElements( face ); + if ( !faceSM ) continue; + + // get nodes shared by volumes (tet and pyra) on the FACE that may be distorted + SMDS_NodeIteratorPtr nodeIt; + if ( faceSM->NbNodes() > 0 ) { + nodeIt = faceSM->GetNodes(); + } + else { + TopExp_Explorer vertex( face, TopAbs_VERTEX ); + SMESHDS_SubMesh* vertexSM = meshDS->MeshElements( vertex.Current() ); + if ( !vertexSM ) continue; + nodeIt = vertexSM->GetNodes(); + } + + // find suspicious volumes adjacent to the FACE + vector< const SMDS_MeshNode* > nOnFace( 4 ); + const SMDS_MeshNode* nInSolid; + //vector< const SMDS_MeshElement* > intersectedFaces; + while ( nodeIt->more() ) + { + const SMDS_MeshNode* n = nodeIt->next(); + SMDS_ElemIteratorPtr volIt = n->GetInverseElementIterator( SMDSAbs_Volume ); + while ( volIt->more() ) + { + const SMDS_MeshElement* vol = volIt->next(); + int nbN = vol->NbCornerNodes(); + if ( ( nbN != 4 && nbN != 5 ) || + !solidSM->Contains( vol ) || + !checkedVols.insert( vol ).second ) + continue; + + // get nodes on FACE and in SOLID of a suspicious volume + nOnFace.clear(); nInSolid = 0; + SMDS_MeshElement::iterator volNode = vol->begin_nodes(); + for ( int nb = nbN; nb > 0; ++volNode, --nb ) + { + n = *volNode; + if ( n->GetPosition()->GetDim() == 3 ) + nInSolid = n; + else + nOnFace.push_back( n ); + } + if ( !nInSolid || nOnFace.size() != nbN - 1 ) + continue; + + // get size of the vol + SMESH_TNodeXYZ pInSolid( nInSolid ), pOnFace0( nOnFace[0] ); + double volLength = pInSolid.SquareDistance( nOnFace[0] ); + for ( size_t i = 1; i < nOnFace.size(); ++i ) + { + volLength = Max( volLength, pOnFace0.SquareDistance( nOnFace[i] )); + } + + // check if vol is close to concaveFaces + const SMDS_MeshElement* closeFace = + faceSearcher->FindClosestTo( pInSolid, SMDSAbs_Face ); + if ( !closeFace || + pInSolid.SquareDistance( closeFace->GetNode(0) ) > 4 * volLength ) + continue; + + // check if vol is distorted, i.e. a medium node is much closer + // to nInSolid than the link middle + bool isDistorted = false; + SMDS_FaceOfNodes onFaceTria( nOnFace[0], nOnFace[1], nOnFace[2] ); + if ( !SMESH_Algo::FaceNormal( &onFaceTria, faceNorm, /*normalized=*/false )) + continue; + theHelper.AddTLinks( static_cast< const SMDS_MeshVolume* > ( vol )); + vector< pair< SMESH_TLink, const SMDS_MeshNode* > > links; + for ( size_t i = 0; i < nOnFace.size(); ++i ) // loop on links between nOnFace + for ( size_t j = i+1; j < nOnFace.size(); ++j ) + { + SMESH_TLink link( nOnFace[i], nOnFace[j] ); + TLinkNodeMap::const_iterator linkIt = + theHelper.GetTLinkNodeMap().find( link ); + if ( linkIt != theHelper.GetTLinkNodeMap().end() ) + { + links.push_back( make_pair( linkIt->first, linkIt->second )); + if ( !isDistorted ) { + // compare projections of nInSolid and nMedium to face normal + gp_Pnt pMedium = SMESH_TNodeXYZ( linkIt->second ); + double hMedium = faceNorm * gp_Vec( pOnFace0, pMedium ).XYZ(); + double hVol = faceNorm * gp_Vec( pOnFace0, pInSolid ).XYZ(); + isDistorted = ( Abs( hMedium ) > Abs( hVol * 0.5 )); + } + } + } + // move medium nodes to link middle + if ( isDistorted ) + { + for ( size_t i = 0; i < links.size(); ++i ) + { + const SMDS_MeshNode* nMedium = links[i].second; + if ( movedNodes.insert( nMedium ).second ) + { + gp_Pnt pMid3D = 0.5 * ( SMESH_TNodeXYZ( links[i].first.node1() ) + + SMESH_TNodeXYZ( links[i].first.node2() )); + meshDS->MoveNode( nMedium, pMid3D.X(), pMid3D.Y(), pMid3D.Z() ); + MSG( "move OUT of solid " << nMedium ); + } + } + theError->myBadElements.push_back( vol ); + } + } // loop on volumes sharing a node on FACE + } // loop on nodes on FACE + } // loop on FACEs of a SOLID + + if ( !theError->myBadElements.empty() ) + theError->myName = EDITERR_NO_MEDIUM_ON_GEOM; + } // 3D case + } + } //namespace //======================================================================= /*! * \brief Move medium nodes of faces and volumes to fix distorted elements + * \param error - container of fixed distorted elements * \param volumeOnly - to fix nodes on faces or not, if the shape is solid * * Issue 0020307: EDF 992 SMESH : Linea/Quadratic with Medium Node on Geometry */ //======================================================================= -void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) +void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, + bool volumeOnly) { - // 0. 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 FACEs // ---------------------------------------------- if ( myShape.IsNull() ) { if ( !myMesh->HasShapeToMesh() ) return; @@ -2603,23 +3701,30 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) 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); + h.ToFixNodeParameters(true); + h.FixQuadraticElements( compError, false ); } } // fix nodes on geom faces #ifdef _DEBUG_ - int nbfaces = faces.Extent(); + 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); + h.ToFixNodeParameters(true); + h.FixQuadraticElements( compError, true); } //perf_print_all_meters(1); + if ( compError && compError->myName == EDITERR_NO_MEDIUM_ON_GEOM ) + compError->myComment = "during conversion to quadratic, " + "some medium nodes were not placed on geometry to avoid distorting elements"; return; } @@ -2653,15 +3758,24 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) bool isCurved = false; //bool hasRectFaces = false; //set nbElemNodeSet; + SMDS_VolumeTool volTool; + + TIDSortedNodeSet apexOfPyramid; + const int apexIndex = 4; + + // Issue 0020982 + // Move medium nodes to the link middle for elements whose corner nodes + // are out of geometrical boundary to fix distorted elements. + force3DOutOfBoundary( *this, compError ); 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 ); @@ -2674,10 +3788,17 @@ 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; @@ -2695,6 +3816,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 ) @@ -2716,7 +3840,9 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) QLink link( face->GetNode(iN), face->GetNode((iN+1)%nbN), face->GetNode(iN+nbN) ); pLink = links.insert( link ).first; faceLinks[ iN ] = & *pLink; - if ( !isCurved ) + if ( !isCurved && + link.node1()->GetPosition()->GetTypeOfPosition() < 2 && + link.node2()->GetPosition()->GetTypeOfPosition() < 2 ) isCurved = !link.IsStraight(); } // store QFace @@ -2729,9 +3855,9 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) return; // no curved edges of faces // 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; @@ -2743,7 +3869,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) 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 @@ -2776,12 +3902,12 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) { TChain& chain = chains[iC]; if ( chain.empty() ) continue; - if ( chain.front()->IsStraight() && chain.back()->IsStraight() ) { + if ( chain.front().IsStraight() && chain.back().IsStraight() ) { MSG("3D straight - ignore"); continue; } if ( chain.front()->MediumPos() > bndPos || - chain.back()->MediumPos() > bndPos ) { + chain.back() ->MediumPos() > bndPos ) { MSG("Internal chain - ignore"); continue; } @@ -2811,9 +3937,11 @@ 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 ) { @@ -2832,14 +3960,24 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) ( 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(),uvMove.SquareModulus()); + isStraight[is1] = isStraightLink( (uv2-uv1).SquareModulus(), + 10 * uvMove.SquareModulus()); } -// if ( move0.SquareMagnitude() < straightTol2 && -// move1.SquareMagnitude() < straightTol2 ) { 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; @@ -2906,59 +4044,73 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) } // loop on chains of links } // loop on 2 directions of propagation from quadrangle } // loop on faces - } + } // fix faces and/or volumes // 4. Move nodes - // ----------- + // ------------- 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()); } } -} - -//======================================================================= -/*! - * \brief Iterator on ancestors of the given type - */ -//======================================================================= - -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() - { - 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; - } -}; -//======================================================================= -/*! - * \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)); + // Issue 0020982 + // Move the apex of pyramid together with the most curved link. + // TIDSortedNodeSet::iterator apexIt = apexOfPyramid.begin(); + // for ( ; apexIt != apexOfPyramid.end(); ++apexIt ) + // { + // SMESH_TNodeXYZ apex = *apexIt; + + // gp_Vec maxMove( 0,0,0 ); + // double maxMoveSize2 = 0; + + // // 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()); + // } + // } + // } } +