X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=6ed52110d8303f0353d0be27cc2837d2b27c6173;hp=2ba354d20866b55266980d59fa561136e96bda5b;hb=44cfe3f661cf16b7253c2635368807ea332a4cc5;hpb=e29c1bb636a70d56d65c28211dfeaaae4654b3d0 diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index 2ba354d20..6ed52110d 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -32,6 +32,7 @@ #include "SMDS_IteratorOnIterators.hxx" #include "SMDS_VolumeTool.hxx" #include "SMESH_Block.hxx" +#include "SMESH_HypoFilter.hxx" #include "SMESH_MeshAlgos.hxx" #include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" @@ -190,12 +191,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; @@ -647,10 +649,10 @@ 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 ) { + if ( !uvOK ) + { + if ( !IsSubShape( V, F )) + { MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID << " not in face " << GetMeshDS()->ShapeToIndex( F ) ); // get UV of a vertex closest to the node @@ -667,7 +669,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 +679,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 )) @@ -805,6 +818,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 +844,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 +885,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 +920,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 ); } //======================================================================= @@ -1266,14 +1324,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 == 4 ) + shapeType = TopAbs_FACE; + faceID = myShapeID; + } + else + { + 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,11 +1359,20 @@ 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)); + 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] ); + TopLoc_Location loc; Handle( Geom_Surface ) S = BRep_Tool::Surface( F, loc ); P = S->Value( uvAvg.X(), uvAvg.Y() ).Transformed( loc ); @@ -1307,7 +1384,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 +1488,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 +1517,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; } @@ -1518,6 +1610,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 ) @@ -1934,26 +2048,28 @@ SMDS_MeshFace* SMESH_MesherHelper::AddPolygonalFace (const vectorAddPolygonalFaceWithID(nodes, id); else elem = meshDS->AddPolygonalFace(nodes); } - else { - vector newNodes; + else + { + vector newNodes( nodes.size() * 2 ); + newNodes = nodes; for ( int i = 0; i < nodes.size(); ++i ) { const SMDS_MeshNode* n1 = nodes[i]; const SMDS_MeshNode* n2 = nodes[(i+1)%nodes.size()]; const SMDS_MeshNode* n12 = GetMediumNode( n1, n2, force3d, 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 ); @@ -2135,7 +2251,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 +2286,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 +2298,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 ) @@ -2581,7 +2697,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(); @@ -3145,6 +3261,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, @@ -3311,13 +3449,13 @@ namespace { // Structures used by FixQuadraticElements() } 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; } @@ -4824,12 +4962,14 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // mesure 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; @@ -4839,9 +4979,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; @@ -4857,7 +5004,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(); @@ -4866,17 +5013,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 ); @@ -4920,8 +5067,12 @@ 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 ); } @@ -4929,7 +5080,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // 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 = @@ -4985,9 +5136,9 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, 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 ); - faceHlp.CheckNodeUV( TopoDS::Face( faceHlp.GetSubShape() ), pLink->_mediumNode, - uv, /*tol=*/pLink->Move().Modulus(), /*force=*/true, distXYZ ); - p.SetCoord( distXYZ[1], distXYZ[2], distXYZ[3] ); + 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()); @@ -5038,6 +5189,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 ); @@ -5063,7 +5215,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 ); @@ -5072,6 +5224,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 ) @@ -5082,6 +5235,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 ); }