X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=eceb92e3657040b6f37d0ec777c370f534dce5bf;hp=eb3ba5d7194236d9fbb5beb6c83573ccf78c0e6e;hb=1fb3ff621ac19da15742a5e8b9253b594977d14a;hpb=a713f0b919f062a1c9839670d4b51639c8133940 diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index eb3ba5d71..eceb92e36 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2019 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 @@ -31,8 +31,12 @@ #include "SMDS_FacePosition.hxx" #include "SMDS_IteratorOnIterators.hxx" #include "SMDS_VolumeTool.hxx" +#include "SMESHDS_Mesh.hxx" #include "SMESH_Block.hxx" +#include "SMESH_HypoFilter.hxx" +#include "SMESH_Mesh.hxx" #include "SMESH_MeshAlgos.hxx" +#include "SMESH_MeshEditor.hxx" #include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" @@ -70,7 +74,7 @@ using namespace std; namespace { - gp_XYZ XYZ(const SMDS_MeshNode* n) { return gp_XYZ(n->X(), n->Y(), n->Z()); } + inline SMESH_NodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_NodeXYZ(n); } enum { U_periodic = 1, V_periodic = 2 }; } @@ -112,11 +116,33 @@ SMESH_MesherHelper::~SMESH_MesherHelper() } } +//================================================================================ +/*! + * \brief Return SMESH_Gen + */ +//================================================================================ + +SMESH_Gen* SMESH_MesherHelper::GetGen() const +{ + return GetMesh()->GetGen(); +} + +//================================================================================ +/*! + * \brief Return mesh DS + */ +//================================================================================ + +SMESHDS_Mesh* SMESH_MesherHelper::GetMeshDS() const +{ + return GetMesh()->GetMeshDS(); +} + //======================================================================= //function : IsQuadraticSubMesh -//purpose : Check submesh for given shape: if all elements on this shape +//purpose : Check sub-meshes of a given shape: if all elements on sub-shapes // are quadratic, quadratic elements will be created. -// Also fill myTLinkNodeMap +// Fill myTLinkNodeMap //======================================================================= bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) @@ -124,7 +150,6 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) SMESHDS_Mesh* meshDS = GetMeshDS(); // we can create quadratic elements only if all elements // created on sub-shapes of given shape are quadratic - // also we have to fill myTLinkNodeMap myCreateQuadratic = true; mySeamShapeIds.clear(); myDegenShapeIds.clear(); @@ -137,9 +162,6 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) } SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge ); - - int nbOldLinks = myTLinkNodeMap.size(); - if ( !myMesh->HasShapeToMesh() ) { if (( myCreateQuadratic = myMesh->NbFaces( ORDER_QUADRATIC ))) @@ -190,12 +212,13 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) } } - if ( nbOldLinks == myTLinkNodeMap.size() ) + // if ( nbOldLinks == myTLinkNodeMap.size() ) -- 0023068 + if ( myTLinkNodeMap.empty() ) myCreateQuadratic = false; - if(!myCreateQuadratic) { + if ( !myCreateQuadratic ) myTLinkNodeMap.clear(); - } + SetSubShape( aSh ); return myCreateQuadratic; @@ -278,6 +301,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) double u2 = uv1.Coord(1); myPar1[0] = Min( u1, u2 ); myPar2[0] = Max( u1, u2 ); + myParIndex |= U_periodic; } else { @@ -287,6 +311,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) double v2 = uv1.Coord(2); myPar1[1] = Min( v1, v2 ); myPar2[1] = Max( v1, v2 ); + myParIndex |= V_periodic; } } else //if ( !isSeam ) @@ -306,11 +331,18 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) { double f,l, r = 0.2345; Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( edge, face, f, l ); - uv2 = C2d->Value( f * r + l * ( 1.-r )); - if ( du < Precision::PConfusion() ) - isSeam = ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Precision::PConfusion() ); + if ( C2d.IsNull() ) + { + isSeam = false; + } else - isSeam = ( Abs( uv1.Coord(2) - uv2.Coord(2) ) < Precision::PConfusion() ); + { + uv2 = C2d->Value( f * r + l * ( 1.-r )); + if ( du < Precision::PConfusion() ) + isSeam = ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Precision::PConfusion() ); + else + isSeam = ( Abs( uv1.Coord(2) - uv2.Coord(2) ) < Precision::PConfusion() ); + } } } if ( isSeam ) @@ -338,6 +370,37 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) } } +//======================================================================= +/*! + * \brief Copy shape information from another helper. Used to improve performance + * since SetSubShape() can be time consuming if there are many edges + */ +//======================================================================= + +void SMESH_MesherHelper::CopySubShapeInfo(const SMESH_MesherHelper& other) +{ + this->myShape = other.myShape; + this->myShapeID = other.myShapeID; + this->myDegenShapeIds = other.myDegenShapeIds; + this->mySeamShapeIds = other.mySeamShapeIds; + this->myPar1[0] = other.myPar1[0]; + this->myPar1[1] = other.myPar1[1]; + this->myPar2[0] = other.myPar2[0]; + this->myPar2[1] = other.myPar2[1]; + this->myParIndex = other.myParIndex; + this->myFace2Surface = other.myFace2Surface; +} + +//======================================================================= +//function : ShapeToIndex +//purpose : Convert a shape to its index in the SMESHDS_Mesh +//======================================================================= + +int SMESH_MesherHelper::ShapeToIndex( const TopoDS_Shape& S ) const +{ + return GetMeshDS()->ShapeToIndex( S ); +} + //======================================================================= //function : GetNodeUVneedInFaceNode //purpose : Check if inFaceNode argument is necessary for call GetNodeUV(F,..) @@ -581,22 +644,22 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, { gp_Pnt2d uv( Precision::Infinite(), Precision::Infinite() ); - const SMDS_PositionPtr Pos = n->GetPosition(); + SMDS_PositionPtr pos = n->GetPosition(); bool uvOK = false; - if ( Pos->GetTypeOfPosition() == SMDS_TOP_FACE ) + if ( pos->GetTypeOfPosition() == SMDS_TOP_FACE ) { // node has position on face - const SMDS_FacePosition* fpos = static_cast( Pos ); + SMDS_FacePositionPtr fpos = pos; uv.SetCoord( fpos->GetUParameter(), fpos->GetVParameter() ); if ( check ) uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 2.*getFaceMaxTol( F )); // 2. from 22830 } - else if ( Pos->GetTypeOfPosition() == SMDS_TOP_EDGE ) + else if ( pos->GetTypeOfPosition() == SMDS_TOP_EDGE ) { // node has position on EDGE => it is needed to find // corresponding EDGE from FACE, get pcurve for this // EDGE and retrieve value from this pcurve - const SMDS_EdgePosition* epos = static_cast( Pos ); + SMDS_EdgePositionPtr epos = pos; const int edgeID = n->getshapeId(); const TopoDS_Edge& E = TopoDS::Edge( GetMeshDS()->IndexToShape( edgeID )); double f, l, u = epos->GetUParameter(); @@ -637,7 +700,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, uv = newUV; } } - else if ( Pos->GetTypeOfPosition() == SMDS_TOP_VERTEX ) + else if ( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX ) { if ( int vertexID = n->getshapeId() ) { const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID)); @@ -647,12 +710,11 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, } catch (Standard_Failure& exc) { } - if ( !uvOK ) { - for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !uvOK && vert.More(); vert.Next() ) - uvOK = ( V == vert.Current() ); - if ( !uvOK ) { - MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID - << " not in face " << GetMeshDS()->ShapeToIndex( F ) ); + if ( !uvOK ) + { + if ( !IsSubShape( V, F )) + { + MESSAGE("GetNodeUV() Vertex "<< vertexID <<" not in face "<< GetMeshDS()->ShapeToIndex(F)); // get UV of a vertex closest to the node double dist = 1e100; gp_Pnt pn = XYZ( n ); @@ -667,7 +729,8 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, } } } - else { + else + { uvOK = false; TopTools_ListIteratorOfListOfShape it( myMesh->GetAncestors( V )); for ( ; it.More(); it.Next() ) { @@ -676,13 +739,23 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, double f,l; Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(edge, F, f, l); if ( !C2d.IsNull() ) { - double u = ( V == TopExp::FirstVertex( edge ) ) ? f : l; + double u = ( V == IthVertex( 0, edge )) ? f : l; uv = C2d->Value( u ); uvOK = true; break; } } } + if ( !uvOK && V.Orientation() == TopAbs_INTERNAL ) + { + Handle(ShapeAnalysis_Surface) projector = GetSurface( F ); + if ( n2 ) uv = GetNodeUV( F, n2 ); + if ( Precision::IsInfinite( uv.X() )) + uv = projector->NextValueOfUV( uv, BRep_Tool::Pnt( V ), BRep_Tool::Tolerance( F )); + else + uv = projector->ValueOfUV( BRep_Tool::Pnt( V ), BRep_Tool::Tolerance( F )); + uvOK = ( projector->Gap() < getFaceMaxTol( F )); + } } } if ( n2 && IsSeamShape( vertexID )) @@ -732,33 +805,20 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, // check that uv is correct TopLoc_Location loc; Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc ); - gp_Pnt nodePnt = XYZ( n ), surfPnt(0,0,0); + SMESH_NodeXYZ nXYZ( n ); + gp_Pnt nodePnt = nXYZ, surfPnt(0,0,0); double dist = 0; if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() ); 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 ); - if ( !projector.IsDone() || projector.NbPoints() < 1 ) - { - MESSAGE( "SMESH_MesherHelper::CheckNodeUV() failed to project" ); - return false; - } - Quantity_Parameter U,V; - projector.LowerDistanceParameters(U,V); - uv.SetCoord( U,V ); - surfPnt = surface->Value( U, V ); - dist = nodePnt.Distance( surfPnt ); + Handle(ShapeAnalysis_Surface) sprojector = GetSurface( F ); + uv = sprojector->ValueOfUV( nXYZ, tol ).XY(); + surfPnt = sprojector->Value( uv ); + dist = surfPnt.Distance( nXYZ ); if ( distXYZ ) { - surfPnt.Transform( loc ); distXYZ[0] = dist; distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z(); } @@ -770,7 +830,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, // store the fixed UV on the face if ( myShape.IsSame(F) && shapeID == myShapeID && myFixNodeParameters ) const_cast(n)->SetPosition - ( SMDS_PositionPtr( new SMDS_FacePosition( U, V ))); + ( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() ))); } else if ( myShape.IsSame(F) && uv.Modulus() > numeric_limits::min() ) { @@ -782,7 +842,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, //======================================================================= //function : GetProjector -//purpose : Return projector intitialized by given face without location, which is returned +//purpose : Return projector initialized by given face without location, which is returned //======================================================================= GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face& F, @@ -805,6 +865,24 @@ GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face& return *( i_proj->second ); } +//======================================================================= +//function : GetSurface +//purpose : Return a cached ShapeAnalysis_Surface of a FACE +//======================================================================= + +Handle(ShapeAnalysis_Surface) SMESH_MesherHelper::GetSurface(const TopoDS_Face& F ) const +{ + Handle(Geom_Surface) surface = BRep_Tool::Surface( F ); + int faceID = GetMeshDS()->ShapeToIndex( F ); + TID2Surface::iterator i_surf = myFace2Surface.find( faceID ); + if ( i_surf == myFace2Surface.end() && faceID ) + { + Handle(ShapeAnalysis_Surface) surf( new ShapeAnalysis_Surface( surface )); + i_surf = myFace2Surface.insert( make_pair( faceID, surf )).first; + } + return i_surf->second; +} + namespace { gp_XY AverageUV(const gp_XY& uv1, const gp_XY& uv2) { return ( uv1 + uv2 ) / 2.; } @@ -813,18 +891,20 @@ namespace } //======================================================================= -//function : applyIn2D +//function : ApplyIn2D //purpose : Perform given operation on two 2d points in parameric space of given surface. // It takes into account period of the surface. Use gp_XY_FunPtr macro // to easily define pointer to function of gp_XY class. //======================================================================= -gp_XY SMESH_MesherHelper::applyIn2D(const Handle(Geom_Surface)& surface, - const gp_XY& uv1, - const gp_XY& uv2, - xyFunPtr fun, - const bool resultInPeriod) +gp_XY SMESH_MesherHelper::ApplyIn2D(Handle(Geom_Surface) surface, + const gp_XY& uv1, + const gp_XY& uv2, + xyFunPtr fun, + const bool resultInPeriod) { + if ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface ))) + surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface(); Standard_Boolean isUPeriodic = surface.IsNull() ? false : surface->IsUPeriodic(); Standard_Boolean isVPeriodic = surface.IsNull() ? false : surface->IsVPeriodic(); if ( !isUPeriodic && !isVPeriodic ) @@ -852,6 +932,31 @@ gp_XY SMESH_MesherHelper::applyIn2D(const Handle(Geom_Surface)& surface, return res; } + +//======================================================================= +//function : AdjustByPeriod +//purpose : Move node positions on a FACE within surface period +//======================================================================= + +void SMESH_MesherHelper::AdjustByPeriod( const TopoDS_Face& face, gp_XY uv[], const int nbUV ) +{ + SMESH_MesherHelper h( *myMesh ), *ph = face.IsSame( myShape ) ? this : &h; + ph->SetSubShape( face ); + + for ( int iCoo = U_periodic; iCoo <= V_periodic; ++iCoo ) + if ( ph->GetPeriodicIndex() & iCoo ) + { + const double period = ( ph->myPar2[iCoo-1] - ph->myPar1[iCoo-1] ); + const double xRef = uv[0].Coord( iCoo ); + for ( int i = 1; i < nbUV; ++i ) + { + double x = uv[i].Coord( iCoo ); + double dx = ShapeAnalysis::AdjustByPeriod( x, xRef, period ); + uv[i].SetCoord( iCoo, x + dx ); + } + } +} + //======================================================================= //function : GetMiddleUV //purpose : Return middle UV taking in account surface period @@ -862,13 +967,13 @@ gp_XY SMESH_MesherHelper::GetMiddleUV(const Handle(Geom_Surface)& surface, const gp_XY& p2) { // NOTE: - // the proper place of getting basic surface seems to be in applyIn2D() + // 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(); + // 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 ); + return ApplyIn2D( surface, p1, p2, & AverageUV ); } //======================================================================= @@ -914,8 +1019,7 @@ double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge& E, const SMDS_PositionPtr pos = n->GetPosition(); if ( pos->GetTypeOfPosition()==SMDS_TOP_EDGE ) { - const SMDS_EdgePosition* epos = static_cast( pos ); - param = epos->GetUParameter(); + param = pos->GetParameters()[0]; } else if( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX ) { @@ -932,6 +1036,16 @@ double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge& E, int vertexID = n->getshapeId(); const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID)); param = BRep_Tool::Parameter( V, E ); + + if ( inEdgeNode ) + { + BRepAdaptor_Curve curve( E ); + if ( curve.IsPeriodic() ) + { + double uInEdge = GetNodeU( E, inEdgeNode ); + param += ShapeAnalysis::AdjustByPeriod( param, uInEdge, curve.Period() ); + } + } } } if ( check ) @@ -1013,9 +1127,8 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, MESSAGE( "SMESH_MesherHelper::CheckNodeU() failed to project" ); return false; } - Quantity_Parameter U = projector->LowerDistanceParameter(); + Standard_Real U = projector->LowerDistanceParameter(); u = double( U ); - MESSAGE(" f " << f << " l " << l << " u " << u); curvPnt = curve->Value( u ); dist = nodePnt.Distance( curvPnt ); if ( distXYZ ) { @@ -1025,8 +1138,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, } if ( dist > tol ) { - MESSAGE( "SMESH_MesherHelper::CheckNodeU(), invalid projection" ); - MESSAGE("distance " << dist << " " << tol ); + MESSAGE( "CheckNodeU(), invalid projection; distance " << dist << "; tol " << tol ); return false; } // store the fixed U on the edge @@ -1266,14 +1378,24 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1, if ( solidID < 1 && !faceId2nbNodes.empty() ) // SOLID not found { // find ID of the FACE the four corner nodes belong to - itMapWithIdFace = faceId2nbNodes.begin(); - for ( ; itMapWithIdFace != faceId2nbNodes.end(); ++itMapWithIdFace) + itMapWithIdFace = faceId2nbNodes.find( myShapeID ); // IPAL52698 + if ( itMapWithIdFace != faceId2nbNodes.end() && + itMapWithIdFace->second == 4 ) + { + shapeType = TopAbs_FACE; + faceID = myShapeID; + } + else { - if ( itMapWithIdFace->second == 4 ) + itMapWithIdFace = faceId2nbNodes.begin(); + for ( ; itMapWithIdFace != faceId2nbNodes.end(); ++itMapWithIdFace) { - shapeType = TopAbs_FACE; - faceID = (*itMapWithIdFace).first; - break; + if ( itMapWithIdFace->second == 4 ) + { + shapeType = TopAbs_FACE; + faceID = (*itMapWithIdFace).first; + break; + } } } } @@ -1291,14 +1413,34 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1, bool toCheck = true; if ( !F.IsNull() && !force3d ) { - uvAvg = calcTFI (0.5, 0.5, - GetNodeUV(F,n1,n3,&toCheck), GetNodeUV(F,n2,n4,&toCheck), - GetNodeUV(F,n3,n1,&toCheck), GetNodeUV(F,n4,n2,&toCheck), - GetNodeUV(F,n12,n3), GetNodeUV(F,n23,n4), - GetNodeUV(F,n34,n2), GetNodeUV(F,n41,n2)); - TopLoc_Location loc; - Handle( Geom_Surface ) S = BRep_Tool::Surface( F, loc ); - P = S->Value( uvAvg.X(), uvAvg.Y() ).Transformed( loc ); + Handle(ShapeAnalysis_Surface) surface = GetSurface( F ); + if ( HasDegeneratedEdges() || surface->HasSingularities( 1e-7 )) + { + gp_Pnt center = calcTFI (0.5, 0.5, // IPAL0052863 + SMESH_TNodeXYZ(n1), SMESH_TNodeXYZ(n2), + SMESH_TNodeXYZ(n3), SMESH_TNodeXYZ(n4), + SMESH_TNodeXYZ(n12), SMESH_TNodeXYZ(n23), + SMESH_TNodeXYZ(n34), SMESH_TNodeXYZ(n41)); + gp_Pnt2d uv12 = GetNodeUV( F, n12, n3, &toCheck ); + uvAvg = surface->NextValueOfUV( uv12, center, BRep_Tool::Tolerance( F )).XY(); + } + else + { + gp_XY uv[8] = { + GetNodeUV( F,n1, n3, &toCheck ), + GetNodeUV( F,n2, n4, &toCheck ), + GetNodeUV( F,n3, n1, &toCheck ), + GetNodeUV( F,n4, n2, &toCheck ), + GetNodeUV( F,n12, n3 ), + GetNodeUV( F,n23, n4 ), + GetNodeUV( F,n34, n2 ), + GetNodeUV( F,n41, n2 ) + }; + AdjustByPeriod( F, uv, 8 ); // put uv[] within a period (IPAL52698) + + uvAvg = calcTFI (0.5, 0.5, uv[0],uv[1],uv[2],uv[3], uv[4],uv[5],uv[6],uv[7] ); + } + P = surface->Value( uvAvg ); centralNode = meshDS->AddNode( P.X(), P.Y(), P.Z() ); // if ( mySetElemOnShape ) node is not elem! meshDS->SetNodeOnFace( centralNode, faceID, uvAvg.X(), uvAvg.Y() ); @@ -1307,7 +1449,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1, { P = calcTFI (0.5, 0.5, SMESH_TNodeXYZ(n1), SMESH_TNodeXYZ(n2), - SMESH_TNodeXYZ(n3), SMESH_TNodeXYZ(n4), + SMESH_TNodeXYZ(n3), SMESH_TNodeXYZ(n4), SMESH_TNodeXYZ(n12), SMESH_TNodeXYZ(n23), SMESH_TNodeXYZ(n34), SMESH_TNodeXYZ(n41)); centralNode = meshDS->AddNode( P.X(), P.Y(), P.Z() ); @@ -1411,14 +1553,24 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1, if ( solidID < 1 && !faceId2nbNodes.empty() ) // SOLID not found { // find ID of the FACE the four corner nodes belong to - itMapWithIdFace = faceId2nbNodes.begin(); - for ( ; itMapWithIdFace != faceId2nbNodes.end(); ++itMapWithIdFace) + itMapWithIdFace = faceId2nbNodes.find( myShapeID ); // IPAL52698 + if ( itMapWithIdFace != faceId2nbNodes.end() && + itMapWithIdFace->second == 4 ) { - if ( itMapWithIdFace->second == 3 ) + shapeType = TopAbs_FACE; + faceID = myShapeID; + } + else + { + itMapWithIdFace = faceId2nbNodes.begin(); + for ( ; itMapWithIdFace != faceId2nbNodes.end(); ++itMapWithIdFace) { - shapeType = TopAbs_FACE; - faceID = (*itMapWithIdFace).first; - break; + if ( itMapWithIdFace->second == 3 ) + { + shapeType = TopAbs_FACE; + faceID = (*itMapWithIdFace).first; + break; + } } } } @@ -1430,13 +1582,18 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1, { F = TopoDS::Face( meshDS->IndexToShape( faceID )); bool checkOK = true, badTria = false; - gp_XY uv1 = GetNodeUV( F, n1, n23, &checkOK ); - gp_XY uv2 = GetNodeUV( F, n2, n31, &checkOK ); - gp_XY uv3 = GetNodeUV( F, n3, n12, &checkOK ); - gp_XY uv12 = GetNodeUV( F, n12, n3, &checkOK ); - gp_XY uv23 = GetNodeUV( F, n23, n1, &checkOK ); - gp_XY uv31 = GetNodeUV( F, n31, n2, &checkOK ); - uvAvg = GetCenterUV( uv1,uv2,uv3, uv12,uv23,uv31, &badTria ); + gp_XY uv[6] = { + GetNodeUV( F, n1, n23, &checkOK ), + GetNodeUV( F, n2, n31, &checkOK ), + GetNodeUV( F, n3, n12, &checkOK ), + GetNodeUV( F, n12, n3, &checkOK ), + GetNodeUV( F, n23, n1, &checkOK ), + GetNodeUV( F, n31, n2, &checkOK ) + }; + AdjustByPeriod( F, uv, 6 ); // put uv[] within a period (IPAL52698) + + uvAvg = GetCenterUV( uv[0],uv[1],uv[2], uv[3],uv[4],uv[5], &badTria ); + if ( badTria || !checkOK ) force3d = true; } @@ -1506,7 +1663,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, // get type of shape for the new medium node int faceID = -1, edgeID = -1; - TopoDS_Edge E; double u [2]; + TopoDS_Edge E; double u [2] = {0.,0.}; TopoDS_Face F; gp_XY uv[2]; bool uvOK[2] = { true, true }; const bool useCurSubShape = ( !myShape.IsNull() && myShape.ShapeType() == TopAbs_EDGE ); @@ -1518,6 +1675,28 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, { F = TopoDS::Face(meshDS->IndexToShape( faceID = pos.first )); uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]); + if (( !force3d ) && + ( HasDegeneratedEdges() || GetSurface( F )->HasSingularities( 1e-7 ))) + { + // IPAL52850 (degen VERTEX not at singularity) + // project middle point to a surface + SMESH_TNodeXYZ p1( n1 ), p2( n2 ); + gp_Pnt pMid = 0.5 * ( p1 + p2 ); + Handle(ShapeAnalysis_Surface) projector = GetSurface( F ); + gp_Pnt2d uvMid; + if ( uvOK[0] ) + uvMid = projector->NextValueOfUV( uv[0], pMid, BRep_Tool::Tolerance( F )); + else + uvMid = projector->ValueOfUV( pMid, getFaceMaxTol( F )); + if ( projector->Gap() * projector->Gap() < ( p1 - p2 ).SquareModulus() / 4 ) + { + gp_Pnt pProj = projector->Value( uvMid ); + n12 = meshDS->AddNode( pProj.X(), pProj.Y(), pProj.Z() ); + meshDS->SetNodeOnFace( n12, faceID, uvMid.X(), uvMid.Y() ); + myTLinkNodeMap.insert( make_pair ( link, n12 )); + return n12; + } + } uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]); } else if ( pos.second == TopAbs_EDGE ) @@ -1724,7 +1903,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_ if ( !bestEdge.IsNull() ) { - // move n12 to position of a successfull projection + // move n12 to position of a successful projection //double tol = BRep_Tool::Tolerance(edges[ iOkEdge ]); if ( !force3d /*&& distMiddleProj > 2*tol*/ ) { @@ -1934,26 +2113,28 @@ SMDS_MeshFace* SMESH_MesherHelper::AddPolygonalFace (const vectorAddPolygonalFaceWithID(nodes, id); else elem = meshDS->AddPolygonalFace(nodes); } - else { - vector newNodes; - for ( int i = 0; i < nodes.size(); ++i ) + else + { + vector newNodes( nodes.size() * 2 ); + newNodes = nodes; + for ( size_t 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, TopAbs_FACE ); - newNodes.push_back( n1 ); newNodes.push_back( n12 ); } if(id) - elem = meshDS->AddPolygonalFaceWithID(newNodes, id); + elem = meshDS->AddQuadPolygonalFaceWithID(newNodes, id); else - elem = meshDS->AddPolygonalFace(newNodes); + elem = meshDS->AddQuadPolygonalFace(newNodes); } if ( mySetElemOnShape && myShapeID > 0 ) meshDS->SetMeshElementOnShape( elem, myShapeID ); @@ -1995,13 +2176,30 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, const SMDS_MeshNode* n14 = GetMediumNode( n1, n4, force3d, TopAbs_SOLID ); const SMDS_MeshNode* n25 = GetMediumNode( n2, n5, force3d, TopAbs_SOLID ); const SMDS_MeshNode* n36 = GetMediumNode( n3, n6, force3d, TopAbs_SOLID ); + if ( myCreateBiQuadratic ) + { + const SMDS_MeshNode* n1245 = GetCentralNode( n1,n2,n4,n5,n12,n25,n45,n14,force3d ); + const SMDS_MeshNode* n1346 = GetCentralNode( n1,n3,n4,n6,n31,n36,n64,n14,force3d ); + const SMDS_MeshNode* n2356 = GetCentralNode( n2,n3,n6,n5,n23,n36,n56,n25,force3d ); - if(id) - elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, - n12, n23, n31, n45, n56, n64, n14, n25, n36, id); + if(id) + elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36, + n1245, n2356, n1346, id); + else + elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36, + n1245, n2356, n1346); + } else - elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, - n12, n23, n31, n45, n56, n64, n14, n25, n36); + { + if(id) + elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36, id); + else + elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36); + } } if ( mySetElemOnShape && myShapeID > 0 ) meshDS->SetMeshElementOnShape( elem, myShapeID ); @@ -2135,7 +2333,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, const SMDS_MeshNode* n26 = GetMediumNode( n2, n6, force3d, TopAbs_SOLID ); const SMDS_MeshNode* n37 = GetMediumNode( n3, n7, force3d, TopAbs_SOLID ); const SMDS_MeshNode* n48 = GetMediumNode( n4, n8, force3d, TopAbs_SOLID ); - if(myCreateBiQuadratic) + if ( myCreateBiQuadratic ) { const SMDS_MeshNode* n1234 = GetCentralNode( n1,n2,n3,n4,n12,n23,n34,n41,force3d ); const SMDS_MeshNode* n1256 = GetCentralNode( n1,n2,n5,n6,n12,n26,n56,n15,force3d ); @@ -2170,9 +2368,9 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = SMESH_TNodeXYZ( n3478 ); pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = SMESH_TNodeXYZ( n1256 ); - pointsOnShapes[ SMESH_Block::ID_Fx0z ] = SMESH_TNodeXYZ( n1458 ); - pointsOnShapes[ SMESH_Block::ID_Fx1z ] = SMESH_TNodeXYZ( n2367 ); - pointsOnShapes[ SMESH_Block::ID_F0yz ] = SMESH_TNodeXYZ( n1234 ); + pointsOnShapes[ SMESH_Block::ID_Fx0z ] = SMESH_TNodeXYZ( n1458 ); + pointsOnShapes[ SMESH_Block::ID_Fx1z ] = SMESH_TNodeXYZ( n2367 ); + pointsOnShapes[ SMESH_Block::ID_F0yz ] = SMESH_TNodeXYZ( n1234 ); pointsOnShapes[ SMESH_Block::ID_F1yz ] = SMESH_TNodeXYZ( n5678 ); gp_XYZ centerCube(0.5, 0.5, 0.5); @@ -2182,27 +2380,27 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, meshDS->AddNode( nCenterElem.X(), nCenterElem.Y(), nCenterElem.Z() ); meshDS->SetNodeInVolume( nCenter, myShapeID ); - if(id) + if(id) elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, - n12, n23, n34, n41, n56, n67, - n78, n85, n15, n26, n37, n48, - n1234, n1256, n2367, n3478, n1458, n5678, nCenter, id); + n12, n23, n34, n41, n56, n67, + n78, n85, n15, n26, n37, n48, + n1234, n1256, n2367, n3478, n1458, n5678, nCenter, id); else elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8, - n12, n23, n34, n41, n56, n67, - n78, n85, n15, n26, n37, n48, - n1234, n1256, n2367, n3478, n1458, n5678, nCenter); + n12, n23, n34, n41, n56, n67, + n78, n85, n15, n26, n37, n48, + n1234, n1256, n2367, n3478, n1458, n5678, nCenter); } else { if(id) elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, - n12, n23, n34, n41, n56, n67, - n78, n85, n15, n26, n37, n48, id); + n12, n23, n34, n41, n56, n67, + n78, n85, n15, n26, n37, n48, id); else elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8, - n12, n23, n34, n41, n56, n67, - n78, n85, n15, n26, n37, n48); + n12, n23, n34, n41, n56, n67, + n78, n85, n15, n26, n37, n48); } } if ( mySetElemOnShape && myShapeID > 0 ) @@ -2266,7 +2464,7 @@ SMESH_MesherHelper::AddPolyhedralVolume (const std::vector { vector newNodes; vector newQuantities; - for ( int iFace=0, iN=0; iFace < quantities.size(); ++iFace) + for ( size_t iFace = 0, iN = 0; iFace < quantities.size(); ++iFace ) { int nbNodesInFace = quantities[iFace]; newQuantities.push_back(0); @@ -2275,10 +2473,10 @@ SMESH_MesherHelper::AddPolyhedralVolume (const std::vector 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 ) + // if ( n1->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE && + // n2->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE ) { const SMDS_MeshNode* n12 = GetMediumNode( n1, n2, force3d, TopAbs_SOLID ); newNodes.push_back( n12 ); @@ -2401,6 +2599,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 } // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them + const SMDS_MeshNode* prevEndNodes[2] = { 0, 0 }; edge = theBaseSide.begin(); for ( int iE = 0; edge != theBaseSide.end(); ++edge, ++iE ) { @@ -2468,19 +2667,24 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 const double prevPar = theParam2ColumnMap.empty() ? 0 : theParam2ColumnMap.rbegin()->first; for ( u_n = sortedBaseNN.begin(); u_n != sortedBaseNN.end(); u_n++ ) { + if ( u_n->second == prevEndNodes[0] || + u_n->second == prevEndNodes[1] ) + continue; 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 ); } + prevEndNodes[0] = sortedBaseNN.begin()->second; + prevEndNodes[1] = sortedBaseNN.rbegin()->second; } if ( theParam2ColumnMap.size() < 2 ) 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 + size_t prevNbRows = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here + size_t expectNbRows = 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 @@ -2493,10 +2697,10 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 { vector& nCol1 = par_nVec_1->second; vector& nCol2 = par_nVec_2->second; - nCol1.resize( prevNbRows + expectedNbRows ); - nCol2.resize( prevNbRows + expectedNbRows ); + nCol1.resize( prevNbRows + expectNbRows ); + nCol2.resize( prevNbRows + expectNbRows ); - int i1, i2, foundNbRows = 0; + int i1, i2; size_t 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 @@ -2508,7 +2712,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 int nbNodes = face->NbCornerNodes(); if ( nbNodes != 4 ) return false; - if ( foundNbRows + 1 > expectedNbRows ) + if ( foundNbRows + 1 > expectNbRows ) return false; n1 = face->GetNode( (i2+2) % 4 ); // opposite corner of quadrangle face n2 = face->GetNode( (i1+2) % 4 ); @@ -2518,12 +2722,12 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 } avoidSet.insert( face ); } - if ( foundNbRows != expectedNbRows ) + if ((size_t) foundNbRows != expectNbRows ) return false; avoidSet.clear(); } return ( theParam2ColumnMap.size() > 1 && - theParam2ColumnMap.begin()->second.size() == prevNbRows + expectedNbRows ); + theParam2ColumnMap.begin()->second.size() == prevNbRows + expectNbRows ); } namespace @@ -2581,7 +2785,7 @@ bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM ) faceAnalyser.SetSubShape( faceSM->GetSubShape() ); // rotate edges to get the first node being at corner - // (in principle it's not necessary but so far none SALOME algo can make + // (in principle it's not necessary because 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(); @@ -2652,8 +2856,9 @@ bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM ) //purpose : Return true if 2D mesh on FACE is ditorted //======================================================================= -bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, - bool checkUV) +bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, + bool checkUV, + SMESH_MesherHelper* faceHelper) { if ( !faceSM || faceSM->GetSubShape().ShapeType() != TopAbs_FACE ) return false; @@ -2661,12 +2866,26 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, bool haveBadFaces = false; SMESH_MesherHelper helper( *faceSM->GetFather() ); + if ( faceHelper ) + helper.CopySubShapeInfo( *faceHelper ); helper.SetSubShape( faceSM->GetSubShape() ); const TopoDS_Face& F = TopoDS::Face( faceSM->GetSubShape() ); SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( F ); if ( !smDS || smDS->NbElements() == 0 ) return false; + bool subIdsValid = true; // shape ID of nodes is OK + if ( helper.HasSeam() ) + { + // check if nodes are bound to seam edges + SMESH_subMeshIteratorPtr smIt = faceSM->getDependsOnIterator(/*includeSelf=*/false); + while ( smIt->more() && subIdsValid ) + { + SMESH_subMesh* sm = smIt->next(); + if ( helper.IsSeamShape( sm->GetId() ) && sm->IsEmpty() ) + subIdsValid = false; + } + } SMDS_ElemIteratorPtr faceIt = smDS->GetElements(); double prevArea = 0; vector< const SMDS_MeshNode* > nodes; @@ -2682,7 +2901,7 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, for ( size_t i = 0; i < nodes.size(); ++n, ++i ) nodes[ i ] = *n; - // avoid elems on degenarate shapes as UV on them can be wrong + // avoid elems on degenerate shapes as UV on them can be wrong if ( helper.HasDegeneratedEdges() ) { bool isOnDegen = false; @@ -2691,12 +2910,20 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, if ( isOnDegen ) continue; } - // prepare to getting UVs + // prepare for getting UVs const SMDS_MeshNode* inFaceNode = 0; if ( helper.HasSeam() ) { for ( size_t i = 0; ( i < nodes.size() && !inFaceNode ); ++i ) if ( !helper.IsSeamShape( nodes[ i ]->getshapeId() )) + { inFaceNode = nodes[ i ]; + if ( !subIdsValid ) + { + gp_XY uv = helper.GetNodeUV( F, inFaceNode ); + if ( helper.IsOnSeam( uv )) + inFaceNode = NULL; + } + } if ( !inFaceNode ) continue; } @@ -2705,6 +2932,14 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, for ( size_t i = 0; i < nodes.size(); ++i ) uv[ i ] = helper.GetNodeUV( F, nodes[ i ], inFaceNode, toCheckUV ); + if ( !subIdsValid ) // fix uv on seam + { + gp_XY uvInFace = helper.GetNodeUV( F, inFaceNode ); + for ( size_t i = 0; i < uv.size(); ++i ) + if ( helper.IsOnSeam( uv[i] )) + uv[i] = helper.getUVOnSeam( uv[i], uvInFace ).XY(); + } + // compare orientation of triangles double faceArea = 0; for ( int iT = 0, nbT = nodes.size()-2; iT < nbT; ++iT ) @@ -2725,7 +2960,7 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, * \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 + * in the corresponding submesh point in different directions */ //================================================================================ @@ -2743,38 +2978,90 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace) if ( !aSubMeshDSFace ) return isReversed; - // find an element with a good normal - gp_Vec Ne; - bool normalOK = false; - gp_XY uv; + // find an element on a bounday of theFace SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); - while ( !normalOK && iteratorElem->more() ) // loop on elements on theFace + const SMDS_MeshNode* nn[2]; + while ( iteratorElem->more() ) // loop on elements on theFace { const SMDS_MeshElement* elem = iteratorElem->next(); - if ( elem && elem->NbCornerNodes() > 2 ) + if ( ! elem ) continue; + + // look for 2 nodes on EDGE + int nbNodes = elem->NbCornerNodes(); + nn[0] = elem->GetNode( nbNodes-1 ); + for ( int iN = 0; iN < nbNodes; ++iN ) { - SMESH_TNodeXYZ nPnt[3]; - SMDS_ElemIteratorPtr nodesIt = elem->nodesIterator(); - int iNodeOnFace = 0, iPosDim = SMDS_TOP_VERTEX; - for ( int iN = 0; nodesIt->more() && iN < 3; ++iN) // loop on nodes + nn[1] = elem->GetNode( iN ); + if ( nn[0]->GetPosition()->GetDim() < 2 && + nn[1]->GetPosition()->GetDim() < 2 ) { - nPnt[ iN ] = nodesIt->next(); - if ( nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition() > iPosDim ) + TopoDS_Shape s0 = GetSubShapeByNode( nn[0], GetMeshDS() ); + TopoDS_Shape s1 = GetSubShapeByNode( nn[1], GetMeshDS() ); + TopoDS_Shape E = GetCommonAncestor( s0, s1, *myMesh, TopAbs_EDGE ); + if ( !E.IsNull() && !s0.IsSame( s1 ) && E.Orientation() != TopAbs_INTERNAL ) { - iNodeOnFace = iN; - iPosDim = nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition(); + // is E seam edge? + int nb = 0; + for ( TopExp_Explorer exp( theFace, TopAbs_EDGE ); exp.More(); exp.Next() ) + if ( E.IsSame( exp.Current() )) { + ++nb; + E = exp.Current(); // to know orientation + } + if ( nb == 1 ) + { + bool ok = true; + double u0 = GetNodeU( TopoDS::Edge( E ), nn[0], nn[1], &ok ); + double u1 = GetNodeU( TopoDS::Edge( E ), nn[1], nn[0], &ok ); + if ( ok ) + { + // check that the 2 nodes are connected with a segment (IPAL53055) + ok = false; + const SMDS_MeshElement* seg; + if ( SMESHDS_SubMesh* sm = GetMeshDS()->MeshElements( E )) + if (( sm->NbElements() > 0 ) && + ( seg = GetMeshDS()->FindEdge( nn[0], nn[1] ))) + ok = sm->Contains( seg ); + } + if ( ok ) + { + isReversed = ( u0 > u1 ); + if ( E.Orientation() == TopAbs_REVERSED ) + isReversed = !isReversed; + return isReversed; + } + } } } - // compute normal - gp_Vec v01( nPnt[0], nPnt[1] ), v02( nPnt[0], nPnt[2] ); - if ( v01.SquareMagnitude() > RealSmall() && - v02.SquareMagnitude() > RealSmall() ) + nn[0] = nn[1]; + } + } + + // find an element with a good normal + gp_Vec Ne; + bool normalOK = false; + gp_XY uv; + iteratorElem = aSubMeshDSFace->GetElements(); + while ( !normalOK && iteratorElem->more() ) // loop on elements on theFace + { + const SMDS_MeshElement* elem = iteratorElem->next(); + if ( ! SMESH_MeshAlgos::FaceNormal( elem, const_cast( Ne.XYZ() ), /*normalized=*/0 )) + continue; + normalOK = true; + + // get UV of a node inside theFACE + SMDS_ElemIteratorPtr nodesIt = elem->nodesIterator(); + const SMDS_MeshNode* nInFace = 0; + int iPosDim = SMDS_TOP_VERTEX; + while ( nodesIt->more() ) // loop on nodes + { + const SMDS_MeshNode* n = static_cast( nodesIt->next() ); + if ( n->GetPosition()->GetTypeOfPosition() >= iPosDim ) { - Ne = v01 ^ v02; - if (( normalOK = ( Ne.SquareMagnitude() > RealSmall() ))) - uv = GetNodeUV( theFace, nPnt[iNodeOnFace]._node, 0, &normalOK ); + nInFace = n; + iPosDim = n->GetPosition()->GetTypeOfPosition(); } } + uv = GetNodeUV( theFace, nInFace, 0, &normalOK ); } if ( !normalOK ) return isReversed; @@ -2785,11 +3072,8 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace) // 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; - } + 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 ); @@ -2865,7 +3149,7 @@ TopAbs_Orientation SMESH_MesherHelper::GetSubShapeOri(const TopoDS_Shape& shape, //======================================================================= //function : IsSubShape -//purpose : +//purpose : //======================================================================= bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, @@ -2879,8 +3163,6 @@ bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, if ( shape.IsSame( exp.Current() )) return true; } - SCRUTE((shape.IsNull())); - SCRUTE((mainShape.IsNull())); return false; } @@ -2965,7 +3247,7 @@ double SMESH_MesherHelper::getFaceMaxTol( const TopoDS_Shape& face ) const * of the FACE normal * \return double - the angle (between -Pi and Pi), negative if the angle is concave, * 1e100 in case of failure - * \waring Care about order of the EDGEs and their orientation to be as they are + * \warning Care about order of the EDGEs and their orientation to be as they are * within the FACE! Don't pass degenerated EDGEs neither! */ //================================================================================ @@ -3106,6 +3388,28 @@ TopAbs_ShapeEnum SMESH_MesherHelper::GetGroupType(const TopoDS_Shape& group, return TopAbs_SHAPE; } +//================================================================================ +/*! + * \brief Returns a shape, to which a hypothesis used to mesh a given shape is assigned + * \param [in] hyp - the hypothesis + * \param [in] shape - the shape, for meshing which the \a hyp is used + * \param [in] mesh - the mesh + * \return TopoDS_Shape - the shape the \a hyp is assigned to + */ +//================================================================================ + +TopoDS_Shape SMESH_MesherHelper::GetShapeOfHypothesis( const SMESHDS_Hypothesis * hyp, + const TopoDS_Shape& shape, + SMESH_Mesh* mesh) +{ + const SMESH_Hypothesis* h = static_cast( hyp ); + SMESH_HypoFilter hypFilter( SMESH_HypoFilter::Is( h )); + + TopoDS_Shape shapeOfHyp; + mesh->GetHypothesis( shape, hypFilter, /*checkAncestors=*/true, &shapeOfHyp ); + return shapeOfHyp; +} + //======================================================================= //function : IsQuadraticMesh //purpose : Check mesh without geometry for: if all elements on this shape are quadratic, @@ -3153,6 +3457,43 @@ double SMESH_MesherHelper::GetOtherParam(const double param) const return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i]; } +//======================================================================= +//function : NbRealSeam +//purpose : Return a number of real seam edges in the shape set through +// IsQuadraticSubMesh() or SetSubShape(). A real seam edge encounters twice in a wire +//======================================================================= + +size_t SMESH_MesherHelper::NbRealSeam() const +{ + size_t nb = 0; + + std::set< int >::const_iterator id = mySeamShapeIds.begin(); + for ( ; id != mySeamShapeIds.end(); ++id ) + if ( *id < 0 ) ++nb; + else break; + + return nb; +} + +//======================================================================= +//function : IsOnSeam +//purpose : Check if UV is on seam. Return 0 if not, 1 for U seam, 2 for V seam +//======================================================================= + +int SMESH_MesherHelper::IsOnSeam(const gp_XY& uv) const +{ + for ( int i = U_periodic; i <= V_periodic ; ++i ) + if ( myParIndex & i ) + { + double p = uv.Coord( i ); + double tol = ( myPar2[i-1] - myPar1[i-1] ) / 100.; + if ( Abs( p - myPar1[i-1] ) < tol || + Abs( p - myPar2[i-1] ) < tol ) + return i; + } + return 0; +} + namespace { //======================================================================= @@ -3166,11 +3507,16 @@ namespace { TopTools_ListIteratorOfListOfShape _ancIter; TopAbs_ShapeEnum _type; TopTools_MapOfShape _encountered; - TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type) + TopTools_IndexedMapOfShape _allowed; + TAncestorsIterator( const TopTools_ListOfShape& ancestors, + TopAbs_ShapeEnum type, + const TopoDS_Shape* container/* = 0*/) : _ancIter( ancestors ), _type( type ) { + if ( container && !container->IsNull() ) + TopExp::MapShapes( *container, type, _allowed); if ( _ancIter.More() ) { - if ( _ancIter.Value().ShapeType() != _type ) next(); + if ( !isCurrentAllowed() ) next(); else _encountered.Add( _ancIter.Value() ); } } @@ -3183,25 +3529,32 @@ namespace { 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() )) + if ( isCurrentAllowed() && _encountered.Add( _ancIter.Value() )) break; return s; } + bool isCurrentAllowed() + { + return (( _ancIter.Value().ShapeType() == _type ) && + ( _allowed.IsEmpty() || _allowed.Contains( _ancIter.Value() ))); + } }; } // namespace //======================================================================= /*! - * \brief Return iterator on ancestors of the given type + * \brief Return iterator on ancestors of the given type, included into a container shape */ //======================================================================= PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape, const SMESH_Mesh& mesh, - TopAbs_ShapeEnum ancestorType) + TopAbs_ShapeEnum ancestorType, + const TopoDS_Shape* container) { - return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType)); + return PShapeIteratorPtr + ( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType, container)); } //======================================================================= @@ -3217,6 +3570,11 @@ TopoDS_Shape SMESH_MesherHelper::GetCommonAncestor(const TopoDS_Shape& shape1, TopoDS_Shape commonAnc; if ( !shape1.IsNull() && !shape2.IsNull() ) { + if ( shape1.ShapeType() == ancestorType && IsSubShape( shape2, shape1 )) + return shape1; + if ( shape2.ShapeType() == ancestorType && IsSubShape( shape1, shape2 )) + return shape2; + PShapeIteratorPtr ancIt = GetAncestors( shape1, mesh, ancestorType ); while ( const TopoDS_Shape* anc = ancIt->next() ) if ( IsSubShape( shape2, *anc )) @@ -3257,30 +3615,32 @@ namespace { // Structures used by FixQuadraticElements() mutable vector _faces; mutable gp_Vec _nodeMove; mutable int _nbMoves; + mutable bool _is2dFixed; // is moved along surface or in 3D QLink(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshNode* nm): SMESH_TLink( n1,n2 ), _mediumNode(nm), _nodeMove(0,0,0), _nbMoves(0) { _faces.reserve(4); - //if ( MediumPos() != SMDS_TOP_3DSPACE ) - _nodeMove = MediumPnt() - MiddlePnt(); + _nodeMove = MediumPnt() - MiddlePnt(); + _is2dFixed = ( MediumPos() != SMDS_TOP_FACE ); } void SetContinuesFaces() const; const QFace* GetContinuesFace( const QFace* face ) const; - bool OnBoundary() const; + bool OnBoundary() const; gp_XYZ MiddlePnt() const { return ( XYZ( node1() ) + XYZ( node2() )) / 2.; } gp_XYZ MediumPnt() const { return XYZ( _mediumNode ); } - SMDS_TypeOfPosition MediumPos() const + SMDS_TypeOfPosition MediumPos() const { return _mediumNode->GetPosition()->GetTypeOfPosition(); } - SMDS_TypeOfPosition EndPos(bool isSecond) const + SMDS_TypeOfPosition EndPos(bool isSecond) const { return (isSecond ? node2() : node1())->GetPosition()->GetTypeOfPosition(); } const SMDS_MeshNode* EndPosNode(SMDS_TypeOfPosition pos) const { return EndPos(0) == pos ? node1() : EndPos(1) == pos ? node2() : 0; } - void Move(const gp_Vec& move, bool sum=false) const - { _nodeMove += move; _nbMoves += sum ? (_nbMoves==0) : 1; } + void Move(const gp_Vec& move, bool sum=false, bool is2dFixed=false) const + { _nodeMove += move; _nbMoves += sum ? (_nbMoves==0) : 1; _is2dFixed |= is2dFixed; } gp_XYZ Move() const { return _nodeMove.XYZ() / _nbMoves; } bool IsMoved() const { return (_nbMoves > 0 /*&& !IsStraight()*/); } + bool IsFixedOnSurface() const { return _is2dFixed; } bool IsStraight() const { return isStraightLink( (XYZ(node1())-XYZ(node2())).SquareModulus(), _nodeMove.SquareMagnitude()); @@ -3358,11 +3718,11 @@ namespace { // Structures used by FixQuadraticElements() int NbVolumes() const { return !_volumes[0] ? 0 : !_volumes[1] ? 1 : 2; } void AddSelfToLinks() const { - for ( int i = 0; i < _sides.size(); ++i ) + for ( size_t i = 0; i < _sides.size(); ++i ) _sides[i]->_faces.push_back( this ); } int LinkIndex( const QLink* side ) const { - for (int i=0; i<_sides.size(); ++i ) if ( _sides[i] == side ) return i; + for (size_t i = 0; i<_sides.size(); ++i ) if ( _sides[i] == side ) return i; return -1; } bool GetLinkChain( int iSide, TChain& chain, SMDS_TypeOfPosition pos, int& err) const; @@ -3394,7 +3754,7 @@ namespace { // Structures used by FixQuadraticElements() const SMDS_MeshNode* nodeToContain) const; const SMDS_MeshNode* GetNodeInFace() const { - for ( int iL = 0; iL < _sides.size(); ++iL ) + for ( size_t iL = 0; iL < _sides.size(); ++iL ) if ( _sides[iL]->MediumPos() == SMDS_TOP_FACE ) return _sides[iL]->_mediumNode; return 0; } @@ -3447,7 +3807,7 @@ namespace { // Structures used by FixQuadraticElements() _sides = links; _sideIsAdded[0]=_sideIsAdded[1]=_sideIsAdded[2]=_sideIsAdded[3]=false; _normal.SetCoord(0,0,0); - for ( int i = 1; i < _sides.size(); ++i ) { + for ( size_t i = 1; i < _sides.size(); ++i ) { const QLink *l1 = _sides[i-1], *l2 = _sides[i]; insert( l1->node1() ); insert( l1->node2() ); // compute normal @@ -3472,7 +3832,7 @@ namespace { // Structures used by FixQuadraticElements() * \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 + * \param pos - position of medium nodes the links should have * \param error - out, specifies what is wrong * \retval bool - false if valid chain can't be built; "valid" means that links * of the chain belongs to rectangles bounding hexahedrons @@ -3481,18 +3841,18 @@ namespace { // Structures used by FixQuadraticElements() bool QFace::GetLinkChain( int iSide, TChain& chain, SMDS_TypeOfPosition pos, int& error) const { - if ( iSide >= _sides.size() ) // wrong argument iSide + if ( iSide >= (int)_sides.size() ) // wrong argument iSide return false; if ( _sideIsAdded[ iSide ]) // already in chain return true; - if ( _sides.size() != 4 ) { // triangle - visit all my continous faces + if ( _sides.size() != 4 ) { // triangle - visit all my continuous faces MSGBEG( *this ); TLinkSet links; list< const QFace* > faces( 1, this ); while ( !faces.empty() ) { const QFace* face = faces.front(); - for ( int i = 0; i < face->_sides.size(); ++i ) { + for ( size_t i = 0; i < face->_sides.size(); ++i ) { if ( !face->_sideIsAdded[i] && face->_sides[i] ) { face->_sideIsAdded[i] = true; // find a face side in the chain @@ -3532,7 +3892,7 @@ namespace { // Structures used by FixQuadraticElements() if ( link->MediumPos() >= pos ) { int nbLinkFaces = link->_faces.size(); if ( nbLinkFaces == 4 || (/*nbLinkFaces < 4 && */link->OnBoundary())) { - // hexahedral mesh or boundary quadrangles - goto a continous face + // hexahedral mesh or boundary quadrangles - goto a continuous face if ( const QFace* f = link->GetContinuesFace( this )) if ( f->_sides.size() == 4 ) return f->GetLinkChain( *chLink, chain, pos, error ); @@ -3575,7 +3935,7 @@ namespace { // Structures used by FixQuadraticElements() typedef list< pair< const QFace*, TLinkInSet > > TFaceLinkList; TFaceLinkList adjacentFaces; - for ( int iL = 0; iL < _sides.size(); ++iL ) + for ( size_t iL = 0; iL < _sides.size(); ++iL ) { if ( avoidLink._qlink == _sides[iL] ) continue; @@ -3628,10 +3988,10 @@ namespace { // Structures used by FixQuadraticElements() const TChainLink& avoidLink, const SMDS_MeshNode* nodeToContain) const { - for ( int i = 0; i < _sides.size(); ++i ) + for ( size_t i = 0; i < _sides.size(); ++i ) if ( avoidLink._qlink != _sides[i] && (_sides[i]->node1() == nodeToContain || _sides[i]->node2() == nodeToContain )) - return links.find( _sides[ i ]); + return links.find( _sides[i] ); return links.end(); } @@ -3660,7 +4020,7 @@ namespace { // Structures used by FixQuadraticElements() * \brief Move medium node of theLink according to its distance from boundary * \param theLink - link to fix * \param theRefVec - movement of boundary - * \param theLinks - all adjacent links of continous triangles + * \param theLinks - all adjacent links of continuous triangles * \param theFaceHelper - helper is not used so far * \param thePrevLen - distance from the boundary * \param theStep - number of steps till movement propagation limit @@ -3682,7 +4042,7 @@ namespace { // Structures used by FixQuadraticElements() if ( !theStep ) return thePrevLen; // propagation limit reached - int iL; // index of theLink + size_t iL; // index of theLink for ( iL = 0; iL < _sides.size(); ++iL ) if ( theLink._qlink == _sides[ iL ]) break; @@ -3759,7 +4119,7 @@ namespace { // Structures used by FixQuadraticElements() double r = thePrevLen / fullLen; gp_Vec move = linkNorm * refProj * ( 1 - r ); - theLink->Move( move, true ); + theLink->Move( move, /*sum=*/true ); MSG(string(theStep,'.')<<" Move "<< theLink->_mediumNode->GetID()<< " by " << refProj * ( 1 - r ) << " following " << @@ -3780,9 +4140,9 @@ namespace { // Structures used by FixQuadraticElements() { // 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 ) + for ( TIDSortedNodeSet::const_iterator n = begin(); n != end(); ++n ) + gc += XYZ( *n ) / double( size() ); + for ( size_t 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())); @@ -3798,12 +4158,12 @@ namespace { // Structures used by FixQuadraticElements() return true; } return false; - + } //================================================================================ /*! - * \brief Find pairs of continues faces + * \brief Find pairs of continues faces */ //================================================================================ @@ -3814,7 +4174,7 @@ namespace { // Structures used by FixQuadraticElements() // | Between _faces of link x2 two vertical faces are continues // x1----x2-----x3 and two horizontal faces are continues. We set vertical faces // | to _faces[0] and _faces[1] and horizontal faces to - // v2 | v3 _faces[2] and _faces[3] (or vise versa). + // v2 | v3 _faces[2] and _faces[3] (or vice versa). // x4 if ( _faces.empty() ) @@ -3822,7 +4182,7 @@ namespace { // Structures used by FixQuadraticElements() 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 ) + for ( size_t iF = 1; iFaceCont < 0 && iF < _faces.size(); ++iF ) { // look for a face bounding none of volumes bound by _faces[0] bool sameVol = false; @@ -3864,12 +4224,13 @@ namespace { // Structures used by FixQuadraticElements() const QFace* QLink::GetContinuesFace( const QFace* face ) const { - for ( int i = 0; i < _faces.size(); ++i ) { - if ( _faces[i] == face ) { - int iF = i < 2 ? 1-i : 5-i; - return iF < _faces.size() ? _faces[iF] : 0; + if ( _faces.size() <= 4 ) + for ( size_t i = 0; i < _faces.size(); ++i ) { + if ( _faces[i] == face ) { + int iF = i < 2 ? 1-i : 5-i; + return iF < (int)_faces.size() ? _faces[iF] : 0; + } } - } return 0; } //================================================================================ @@ -3880,7 +4241,7 @@ namespace { // Structures used by FixQuadraticElements() bool QLink::OnBoundary() const { - for ( int i = 0; i < _faces.size(); ++i ) + for ( size_t i = 0; i < _faces.size(); ++i ) if (_faces[i] && _faces[i]->IsBoundary()) return true; return false; } @@ -3949,7 +4310,7 @@ namespace { // Structures used by FixQuadraticElements() for ( ; bnd != bndEnd; ++bnd ) { const QLink* bndLink = *bnd; - for ( int i = 0; i < bndLink->_faces.size(); ++i ) // loop on faces of bndLink + for ( size_t i = 0; i < bndLink->_faces.size(); ++i ) // loop on faces of bndLink { const QFace* face = bndLink->_faces[i]; // quadrange lateral face of a prism if ( !face ) continue; @@ -4014,9 +4375,9 @@ namespace { // Structures used by FixQuadraticElements() vector< TChain> & resultChains, SMDS_TypeOfPosition pos ) { - // put links in the set and evalute number of result chains by number of boundary links + // put links in the set and evaluate number of result chains by number of boundary links TLinkSet linkSet; - int nbBndLinks = 0; + size_t nbBndLinks = 0; for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) { linkSet.insert( *lnk ); nbBndLinks += lnk->IsBoundary(); @@ -4065,7 +4426,7 @@ namespace { // Structures used by FixQuadraticElements() TLinkInSet botLink = startLink; // current horizontal link to go up from corner = startCorner; // current corner the botLink ends at - int iRow = 0; + size_t iRow = 0; while ( botLink != linksEnd ) // loop on rows { // add botLink to the columnChain @@ -4162,7 +4523,7 @@ namespace { // Structures used by FixQuadraticElements() // In the linkSet, there must remain the last links of rowChains; add them if ( linkSet.size() != rowChains.size() ) return _BAD_SET_SIZE; - for ( int iRow = 0; iRow < rowChains.size(); ++iRow ) { + for ( size_t iRow = 0; iRow < rowChains.size(); ++iRow ) { // find the link (startLink) ending at startCorner corner = 0; for ( startLink = linkSet.begin(); startLink != linksEnd; ++startLink ) { @@ -4198,8 +4559,16 @@ namespace { // Structures used by FixQuadraticElements() TopoDS_Shape shape = theHelper.GetSubShape().Oriented( TopAbs_FORWARD ); if ( shape.IsNull() ) return; - if ( !theError ) theError = SMESH_ComputeError::New(); - + if ( !dynamic_cast( theError.get() )) + { + if ( !theError ) + theError.reset( new SMESH_BadInputElements( meshDS )); + else + theError.reset( new SMESH_BadInputElements( meshDS, + theError->myName, + theError->myComment, + theError->myAlgo)); + } gp_XYZ faceNorm; if ( shape.ShapeType() == TopAbs_FACE ) // 2D @@ -4254,6 +4623,7 @@ namespace { // Structures used by FixQuadraticElements() { continue; } + default:; } // get nodes shared by faces that may be distorted SMDS_NodeIteratorPtr nodeIt; @@ -4316,13 +4686,13 @@ namespace { // Structures used by FixQuadraticElements() 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 ); + static_cast( theError.get() )->add( f ); } } } } } - if ( !theError->myBadElements.empty() ) + if ( theError->HasBadElems() ) theError->myName = EDITERR_NO_MEDIUM_ON_GEOM; return; @@ -4367,6 +4737,7 @@ namespace { // Structures used by FixQuadraticElements() { concaveFaces.push_back( face ); } + default:; } } if ( concaveFaces.empty() ) @@ -4388,8 +4759,8 @@ namespace { // Structures used by FixQuadraticElements() < 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 + // search to check if a volume is close to a concave face + SMESHUtils::Deleter< SMESH_ElementSearcher > faceSearcher ( SMESH_MeshAlgos::GetElementSearcher( *theHelper.GetMeshDS(), faceIter )); // classifier @@ -4432,7 +4803,7 @@ namespace { // Structures used by FixQuadraticElements() while ( volIt->more() ) { const SMDS_MeshElement* vol = volIt->next(); - int nbN = vol->NbCornerNodes(); + size_t nbN = vol->NbCornerNodes(); if ( ( nbN != 4 && nbN != 5 ) || !solidSM->Contains( vol ) || !checkedVols.insert( vol ).second ) @@ -4491,7 +4862,7 @@ namespace { // Structures used by FixQuadraticElements() 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 )); + isDistorted = ( Abs( hMedium ) > Abs( hVol * 0.75 )); } } } @@ -4509,13 +4880,13 @@ namespace { // Structures used by FixQuadraticElements() MSG( "move OUT of solid " << nMedium ); } } - theError->myBadElements.push_back( vol ); + static_cast( theError.get() )->add( vol ); } } // loop on volumes sharing a node on FACE } // loop on nodes on FACE } // loop on FACEs of a SOLID - if ( !theError->myBadElements.empty() ) + if ( theError->HasBadElems() ) theError->myName = EDITERR_NO_MEDIUM_ON_GEOM; } // 3D case } @@ -4527,7 +4898,7 @@ namespace { // Structures used by FixQuadraticElements() * \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 */ //======================================================================= @@ -4535,6 +4906,7 @@ namespace { // Structures used by FixQuadraticElements() void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, bool volumeOnly) { + //MESSAGE("FixQuadraticElements " << volumeOnly); // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion if ( getenv("NO_FixQuadraticElements") ) return; @@ -4572,9 +4944,11 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, } // fix nodes on geom faces #ifdef _DEBUG_ - int nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; + int nbfaces = nbSolids; + nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; #endif for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) { + MESSAGE("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key())); MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key())); SMESH_MesherHelper h(*myMesh); h.SetSubShape( fIt.Key() ); @@ -4750,7 +5124,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, } else if ( error == ERR_TRI ) { // chain contains continues triangles TSplitTriaResult res = splitTrianglesIntoChains( rawChain, chains, pos ); - if ( res != _OK ) { // not quadrangles split into triangles + if ( res != _OK ) { // not 'quadrangles split into triangles' in chain fixTriaNearBoundary( rawChain, *this ); break; } @@ -4762,7 +5136,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, else { continue; } - for ( int iC = 0; iC < chains.size(); ++iC ) + for ( size_t iC = 0; iC < chains.size(); ++iC ) { TChain& chain = chains[iC]; if ( chain.empty() ) continue; @@ -4775,15 +5149,17 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, MSG("Internal chain - ignore"); continue; } - // mesure chain length and compute link position along the chain + // measure chain length and compute link position along the chain double chainLen = 0; vector< double > linkPos; + TChain savedChain; // backup MSGBEG( "Link medium nodes: "); TChain::iterator link0 = chain.begin(), link1 = chain.begin(), link2; for ( ++link1; link1 != chain.end(); ++link1, ++link0 ) { MSGBEG( (*link0)->_mediumNode->GetID() << "-" <<(*link1)->_mediumNode->GetID()<<" "); double len = ((*link0)->MiddlePnt() - (*link1)->MiddlePnt()).Modulus(); while ( len < numeric_limits::min() ) { // remove degenerated link + if ( savedChain.empty() ) savedChain = chain; link1 = chain.erase( link1 ); if ( link1 == chain.end() ) break; @@ -4793,9 +5169,16 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, linkPos.push_back( chainLen ); } MSG(""); - if ( linkPos.size() < 2 ) - continue; - + if ( linkPos.size() <= 2 && savedChain.size() > 2 ) { + //continue; + linkPos.clear(); + chainLen = 0; + chain = savedChain; + for ( link1 = chain.begin(); link1 != chain.end(); ++link1 ) { + chainLen += 1; + linkPos.push_back( chainLen ); + } + } gp_Vec move0 = chain.front()->_nodeMove; gp_Vec move1 = chain.back ()->_nodeMove; @@ -4811,7 +5194,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, face = TopoDS::Face( f ); faceHlp.SetSubShape( face ); Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc); - bool isStraight[2]; + //bool isStraight[2]; // commented for issue 0023118 for ( int is1 = 0; is1 < 2; ++is1 ) // move0 or move1 { TChainLink& link = is1 ? chain.back() : chain.front(); @@ -4820,17 +5203,17 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, gp_XY uv2 = faceHlp.GetNodeUV( face, link->node2(), nodeOnFace, &checkUV ); gp_XY uv12 = faceHlp.GetMiddleUV( surf, uv1, uv2 ); // uvMove = uvm - uv12 - gp_XY uvMove = applyIn2D(surf, uvm, uv12, gp_XY_Subtracted, /*inPeriod=*/false); + 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 ( isStraight[0] && isStraight[1] ) { - MSG("2D straight - ignore"); - continue; // straight - no need to move nodes of internal links + // isStraight[is1] = isStraightLink( (uv2-uv1).SquareModulus(), + // 10 * uvMove.SquareModulus()); } + // 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 = faceHlp.GetNodeUV( face, linkOnFace->_mediumNode, 0, &checkUV ); @@ -4874,19 +5257,23 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // transform to global gp_Vec x01( (*link0)->MiddlePnt(), (*link1)->MiddlePnt() ); gp_Vec x12( (*link1)->MiddlePnt(), (*link2)->MiddlePnt() ); - gp_Vec x = x01.Normalized() + x12.Normalized(); - trsf.SetTransformation( gp_Ax3( gp::Origin(), link1->Normal(), x), gp_Ax3() ); + try { + gp_Vec x = x01.Normalized() + x12.Normalized(); + trsf.SetTransformation( gp_Ax3( gp::Origin(), link1->Normal(), x), gp_Ax3() ); + } catch ( Standard_Failure ) { + trsf.Invert(); + } move.Transform(trsf); + (*link1)->Move( move, /*sum=*/false, /*is2dFixed=*/false ); } else { // compute 3D displacement by 2D one Handle(Geom_Surface) s = BRep_Tool::Surface(face,loc); gp_XY oldUV = faceHlp.GetNodeUV( face, (*link1)->_mediumNode, 0, &checkUV ); - gp_XY newUV = applyIn2D( s, oldUV, gp_XY( move.X(),move.Y()), gp_XY_Added ); + gp_XY newUV = ApplyIn2D( s, oldUV, gp_XY( move.X(),move.Y()), gp_XY_Added ); gp_Pnt newPnt = s->Value( newUV.X(), newUV.Y()); move = gp_Vec( XYZ((*link1)->_mediumNode), newPnt.Transformed(loc) ); - if ( SMDS_FacePosition* nPos = - dynamic_cast< SMDS_FacePosition* >((*link1)->_mediumNode->GetPosition())) + if ( SMDS_FacePositionPtr nPos = (*link1)->_mediumNode->GetPosition()) nPos->SetParameters( newUV.X(), newUV.Y() ); #ifdef _DEBUG_ if ( (XYZ((*link1)->node1()) - XYZ((*link1)->node2())).SquareModulus() < @@ -4899,10 +5286,11 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, "uv2: "<Move( move, /*sum=*/false, /*is2dFixed=*/true ); } - (*link1)->Move( move ); MSG( "Move " << (*link1)->_mediumNode->GetID() << " following " << chain.front()->_mediumNode->GetID() <<"-" << chain.back ()->_mediumNode->GetID() << @@ -4916,22 +5304,38 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // 4. Move nodes // ------------- - TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa; - const SMDS_MeshElement *biQuadQua, *triQuadHex; + TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa, biQuadPenta; const bool toFixCentralNodes = ( myMesh->NbBiQuadQuadrangles() + myMesh->NbBiQuadTriangles() + - myMesh->NbTriQuadraticHexas() ); + myMesh->NbTriQuadraticHexas() + + myMesh->NbBiQuadPrisms()); + double distXYZ[4]; + faceHlp.ToFixNodeParameters( true ); for ( pLink = links.begin(); pLink != links.end(); ++pLink ) { if ( pLink->IsMoved() ) { gp_Pnt p = pLink->MiddlePnt() + pLink->Move(); + + // put on surface nodes on FACE but moved in 3D (23050) + if ( !pLink->IsFixedOnSurface() ) + { + faceHlp.SetSubShape( pLink->_mediumNode->getshapeId() ); + if ( faceHlp.GetSubShape().ShapeType() == TopAbs_FACE ) + { + const_cast( pLink->_mediumNode )->setXYZ( p.X(), p.Y(), p.Z()); + p.Coord( distXYZ[1], distXYZ[2], distXYZ[3] ); + gp_XY uv( Precision::Infinite(), 0 ); + if ( faceHlp.CheckNodeUV( TopoDS::Face( faceHlp.GetSubShape() ), pLink->_mediumNode, + uv, /*tol=*/pLink->Move().Modulus(), /*force=*/true, distXYZ )) + p.SetCoord( distXYZ[1], distXYZ[2], distXYZ[3] ); + } + } GetMeshDS()->MoveNode( pLink->_mediumNode, p.X(), p.Y(), p.Z()); // collect bi-quadratic elements if ( toFixCentralNodes ) { - biQuadQua = triQuadHex = 0; SMDS_ElemIteratorPtr eIt = pLink->_mediumNode->GetInverseElementIterator(); while ( eIt->more() ) { @@ -4940,6 +5344,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, case SMDSEntity_BiQuad_Quadrangle: biQuadQuas.insert( e ); break; case SMDSEntity_BiQuad_Triangle: biQuadTris.insert( e ); break; case SMDSEntity_TriQuad_Hexa: triQuadHexa.insert( e ); break; + case SMDSEntity_BiQuad_Penta: biQuadPenta.insert( e ); break; default:; } } @@ -4974,6 +5379,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, if ( i > 3 && nodes[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) CheckNodeUV( F, nodes[i], uv[ i ], 2*tol, /*force=*/true ); } + AdjustByPeriod( F, uv, 8 ); // put uv[] within a period (IPAL52698) // move the central node gp_XY uvCent = calcTFI (0.5, 0.5, uv[0],uv[1],uv[2],uv[3],uv[4],uv[5],uv[6],uv[7] ); gp_Pnt p = surf->Value( uvCent.X(), uvCent.Y() ).Transformed( loc ); @@ -4999,7 +5405,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // nodes nodes.assign( tria->begin_nodes(), tria->end_nodes() ); // UV - bool uvOK = true, badTria; + bool uvOK = true, badTria = false; for ( int i = 0; i < 6; ++i ) { uv[ i ] = GetNodeUV( F, nodes[i], nodes[(i+1)%3], &uvOK ); @@ -5008,6 +5414,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, if ( nodes[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) CheckNodeUV( F, nodes[i], uv[ i ], 2*tol, /*force=*/true ); } + // move the central node gp_Pnt p; if ( !uvOK || badTria ) @@ -5018,6 +5425,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, } else { + AdjustByPeriod( F, uv, 6 ); // put uv[] within a period (IPAL52698) gp_XY uvCent = GetCenterUV( uv[0], uv[1], uv[2], uv[3], uv[4], uv[5], &badTria ); p = surf->Value( uvCent.X(), uvCent.Y() ).Transformed( loc ); } @@ -5087,4 +5495,34 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, nCenterCoords.X(), nCenterCoords.Y(), nCenterCoords.Z()); } } + // treat tri-quadratic hexahedra + { + SMDS_VolumeTool volExp; + TIDSortedElemSet::iterator pentIt = biQuadPenta.begin(); + for ( ; pentIt != biQuadPenta.end(); ++pentIt ) + { + MESSAGE("---"); + volExp.Set( *pentIt, /*ignoreCentralNodes=*/false ); + } + } +#ifdef _DEBUG_ + // avoid warning: defined but not used operator<<() + SMESH_Comment() << *links.begin() << *faces.begin(); +#endif } + +//================================================================================ +/*! + * \brief DEBUG + */ +//================================================================================ + +void SMESH_MesherHelper::WriteShape(const TopoDS_Shape& s) +{ + const char* name = "/tmp/shape.brep"; + BRepTools::Write( s, name ); +#ifdef _DEBUG_ + std::cout << name << std::endl; +#endif +} +