X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=19461a7744fc0c6d7acf98e2da70bdd51eb5f962;hp=4d1321497a6d08eb6ad1fe25b8c4b43227dc3647;hb=5552aec787c289730a35843d295811cfdadaacd3;hpb=b21a1e5b25983280805e04c04b166b9c12bba0af diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index 4d1321497..19461a774 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2016 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 @@ -30,6 +30,7 @@ #include "SMDS_FacePosition.hxx" #include "SMDS_MeshElement.hxx" #include "SMDS_MeshNode.hxx" +#include "SMESHDS_Mesh.hxx" #include "SMESH_Block.hxx" #include "SMESH_Comment.hxx" #include "SMESH_Gen.hxx" @@ -50,7 +51,6 @@ #include #include #include -#include #include #include #include @@ -64,26 +64,25 @@ #include "utilities.h" #include "Utils_ExceptHandlers.hxx" -#ifndef StdMeshers_Array2OfNode_HeaderFile -#define StdMeshers_Array2OfNode_HeaderFile -typedef const SMDS_MeshNode* SMDS_MeshNodePtr; -typedef NCollection_Array2 StdMeshers_Array2OfNode; -#endif +#include +#include -using namespace std; +typedef NCollection_Array2 StdMeshers_Array2OfNode; -typedef gp_XY gp_UV; +typedef gp_XY gp_UV; typedef SMESH_Comment TComm; +using namespace std; + //============================================================================= /*! * */ //============================================================================= -StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, +StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, SMESH_Gen* gen) - : SMESH_2D_Algo(hypId, studyId, gen), + : SMESH_2D_Algo(hypId, gen), myQuadranglePreference(false), myTrianglePreference(false), myTriaVertexID(-1), @@ -93,7 +92,6 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, myQuadType(QUAD_STANDARD), myHelper( NULL ) { - MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D"); _name = "Quadrangle_2D"; _shapeType = (1 << TopAbs_FACE); _compatibleHypothesis.push_back("QuadrangleParams"); @@ -110,7 +108,6 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D() { - MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D"); } //============================================================================= @@ -130,10 +127,10 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis myTrianglePreference = false; myHelper = (SMESH_MesherHelper*)NULL; myParams = NULL; + myProxyMesh.reset(); myQuadList.clear(); - bool isOk = true; - aStatus = SMESH_Hypothesis::HYP_OK; + aStatus = SMESH_Hypothesis::HYP_OK; const list & hyps = GetUsedHypothesis(aMesh, aShape, false); @@ -161,7 +158,7 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis } else if (strcmp("TrianglePreference", aHyp->GetName()) == 0){ isFirstParams = false; - myTrianglePreference = true; + myTrianglePreference = true; } else { isFirstParams = false; @@ -174,18 +171,18 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis if (isFirstParams) { if (strcmp("QuadranglePreference", aHyp->GetName()) == 0) { myQuadranglePreference = true; - myTrianglePreference = false; + myTrianglePreference = false; myQuadType = QUAD_STANDARD; } else if (strcmp("TrianglePreference", aHyp->GetName()) == 0){ myQuadranglePreference = false; - myTrianglePreference = true; + myTrianglePreference = true; myQuadType = QUAD_STANDARD; } } - else { - const StdMeshers_QuadrangleParams* aHyp2 = - (const StdMeshers_QuadrangleParams*)aHyp; + else if (const StdMeshers_QuadrangleParams* aHyp2 = + dynamic_cast( aHyp )) + { myTriaVertexID = aHyp2->GetTriaVertex(); if (!myQuadranglePreference && !myTrianglePreference) { // priority of hypos @@ -206,7 +203,7 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis //============================================================================= /*! - * + * */ //============================================================================= @@ -228,10 +225,11 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, myHelper = &helper; _quadraticMesh = myHelper->IsQuadraticSubMesh(aShape); + myHelper->SetElementsOnShape( true ); myNeedSmooth = false; myCheckOri = false; - FaceQuadStruct::Ptr quad = CheckNbEdges( aMesh, F, /*considerMesh=*/true ); + FaceQuadStruct::Ptr quad = CheckNbEdges( aMesh, F, /*considerMesh=*/true, myHelper ); if (!quad) return false; myQuadList.clear(); @@ -249,16 +247,16 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, enum { NOT_COMPUTED = -1, COMPUTE_FAILED = 0, COMPUTE_OK = 1 }; int res = NOT_COMPUTED; - if (myQuadranglePreference) + if ( myQuadranglePreference ) { int nfull = n1+n2+n3+n4; if ((nfull % 2) == 0 && ((n1 != n3) || (n2 != n4))) { - // special path genarating only quandrangle faces + // special path generating only quandrangle faces res = computeQuadPref( aMesh, F, quad ); } } - else if (myQuadType == QUAD_REDUCED) + else if ( myQuadType == QUAD_REDUCED ) { int n13 = n1 - n3; int n24 = n2 - n4; @@ -471,10 +469,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, b = quad->uv_grid[ j * nbhoriz + i + 1].node; c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node; d = quad->uv_grid[(j + 1) * nbhoriz + i ].node; - SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d); - if (face) { - meshDS->SetMeshElementOnShape(face, geomFaceID); - } + myHelper->AddFace(a, b, c, d); } } @@ -495,7 +490,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, int nbright = (int) uv_e1.size(); int nbleft = (int) uv_e3.size(); - if (quad->nbNodeOut(0) && nbvertic == 2) // this should not occure + if (quad->nbNodeOut(0) && nbvertic == 2) // this should not occur { // Down edge is out // @@ -522,7 +517,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, // for each node of the down edge find nearest node // in the first row of the regular grid and link them for (i = 0; i < stop; i++) { - const SMDS_MeshNode *a, *b, *c, *d; + const SMDS_MeshNode *a, *b, *c=0, *d; a = uv_e0[i].node; b = uv_e0[i + 1].node; gp_Pnt pb (b->X(), b->Y(), b->Z()); @@ -536,6 +531,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, } else { // find in the grid node c, nearest to the b + c = 0; double mind = RealLast(); for (int k = g; k <= iup; k++) { @@ -558,8 +554,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = myHelper->AddFace(a, b, c); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, b, c); } else { // make quadrangle if (near - 1 < ilow) @@ -569,8 +564,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); if (!myTrianglePreference){ - SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, b, c, d); } else { splitQuadFace(meshDS, geomFaceID, a, b, c, d); @@ -584,8 +578,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, d = uv_e3[1].node; else d = quad->uv_grid[nbhoriz + k - 1].node; - SMDS_MeshFace* face = myHelper->AddFace(a, c, d); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, c, d); } } g = near; @@ -609,13 +602,13 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, int g = nbhoriz - 1; // last processed node in the regular grid ilow = 0; - iup = nbhoriz - 1; + iup = nbhoriz - 1; int stop = 0; if ( quad->side[3].grid->Edge(0).IsNull() ) // left side is simulated one { - // quad divided at I but not at J, as nbvertic==nbright==2 - stop++; // we stop at a second node + if ( nbright == 2 ) // quad divided at I but not at J (2D_mesh_QuadranglePreference_01/B1) + stop++; // we stop at a second node } else { @@ -642,20 +635,20 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, d = quad->UVPt( g, nbvertic-2 ).node; if ( myTrianglePreference ) { - if ( SMDS_MeshFace* face = myHelper->AddFace(a, d, c)) - meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, d, c); } else { if ( SMDS_MeshFace* face = myHelper->AddFace(a, b, d, c)) { - meshDS->SetMeshElementOnShape(face, geomFaceID); SMESH_ComputeErrorPtr& err = aMesh.GetSubMesh( aFace )->GetComputeError(); if ( !err || err->IsOK() || err->myName < COMPERR_WARNING ) { - err.reset( new SMESH_ComputeError( COMPERR_WARNING, - "Bad quality quad created")); - err->myBadElements.push_back( face ); + SMESH_BadInputElements* badElems = + new SMESH_BadInputElements( meshDS, COMPERR_WARNING, + "Bad quality quad created"); + badElems->add( face ); + err.reset( badElems ); } } --i; @@ -664,10 +657,11 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, } // for each node of the up edge find nearest node // in the first row of the regular grid and link them - for ( ; i > stop; i--) { + for ( ; i > stop; i--) + { a = uv_e2[i].node; b = uv_e2[i - 1].node; - gp_Pnt pb (b->X(), b->Y(), b->Z()); + gp_Pnt pb = SMESH_TNodeXYZ( b ); // find node c in the grid, which will be linked with node b int near = g; @@ -683,7 +677,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, nk = uv_e1[nbright - 2].node; else nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node; - gp_Pnt pnk (nk->X(), nk->Y(), nk->Z()); + gp_Pnt pnk = SMESH_TNodeXYZ( nk ); double dist = pb.Distance(pnk); if (dist < mind - eps) { c = nk; @@ -696,8 +690,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = myHelper->AddFace(a, b, c); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, b, c); } else { // make quadrangle if (near + 1 > iup) @@ -706,8 +699,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node; //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); if (!myTrianglePreference){ - SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, b, c, d); } else { splitQuadFace(meshDS, geomFaceID, a, b, c, d); @@ -720,8 +712,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, d = uv_e1[nbright - 2].node; else d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node; - SMDS_MeshFace* face = myHelper->AddFace(a, c, d); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, c, d); } } g = near; @@ -731,7 +722,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, } // right or left boundary quadrangles - if (quad->nbNodeOut( QUAD_RIGHT_SIDE ) && nbhoriz == 2) // this should not occure + if (quad->nbNodeOut( QUAD_RIGHT_SIDE ) && nbhoriz == 2) // this should not occur { int g = 0; // last processed node in the grid int stop = nbright - 1; @@ -745,8 +736,9 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, gp_Pnt pb (b->X(), b->Y(), b->Z()); // find node c in the grid, nearest to the b + c = 0; int near = g; - if (i == stop - 1) { // up bondary reached + if (i == stop - 1) { // up boundary reached c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node; near = jup; } else { @@ -770,8 +762,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = myHelper->AddFace(a, b, c); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, b, c); } else { // make quadrangle if (near - 1 < jlow) @@ -781,8 +772,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); if (!myTrianglePreference){ - SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, b, c, d); } else { splitQuadFace(meshDS, geomFaceID, a, b, c, d); @@ -795,16 +785,15 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, d = uv_e0[nbdown - 2].node; else d = quad->uv_grid[nbhoriz*k - 2].node; - SMDS_MeshFace* face = myHelper->AddFace(a, c, d); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, c, d); } } g = near; } } } else { - if (quad->nbNodeOut(3) && nbhoriz == 2) { -// MESSAGE("left edge is out"); + if (quad->nbNodeOut(3) && nbhoriz == 2) + { int g = nbvertic - 1; // last processed node in the grid int stop = 0; i = quad->side[ QUAD_LEFT_SIDE ].to-1; // nbleft - 1; @@ -823,20 +812,20 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, d = quad->UVPt( 1, g ).node; if ( myTrianglePreference ) { - if ( SMDS_MeshFace* face = myHelper->AddFace(a, d, c)) - meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, d, c); } else { if ( SMDS_MeshFace* face = myHelper->AddFace(a, b, d, c)) { - meshDS->SetMeshElementOnShape(face, geomFaceID); SMESH_ComputeErrorPtr& err = aMesh.GetSubMesh( aFace )->GetComputeError(); if ( !err || err->IsOK() || err->myName < COMPERR_WARNING ) { - err.reset( new SMESH_ComputeError( COMPERR_WARNING, - "Bad quality quad created")); - err->myBadElements.push_back( face ); + SMESH_BadInputElements* badElems = + new SMESH_BadInputElements( meshDS, COMPERR_WARNING, + "Bad quality quad created"); + badElems->add( face ); + err.reset( badElems ); } } --i; @@ -851,7 +840,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, // find node c in the grid, nearest to the b int near = g; - if (i == stop + 1) { // down bondary reached + if (i == stop + 1) { // down boundary reached c = quad->uv_grid[nbhoriz*jlow + 1].node; near = jlow; } @@ -876,8 +865,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = myHelper->AddFace(a, b, c); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, b, c); } else { // make quadrangle if (near + 1 > jup) @@ -885,8 +873,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, else d = quad->uv_grid[nbhoriz*(near + 1) + 1].node; if (!myTrianglePreference) { - SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, b, c, d); } else { splitQuadFace(meshDS, geomFaceID, a, b, c, d); @@ -899,8 +886,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, d = quad->uv_grid[nbhoriz*jup + 1].node; //uv_e2[1].node; else d = quad->uv_grid[nbhoriz*(k + 1) + 1].node; - SMDS_MeshFace* face = myHelper->AddFace(a, c, d); - if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); + myHelper->AddFace(a, c, d); } } g = near; @@ -1028,48 +1014,473 @@ bool StdMeshers_Quadrangle_2D::IsApplicable( const TopoDS_Shape & aShape, bool t continue; } - int nbNoDegenEdges = 0; + int nbNoDegenEdges = 0, totalNbEdges = 0; TopExp_Explorer eExp( aFace, TopAbs_EDGE ); - for ( ; eExp.More() && nbNoDegenEdges < 3; eExp.Next() ) { + for ( ; eExp.More() && nbNoDegenEdges < 3; eExp.Next(), ++totalNbEdges ) { if ( !SMESH_Algo::isDegenerated( TopoDS::Edge( eExp.Current() ))) ++nbNoDegenEdges; } - if ( toCheckAll && nbNoDegenEdges < 3 ) return false; - if ( !toCheckAll && nbNoDegenEdges >= 3 ) return true; + if ( toCheckAll && ( totalNbEdges < 4 && nbNoDegenEdges < 3 )) return false; + if ( !toCheckAll && ( totalNbEdges >= 4 || nbNoDegenEdges >= 3 )) return true; } return ( toCheckAll && nbFoundFaces != 0 ); } +namespace +{ + //================================================================================ + /*! + * \brief Return true if only two given edges meat at their common vertex + */ + //================================================================================ + + bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1, + const TopoDS_Edge& e2, + SMESH_Mesh & mesh) + { + TopoDS_Vertex v; + if (!TopExp::CommonVertex(e1, e2, v)) + return false; + TopTools_ListIteratorOfListOfShape ancestIt(mesh.GetAncestors(v)); + for (; ancestIt.More() ; ancestIt.Next()) + if (ancestIt.Value().ShapeType() == TopAbs_EDGE) + if (!e1.IsSame(ancestIt.Value()) && !e2.IsSame(ancestIt.Value())) + return false; + return true; + } + + //-------------------------------------------------------------------------------- + /*! + * \brief EDGE of a FACE + */ + struct Edge + { + TopoDS_Edge myEdge; + TopoDS_Vertex my1stVertex; + int myIndex; + double myAngle; // angle at my1stVertex + int myNbSegments; // discretization + Edge* myPrev; // preceding EDGE + Edge* myNext; // next EDGE + + // traits used by boost::intrusive::circular_list_algorithms + typedef Edge node; + typedef Edge * node_ptr; + typedef const Edge * const_node_ptr; + static node_ptr get_next(const_node_ptr n) { return n->myNext; } + static void set_next(node_ptr n, node_ptr next) { n->myNext = next; } + static node_ptr get_previous(const_node_ptr n) { return n->myPrev; } + static void set_previous(node_ptr n, node_ptr prev){ n->myPrev = prev; } + }; + + //-------------------------------------------------------------------------------- + /*! + * \brief Four sides of a quadrangle evaluating its quality + */ + struct QuadQuality + { + typedef std::set< QuadQuality, QuadQuality > set; + + Edge* myCornerE[4]; + int myNbSeg [4]; + + // quality criteria to minimize + int myOppDiff; + double myQuartDiff; + double mySumAngle; + + // Compute quality criateria and add self to the set of variants + // + void AddSelf( QuadQuality::set& theVariants ) + { + if ( myCornerE[2] == myCornerE[1] || // exclude invalid variants + myCornerE[2] == myCornerE[3] || + myCornerE[0] == myCornerE[3] ) + return; + + // count nb segments between corners + mySumAngle = 0; + double totNbSeg = 0; + for ( int i1 = 3, i2 = 0; i2 < 4; i1 = i2++ ) + { + myNbSeg[ i1 ] = 0; + for ( Edge* e = myCornerE[ i1 ]; e != myCornerE[ i2 ]; e = e->myNext ) + myNbSeg[ i1 ] += e->myNbSegments; + mySumAngle -= myCornerE[ i1 ]->myAngle / M_PI; // [-1,1] + totNbSeg += myNbSeg[ i1 ]; + } + + myOppDiff = ( Abs( myNbSeg[0] - myNbSeg[2] ) + + Abs( myNbSeg[1] - myNbSeg[3] )); + + double nbSideIdeal = totNbSeg / 4.; + myQuartDiff = -( Min( Min( myNbSeg[0], myNbSeg[1] ), + Min( myNbSeg[2], myNbSeg[3] )) / nbSideIdeal ); + + theVariants.insert( *this ); + +#ifndef _DEBUG_ + if ( theVariants.size() > 1 ) // erase a worse variant + theVariants.erase( ++theVariants.begin() ); +#endif + }; + + // first criterion - equality of nbSeg of opposite sides + int crit1() const { return myOppDiff; } + + // second criterion - equality of nbSeg of adjacent sides and sharpness of angles + double crit2() const { return myQuartDiff + mySumAngle; } + + bool operator () ( const QuadQuality& q1, const QuadQuality& q2) const + { + if ( q1.crit1() < q2.crit1() ) + return true; + if ( q1.crit1() > q2.crit1() ) + return false; + return q1.crit2() < q2.crit2(); + } + }; + + //================================================================================ + /*! + * \brief Unite EDGEs to get a required number of sides + * \param [in] theNbCorners - the required number of sides + * \param [in] theConsiderMesh - to considered only meshed VERTEXes + * \param [in] theFaceSide - the FACE EDGEs + * \param [out] theVertices - the found corner vertices + */ + //================================================================================ + + void uniteEdges( const int theNbCorners, + const bool theConsiderMesh, + const StdMeshers_FaceSide& theFaceSide, + const TopoDS_Shape& theBaseVertex, + std::vector& theVertices, + bool& theHaveConcaveVertices) + { + // form a circular list of EDGEs + std::vector< Edge > edges( theFaceSide.NbEdges() ); + boost::intrusive::circular_list_algorithms< Edge > circularList; + circularList.init_header( &edges[0] ); + edges[0].myEdge = theFaceSide.Edge( 0 ); + edges[0].myIndex = 0; + edges[0].myNbSegments = 0; + for ( int i = 1; i < theFaceSide.NbEdges(); ++i ) + { + edges[ i ].myEdge = theFaceSide.Edge( i ); + edges[ i ].myIndex = i; + edges[ i ].myNbSegments = 0; + circularList.link_after( &edges[ i-1 ], &edges[ i ] ); + } + // remove degenerated edges + int nbEdges = edges.size(); + Edge* edge0 = &edges[0]; + for ( size_t i = 0; i < edges.size(); ++i ) + if ( SMESH_Algo::isDegenerated( edges[i].myEdge )) + { + edge0 = circularList.unlink( &edges[i] ); + --nbEdges; + } + + // sort edges by angle + std::multimap< double, Edge* > edgeByAngle; + int i, iBase = -1, nbConvexAngles = 0, nbSharpAngles = 0; + const double angTol = 5. / 180 * M_PI; + const double sharpAngle = 0.5 * M_PI - angTol; + Edge* e = edge0; + for ( i = 0; i < nbEdges; ++i, e = e->myNext ) + { + e->my1stVertex = SMESH_MesherHelper::IthVertex( 0, e->myEdge ); + if ( e->my1stVertex.IsSame( theBaseVertex )) + iBase = e->myIndex; + + e->myAngle = -2 * M_PI; + if ( !theConsiderMesh || theFaceSide.VertexNode( e->myIndex )) + { + e->myAngle = SMESH_MesherHelper::GetAngle( e->myPrev->myEdge, e->myEdge, + theFaceSide.Face(), e->my1stVertex ); + if ( e->myAngle > 2 * M_PI ) // GetAngle() failed + e->myAngle *= -1.; + } + edgeByAngle.insert( std::make_pair( e->myAngle, e )); + nbConvexAngles += ( e->myAngle > angTol ); + nbSharpAngles += ( e->myAngle > sharpAngle ); + } + + theHaveConcaveVertices = ( nbConvexAngles < nbEdges ); + + if ((int) theVertices.size() == theNbCorners ) + return; + + theVertices.clear(); + + if ( !theConsiderMesh || theNbCorners < 4 || + nbConvexAngles <= theNbCorners || + nbSharpAngles == theNbCorners ) + { + if ( nbEdges == theNbCorners ) // return all vertices + { + for ( e = edge0; (int) theVertices.size() < theNbCorners; e = e->myNext ) + theVertices.push_back( e->my1stVertex ); + return; + } + + // return corners with maximal angles + + std::set< int > cornerIndices; + if ( iBase != -1 ) + cornerIndices.insert( iBase ); + + std::multimap< double, Edge* >::reverse_iterator a2e = edgeByAngle.rbegin(); + for (; (int) cornerIndices.size() < theNbCorners; ++a2e ) + cornerIndices.insert( a2e->second->myIndex ); + + std::set< int >::iterator i = cornerIndices.begin(); + for ( ; i != cornerIndices.end(); ++i ) + theVertices.push_back( edges[ *i ].my1stVertex ); + + return; + } + + // get nb of segments + int totNbSeg = 0; // tatal nb segments + std::vector nodes; + for ( i = 0, e = edge0; i < nbEdges; ++i, e = e->myNext ) + { + nodes.clear(); + theFaceSide.GetEdgeNodes( e->myIndex, nodes, /*addVertex=*/true, true ); + if ( nodes.size() == 2 && nodes[0] == nodes[1] ) // all nodes merged + { + e->myAngle = -1; // to remove + } + else + { + e->myNbSegments += nodes.size() - 1; + totNbSeg += nodes.size() - 1; + } + + // join with the previous edge those edges with concave angles + if ( e->myAngle <= 0 ) + { + e->myPrev->myNbSegments += e->myNbSegments; + e = circularList.unlink( e )->myPrev; + --nbEdges; + --i; + } + } + + if ( edge0->myNext->myPrev != edge0 ) // edge0 removed, find another edge0 + for ( size_t i = 0; i < edges.size(); ++i ) + if ( edges[i].myNext->myPrev == & edges[i] ) + { + edge0 = &edges[i]; + break; + } + + + // sort different variants by quality + + QuadQuality::set quadVariants; + + // find index of a corner most opposite to corner of edge0 + int iOpposite0, nbHalf = 0; + for ( e = edge0; nbHalf <= totNbSeg / 2; e = e->myNext ) + nbHalf += e->myNbSegments; + iOpposite0 = e->myIndex; + + // compose different variants of quadrangles + QuadQuality quad; + for ( ; edge0->myIndex != iOpposite0; edge0 = edge0->myNext ) + { + quad.myCornerE[ 0 ] = edge0; + + // find opposite corner 2 + for ( nbHalf = 0, e = edge0; nbHalf < totNbSeg / 2; e = e->myNext ) + nbHalf += e->myNbSegments; + if ( e == edge0->myNext ) // no space for corner 1 + e = e->myNext; + quad.myCornerE[ 2 ] = e; + + bool moreVariants2 = ( totNbSeg % 2 || nbHalf != totNbSeg / 2 ); + + // enumerate different variants of corners 1 and 3 + for ( Edge* e1 = edge0->myNext; e1 != quad.myCornerE[ 2 ]; e1 = e1->myNext ) + { + quad.myCornerE[ 1 ] = e1; + + // find opposite corner 3 + for ( nbHalf = 0, e = e1; nbHalf < totNbSeg / 2; e = e->myNext ) + nbHalf += e->myNbSegments; + if ( e == quad.myCornerE[ 2 ] ) + e = e->myNext; + quad.myCornerE[ 3 ] = e; + + bool moreVariants3 = ( totNbSeg % 2 || nbHalf != totNbSeg / 2 ); + + quad.AddSelf( quadVariants ); + + // another variants + if ( moreVariants2 ) + { + quad.myCornerE[ 2 ] = quad.myCornerE[ 2 ]->myPrev; + quad.AddSelf( quadVariants ); + quad.myCornerE[ 2 ] = quad.myCornerE[ 2 ]->myNext; + } + if ( moreVariants3 ) + { + quad.myCornerE[ 3 ] = quad.myCornerE[ 3 ]->myPrev; + quad.AddSelf( quadVariants ); + + if ( moreVariants2 ) + { + quad.myCornerE[ 2 ] = quad.myCornerE[ 2 ]->myPrev; + quad.AddSelf( quadVariants ); + quad.myCornerE[ 2 ] = quad.myCornerE[ 2 ]->myNext; + } + } + } + } + + const QuadQuality& bestQuad = *quadVariants.begin(); + theVertices.resize( 4 ); + theVertices[ 0 ] = bestQuad.myCornerE[ 0 ]->my1stVertex; + theVertices[ 1 ] = bestQuad.myCornerE[ 1 ]->my1stVertex; + theVertices[ 2 ] = bestQuad.myCornerE[ 2 ]->my1stVertex; + theVertices[ 3 ] = bestQuad.myCornerE[ 3 ]->my1stVertex; + + return; + } + +} // namespace + //================================================================================ /*! - * \brief Return true if only two given edges meat at their common vertex + * \brief Finds vertices at the most sharp face corners + * \param [in] theFace - the FACE + * \param [in,out] theWire - the ordered edges of the face. It can be modified to + * have the first VERTEX of the first EDGE in \a vertices + * \param [out] theVertices - the found corner vertices in the order corresponding to + * the order of EDGEs in \a theWire + * \param [out] theNbDegenEdges - nb of degenerated EDGEs in theFace + * \param [in] theConsiderMesh - if \c true, only meshed VERTEXes are considered + * as possible corners + * \return int - number of quad sides found: 0, 3 or 4 */ //================================================================================ -static bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1, - const TopoDS_Edge& e2, - SMESH_Mesh & mesh) +int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, + SMESH_Mesh & theMesh, + std::list& theWire, + std::vector& theVertices, + int & theNbDegenEdges, + const bool theConsiderMesh) { - TopoDS_Vertex v; - if (!TopExp::CommonVertex(e1, e2, v)) - return false; - TopTools_ListIteratorOfListOfShape ancestIt(mesh.GetAncestors(v)); - for (; ancestIt.More() ; ancestIt.Next()) - if (ancestIt.Value().ShapeType() == TopAbs_EDGE) - if (!e1.IsSame(ancestIt.Value()) && !e2.IsSame(ancestIt.Value())) - return false; - return true; + theNbDegenEdges = 0; + + SMESH_MesherHelper helper( theMesh ); + if ( myHelper ) + helper.CopySubShapeInfo( *myHelper ); + + StdMeshers_FaceSide faceSide( theFace, theWire, &theMesh, + /*isFwd=*/true, /*skipMedium=*/true, &helper ); + + // count degenerated EDGEs and possible corner VERTEXes + for ( int iE = 0; iE < faceSide.NbEdges(); ++iE ) + { + if ( SMESH_Algo::isDegenerated( faceSide.Edge( iE ))) + ++theNbDegenEdges; + else if ( !theConsiderMesh || faceSide.VertexNode( iE )) + theVertices.push_back( faceSide.FirstVertex( iE )); + } + + // find out required nb of corners (3 or 4) + int nbCorners = 4; + TopoDS_Shape triaVertex = helper.GetMeshDS()->IndexToShape( myTriaVertexID ); + if ( !triaVertex.IsNull() && + triaVertex.ShapeType() == TopAbs_VERTEX && + helper.IsSubShape( triaVertex, theFace ) && + theVertices.size() != 4 ) + nbCorners = 3; + else + triaVertex.Nullify(); + + // check nb of available EDGEs + if ( faceSide.NbEdges() < nbCorners ) + return error(COMPERR_BAD_SHAPE, + TComm("Face must have 4 sides but not ") << faceSide.NbEdges() ); + + if ( theConsiderMesh ) + { + const int nbSegments = Max( faceSide.NbPoints()-1, faceSide.NbSegments() ); + if ( nbSegments < nbCorners ) + return error(COMPERR_BAD_INPUT_MESH, TComm("Too few boundary nodes: ") << nbSegments); + } + + if ( nbCorners == 3 ) + { + if ( theVertices.size() < 3 ) + return error(COMPERR_BAD_SHAPE, + TComm("Face must have 3 meshed sides but not ") << theVertices.size() ); + } + else // triaVertex not defined or invalid + { + if ( theVertices.size() == 3 && theNbDegenEdges == 0 ) + { + if ( myTriaVertexID < 1 ) + return error(COMPERR_BAD_PARMETERS, + "No Base vertex provided for a trilateral geometrical face"); + + TComm comment("Invalid Base vertex: "); + comment << myTriaVertexID << ", which is not in [ "; + comment << helper.GetMeshDS()->ShapeToIndex( faceSide.FirstVertex(0) ) << ", "; + comment << helper.GetMeshDS()->ShapeToIndex( faceSide.FirstVertex(1) ) << ", "; + comment << helper.GetMeshDS()->ShapeToIndex( faceSide.FirstVertex(2) ) << " ]"; + return error(COMPERR_BAD_PARMETERS, comment ); + } + if ( theVertices.size() + theNbDegenEdges < 4 ) + return error(COMPERR_BAD_SHAPE, + TComm("Face must have 4 meshed sides but not ") << theVertices.size() ); + } + + myCheckOri = false; + if ( theVertices.size() > 3 ) + { + uniteEdges( nbCorners, theConsiderMesh, faceSide, triaVertex, theVertices, myCheckOri ); + } + + if ( nbCorners == 3 && !triaVertex.IsSame( theVertices[0] )) + { + // make theVertices begin from triaVertex + for ( size_t i = 0; i < theVertices.size(); ++i ) + if ( triaVertex.IsSame( theVertices[i] )) + { + theVertices.erase( theVertices.begin(), theVertices.begin() + i ); + break; + } + else + { + theVertices.push_back( theVertices[i] ); + } + } + + // make theWire begin from the 1st corner vertex + while ( !theVertices[0].IsSame( helper.IthVertex( 0, theWire.front() )) || + SMESH_Algo::isDegenerated( theWire.front() )) + theWire.splice( theWire.end(), theWire, theWire.begin() ); + + return nbCorners; } //============================================================================= /*! - * + * */ //============================================================================= FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape, - const bool considerMesh) + const bool considerMesh, + SMESH_MesherHelper* aFaceHelper) { if ( !myQuadList.empty() && myQuadList.front()->face.IsSame( aShape )) return myQuadList.front(); @@ -1088,6 +1499,7 @@ FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & } // find corner vertices of the quad + myHelper = ( aFaceHelper && aFaceHelper->GetSubShape() == aShape ) ? aFaceHelper : NULL; vector corners; int nbDegenEdges, nbSides = getCorners( F, aMesh, edges, corners, nbDegenEdges, considerMesh ); if ( nbSides == 0 ) @@ -1113,7 +1525,7 @@ FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & sideEdges.push_back( *edgeIt++ ); if ( !sideEdges.empty() ) quad->side.push_back( StdMeshers_FaceSide::New(F, sideEdges, &aMesh, iSide < QUAD_TOP_SIDE, - ignoreMediumNodes, myProxyMesh)); + ignoreMediumNodes, myHelper, myProxyMesh)); else --iSide; } @@ -1155,7 +1567,7 @@ FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & } } } - else + else //if ( !myHelper || !myHelper->IsRealSeam( edge )) { sideEdges.push_back( edge ); } @@ -1167,7 +1579,7 @@ FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & { quad->side.push_back ( StdMeshers_FaceSide::New( F, sideEdges, &aMesh, iSide < QUAD_TOP_SIDE, - ignoreMediumNodes, myProxyMesh )); + ignoreMediumNodes, myHelper, myProxyMesh )); ++iSide; } if ( quad->side.size() == 4 ) @@ -1453,6 +1865,8 @@ bool StdMeshers_Quadrangle_2D::setNormalizedGrid (FaceQuadStruct::Ptr quad) int nbhoriz = Min( bSide.NbPoints(), tSide.NbPoints() ); int nbvertic = Min( rSide.NbPoints(), lSide.NbPoints() ); + if ( nbhoriz < 1 || nbvertic < 1 ) + return error("Algo error: empty quad"); if ( myQuadList.size() == 1 ) { @@ -1589,7 +2003,7 @@ bool StdMeshers_Quadrangle_2D::setNormalizedGrid (FaceQuadStruct::Ptr quad) //======================================================================= //function : ShiftQuad -//purpose : auxilary function for computeQuadPref +//purpose : auxiliary function for computeQuadPref //======================================================================= void StdMeshers_Quadrangle_2D::shiftQuad(FaceQuadStruct::Ptr& quad, const int num ) @@ -1599,7 +2013,7 @@ void StdMeshers_Quadrangle_2D::shiftQuad(FaceQuadStruct::Ptr& quad, const int nu //================================================================================ /*! - * \brief Rotate sides of a quad by given nb of quartes + * \brief Rotate sides of a quad CCW by given nb of quartes * \param nb - number of rotation quartes * \param ori - to keep orientation of sides as in an unit quad or not * \param keepGrid - if \c true Side::grid is not changed, Side::from and Side::to @@ -1611,6 +2025,8 @@ void FaceQuadStruct::shift( size_t nb, bool ori, bool keepGrid ) { if ( nb == 0 ) return; + nb = nb % NB_QUAD_SIDES; + vector< Side > newSides( side.size() ); vector< Side* > sidePtrs( side.size() ); for (int i = QUAD_BOTTOM_SIDE; i < NB_QUAD_SIDES; ++i) @@ -1640,12 +2056,38 @@ void FaceQuadStruct::shift( size_t nb, bool ori, bool keepGrid ) } newSides.swap( side ); - uv_grid.clear(); + if ( keepGrid && !uv_grid.empty() ) + { + if ( nb == 2 ) // "PI" + { + std::reverse( uv_grid.begin(), uv_grid.end() ); + } + else + { + FaceQuadStruct newQuad; + newQuad.uv_grid.resize( uv_grid.size() ); + newQuad.iSize = jSize; + newQuad.jSize = iSize; + int i, j, iRev, jRev; + int *iNew = ( nb == 1 ) ? &jRev : &j; + int *jNew = ( nb == 1 ) ? &i : &iRev; + for ( i = 0, iRev = iSize-1; i < iSize; ++i, --iRev ) + for ( j = 0, jRev = jSize-1; j < jSize; ++j, --jRev ) + newQuad.UVPt( *iNew, *jNew ) = UVPt( i, j ); + + std::swap( iSize, jSize ); + std::swap( uv_grid, newQuad.uv_grid ); + } + } + else + { + uv_grid.clear(); + } } //======================================================================= //function : calcUV -//purpose : auxilary function for computeQuadPref +//purpose : auxiliary function for computeQuadPref //======================================================================= static gp_UV calcUV(double x0, double x1, double y0, double y1, @@ -1668,7 +2110,7 @@ static gp_UV calcUV(double x0, double x1, double y0, double y1, //======================================================================= //function : calcUV2 -//purpose : auxilary function for computeQuadPref +//purpose : auxiliary function for computeQuadPref //======================================================================= static gp_UV calcUV2(double x, double y, @@ -1749,7 +2191,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, // | | | | // | |C | | // | L | | R | - // left | |__| | rigth + // left | |__| | right // | / \ | // | / C \ | // |/ \| @@ -1762,7 +2204,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, // | |__| | // | / \ | // | / C \ | - // left |/________\| rigth + // left |/________\| right // | | // | C | // | | @@ -1770,10 +2212,10 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, // 0 bottom 1 - const int bfrom = quad->side[0].from; - const int rfrom = quad->side[1].from; + //const int bfrom = quad->side[0].from; + //const int rfrom = quad->side[1].from; const int tfrom = quad->side[2].from; - const int lfrom = quad->side[3].from; + //const int lfrom = quad->side[3].from; { const vector& uv_eb_vec = quad->side[0].GetUVPtStruct(true,0); const vector& uv_er_vec = quad->side[1].GetUVPtStruct(false,1); @@ -1825,6 +2267,10 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, } sideLCb = StdMeshers_FaceSide::New( pointsLCb, aFace ); p3dom = pointsLCb.back(); + + gp_Pnt xyz = S->Value( p3dom.u, p3dom.v ); + p3dom.node = myHelper->AddNode( xyz.X(), xyz.Y(), xyz.Z(), 0, p3dom.u, p3dom.v ); + pointsLCb.back() = p3dom; } // Make a side separating domains L and Ct StdMeshers_FaceSidePtr sideLCt; @@ -1910,10 +2356,10 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, } // if ( dv != 0 && dh != 0 ) - const int db = quad->side[0].IsReversed() ? -1 : +1; - const int dr = quad->side[1].IsReversed() ? -1 : +1; + //const int db = quad->side[0].IsReversed() ? -1 : +1; + //const int dr = quad->side[1].IsReversed() ? -1 : +1; const int dt = quad->side[2].IsReversed() ? -1 : +1; - const int dl = quad->side[3].IsReversed() ? -1 : +1; + //const int dl = quad->side[3].IsReversed() ? -1 : +1; // Case dv == 0, here possibly myQuadList.size() > 1 // @@ -1976,6 +2422,16 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, sideRCb = StdMeshers_FaceSide::New( pointsRCb, aFace ); pTBL = pointsLCb.back(); pTBR = pointsRCb.back(); + { + gp_Pnt xyz = S->Value( pTBL.u, pTBL.v ); + pTBL.node = myHelper->AddNode( xyz.X(), xyz.Y(), xyz.Z(), 0, pTBL.u, pTBL.v ); + pointsLCb.back() = pTBL; + } + { + gp_Pnt xyz = S->Value( pTBR.u, pTBR.v ); + pTBR.node = myHelper->AddNode( xyz.X(), xyz.Y(), xyz.Z(), 0, pTBR.u, pTBR.v ); + pointsRCb.back() = pTBR; + } } // Make sides separating domains Ct and L and R StdMeshers_FaceSidePtr sideLCt, sideRCt; @@ -2113,7 +2569,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, npl.Append(uv_el[i].normParam); } - int dl,dr; + int dl = 0, dr = 0; if (OldVersion) { // add some params to right and left after the first param // insert to right @@ -2131,7 +2587,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, } int nnn = Min(nr,nl); - // auxilary sequence of XY for creation nodes + // auxiliary sequence of XY for creation nodes // in the bottom part of central domain // Length of UVL and UVR must be == nbv-nnn TColgp_SequenceOfXY UVL, UVR, UVT; @@ -2179,10 +2635,8 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, for (i=1; i<=dl; i++) { for (j=1; jAddFace(NodesL.Value(i,j), NodesL.Value(i+1,j), - NodesL.Value(i+1,j+1), NodesL.Value(i,j+1)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); + myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j), + NodesL.Value(i+1,j+1), NodesL.Value(i,j+1)); } } } @@ -2236,10 +2690,8 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, for (i=1; i<=dr; i++) { for (j=1; jAddFace(NodesR.Value(i,j), NodesR.Value(i+1,j), - NodesR.Value(i+1,j+1), NodesR.Value(i,j+1)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); + myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j), + NodesR.Value(i+1,j+1), NodesR.Value(i,j+1)); } } } @@ -2309,10 +2761,8 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, for (i=1; iAddFace(NodesC.Value(i,j), NodesC.Value(i+1,j), - NodesC.Value(i+1,j+1), NodesC.Value(i,j+1)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); + myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j), + NodesC.Value(i+1,j+1), NodesC.Value(i,j+1)); } } } @@ -2341,10 +2791,8 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, for (j=1; jAddFace(NodesBRD.Value(i,j), NodesBRD.Value(i+1,j), - NodesBRD.Value(i+1,j+1), NodesBRD.Value(i,j+1)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); + myHelper->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i+1,j), + NodesBRD.Value(i+1,j+1), NodesBRD.Value(i,j+1)); } } } @@ -2363,7 +2811,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, TColgp_SequenceOfXY UVtmp; double drparam = npr.Value(nr) - npr.Value(nnn-1); double dlparam = npl.Value(nnn) - npl.Value(nnn-1); - double y0,y1; + double y0 = 0, y1 = 0; for (i=1; i<=drl; i++) { // add existed nodes from right edge NodesC.SetValue(nb,i+1,uv_er[nnn+i-2].node); @@ -2447,10 +2895,8 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, for (j=1; j<=drl+addv; j++) { for (i=1; iAddFace(NodesC.Value(i,j), NodesC.Value(i+1,j), - NodesC.Value(i+1,j+1), NodesC.Value(i,j+1)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); + myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j), + NodesC.Value(i+1,j+1), NodesC.Value(i,j+1)); } } } // end nrAddFace(NodesLast.Value(i,1), NodesLast.Value(i+1,1), - NodesLast.Value(i+1,2), NodesLast.Value(i,2)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); + myHelper->AddFace(NodesLast.Value(i,1), NodesLast.Value(i+1,1), + NodesLast.Value(i+1,2), NodesLast.Value(i,2)); } } } // if ((drl+addv) > 0) @@ -2501,7 +2945,7 @@ bool StdMeshers_Quadrangle_2D::evaluateQuadPref(SMESH_Mesh & aMesh, MapShapeNbElems& aResMap, bool IsQuadratic) { - // Auxilary key in order to keep old variant + // Auxiliary key in order to keep old variant // of meshing after implementation new variant // for bug 0016220 from Mantis. bool OldVersion = false; @@ -2636,21 +3080,16 @@ void StdMeshers_Quadrangle_2D::splitQuadFace(SMESHDS_Mesh * theMeshDS, const SMDS_MeshNode* theNode3, const SMDS_MeshNode* theNode4) { - SMDS_MeshFace* face; if ( SMESH_TNodeXYZ( theNode1 ).SquareDistance( theNode3 ) > SMESH_TNodeXYZ( theNode2 ).SquareDistance( theNode4 ) ) { - face = myHelper->AddFace(theNode2, theNode4 , theNode1); - if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID); - face = myHelper->AddFace(theNode2, theNode3, theNode4); - if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID); + myHelper->AddFace(theNode2, theNode4 , theNode1); + myHelper->AddFace(theNode2, theNode3, theNode4); } else { - face = myHelper->AddFace(theNode1, theNode2 ,theNode3); - if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID); - face = myHelper->AddFace(theNode1, theNode3, theNode4); - if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID); + myHelper->AddFace(theNode1, theNode2 ,theNode3); + myHelper->AddFace(theNode1, theNode3, theNode4); } } @@ -2962,7 +3401,8 @@ bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, const vector& uv_et = quad->side[2].GetUVPtStruct(true,1); const vector& uv_el = quad->side[3].GetUVPtStruct(false,0); - if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl) + if ((int) uv_eb.size() != nb || (int) uv_er.size() != nr || + (int) uv_et.size() != nt || (int) uv_el.size() != nl) return error(COMPERR_BAD_INPUT_MESH); // arrays for normalized params @@ -2987,7 +3427,7 @@ bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, // | | | | // | | | | // | L | | R | - // left | | | | rigth + // left | | | | right // | / \ | // | / C \ | // |/ \| @@ -3014,7 +3454,7 @@ bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, gp_XY a3 (uv_et.front().u, uv_et.front().v); int nnn = Min(nr,nl); - // auxilary sequence of XY for creation of nodes + // auxiliary sequence of XY for creation of nodes // in the bottom part of central domain // it's length must be == nbv-nnn-1 TColgp_SequenceOfXY UVL; @@ -3062,10 +3502,8 @@ bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, // create faces for (i=1; i<=dl; i++) { for (j=1; jAddFace(NodesL.Value(i,j), NodesL.Value(i+1,j), - NodesL.Value(i+1,j+1), NodesL.Value(i,j+1)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); + myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j), + NodesL.Value(i+1,j+1), NodesL.Value(i,j+1)); } } } @@ -3117,10 +3555,8 @@ bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, // create faces for (i=1; i<=dr; i++) { for (j=1; jAddFace(NodesR.Value(i,j), NodesR.Value(i+1,j), - NodesR.Value(i+1,j+1), NodesR.Value(i,j+1)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); + myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j), + NodesR.Value(i+1,j+1), NodesR.Value(i,j+1)); } } } @@ -3181,10 +3617,8 @@ bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, // create faces for (i=1; iAddFace(NodesC.Value(i,j), NodesC.Value(i+1,j), - NodesC.Value(i+1,j+1), NodesC.Value(i,j+1)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); + myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j), + NodesC.Value(i+1,j+1), NodesC.Value(i,j+1)); } } } // end Multiple Reduce implementation @@ -3281,11 +3715,10 @@ bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, const vector& uv_et = quad->side[2].GetUVPtStruct(true,1); const vector& uv_el = quad->side[3].GetUVPtStruct(false,0); - if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl) + if ((int) uv_eb.size() != nb || (int) uv_er.size() != nr || + (int) uv_et.size() != nt || (int) uv_el.size() != nl) return error(COMPERR_BAD_INPUT_MESH); - myHelper->SetElementsOnShape( true ); - gp_UV uv[ UV_SIZE ]; uv[ UV_A0 ].SetCoord( uv_eb.front().u, uv_eb.front().v); uv[ UV_A1 ].SetCoord( uv_eb.back().u, uv_eb.back().v ); @@ -3294,7 +3727,11 @@ bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, vector curr_base = uv_eb, next_base; - UVPtStruct nullUVPtStruct; nullUVPtStruct.node = 0; + UVPtStruct nullUVPtStruct; + nullUVPtStruct.node = 0; + nullUVPtStruct.x = nullUVPtStruct.y = nullUVPtStruct.u = nullUVPtStruct.v = 0; + nullUVPtStruct.param = 0; + int curr_base_len = nb; int next_base_len = 0; @@ -3847,7 +4284,7 @@ void StdMeshers_Quadrangle_2D::updateDegenUV(FaceQuadStruct::Ptr quad) // Set number of nodes on a degenerated side to be same as on an opposite side // ---------------------------------------------------------------------------- - for ( unsigned i = 0; i < quad->side.size(); ++i ) + for ( size_t i = 0; i < quad->side.size(); ++i ) { StdMeshers_FaceSidePtr degSide = quad->side[i]; if ( !myHelper->IsDegenShape( degSide->EdgeID(0) )) @@ -3878,136 +4315,185 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad) { if ( !myNeedSmooth ) return; - // Get nodes to smooth + SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); + const double tol = BRep_Tool::Tolerance( quad->face ); + Handle(ShapeAnalysis_Surface) surface = myHelper->GetSurface( quad->face ); - // TODO: do not smooth fixed nodes + if ( myHelper->HasDegeneratedEdges() && myForcedPnts.empty() ) + { + // "smooth" by computing node positions using 3D TFI and further projection - typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap; - TNo2SmooNoMap smooNoMap; + list< FaceQuadStruct::Ptr >::iterator q = myQuadList.begin(); + for ( ; q != myQuadList.end() ; ++q ) + { + quad = *q; + int nbhoriz = quad->iSize; + int nbvertic = quad->jSize; - const TopoDS_Face& geomFace = TopoDS::Face( myHelper->GetSubShape() ); - Handle(Geom_Surface) surface = BRep_Tool::Surface( geomFace ); - double U1, U2, V1, V2; - surface->Bounds(U1, U2, V1, V2); - GeomAPI_ProjectPointOnSurf proj; - proj.Init( surface, U1, U2, V1, V2, BRep_Tool::Tolerance( geomFace ) ); + SMESH_TNodeXYZ a0( quad->UVPt( 0, 0 ).node ); + SMESH_TNodeXYZ a1( quad->UVPt( nbhoriz-1, 0 ).node ); + SMESH_TNodeXYZ a2( quad->UVPt( nbhoriz-1, nbvertic-1 ).node ); + SMESH_TNodeXYZ a3( quad->UVPt( 0, nbvertic-1 ).node ); - SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); - SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace ); - SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes(); - while ( nIt->more() ) // loop on nodes bound to a FACE - { - const SMDS_MeshNode* node = nIt->next(); - TSmoothNode & sNode = smooNoMap[ node ]; - sNode._uv = myHelper->GetNodeUV( geomFace, node ); - sNode._xyz = SMESH_TNodeXYZ( node ); - - // set sNode._triangles - SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face ); - while ( fIt->more() ) - { - const SMDS_MeshElement* face = fIt->next(); - const int nbN = face->NbCornerNodes(); - const int nInd = face->GetNodeIndex( node ); - const int prevInd = myHelper->WrapIndex( nInd - 1, nbN ); - const int nextInd = myHelper->WrapIndex( nInd + 1, nbN ); - const SMDS_MeshNode* prevNode = face->GetNode( prevInd ); - const SMDS_MeshNode* nextNode = face->GetNode( nextInd ); - sNode._triangles.push_back( TTriangle( & smooNoMap[ prevNode ], - & smooNoMap[ nextNode ])); + for (int i = 1; i < nbhoriz-1; i++) + { + SMESH_TNodeXYZ p0( quad->UVPt( i, 0 ).node ); + SMESH_TNodeXYZ p2( quad->UVPt( i, nbvertic-1 ).node ); + for (int j = 1; j < nbvertic-1; j++) + { + SMESH_TNodeXYZ p1( quad->UVPt( nbhoriz-1, j ).node ); + SMESH_TNodeXYZ p3( quad->UVPt( 0, j ).node ); + + UVPtStruct& uvp = quad->UVPt( i, j ); + + gp_Pnt p = myHelper->calcTFI(uvp.x,uvp.y, a0,a1,a2,a3, p0,p1,p2,p3); + gp_Pnt2d uv = surface->NextValueOfUV( uvp.UV(), p, 10*tol ); + gp_Pnt pnew = surface->Value( uv ); + + meshDS->MoveNode( uvp.node, pnew.X(), pnew.Y(), pnew.Z() ); + uvp.u = uv.X(); + uvp.v = uv.Y(); + } + } } } - // set _uv of smooth nodes on FACE boundary - for ( unsigned i = 0; i < quad->side.size(); ++i ) + else { - const vector& uvVec = quad->side[i].GetUVPtStruct(); - for ( unsigned j = 0; j < uvVec.size(); ++j ) + // Get nodes to smooth + + typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap; + TNo2SmooNoMap smooNoMap; + + // fixed nodes + boost::container::flat_set< const SMDS_MeshNode* > fixedNodes; + for ( size_t i = 0; i < myForcedPnts.size(); ++i ) { - TSmoothNode & sNode = smooNoMap[ uvVec[j].node ]; - sNode._uv = uvVec[j].UV(); - sNode._xyz = SMESH_TNodeXYZ( uvVec[j].node ); + fixedNodes.insert( myForcedPnts[i].node ); + if ( myForcedPnts[i].node->getshapeId() != myHelper->GetSubShapeID() ) + { + TSmoothNode & sNode = smooNoMap[ myForcedPnts[i].node ]; + sNode._uv = myForcedPnts[i].uv; + sNode._xyz = SMESH_TNodeXYZ( myForcedPnts[i].node ); + } } - } + SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( quad->face ); + SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes(); + while ( nIt->more() ) // loop on nodes bound to a FACE + { + const SMDS_MeshNode* node = nIt->next(); + TSmoothNode & sNode = smooNoMap[ node ]; + sNode._uv = myHelper->GetNodeUV( quad->face, node ); + sNode._xyz = SMESH_TNodeXYZ( node ); + if ( fixedNodes.count( node )) + continue; // fixed - no triangles + + // set sNode._triangles + SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face ); + while ( fIt->more() ) + { + const SMDS_MeshElement* face = fIt->next(); + const int nbN = face->NbCornerNodes(); + const int nInd = face->GetNodeIndex( node ); + const int prevInd = myHelper->WrapIndex( nInd - 1, nbN ); + const int nextInd = myHelper->WrapIndex( nInd + 1, nbN ); + const SMDS_MeshNode* prevNode = face->GetNode( prevInd ); + const SMDS_MeshNode* nextNode = face->GetNode( nextInd ); + sNode._triangles.push_back( TTriangle( & smooNoMap[ prevNode ], + & smooNoMap[ nextNode ])); + } + } + // set _uv of smooth nodes on FACE boundary + set< StdMeshers_FaceSide* > sidesOnEdge; + list< FaceQuadStruct::Ptr >::iterator q = myQuadList.begin(); + for ( ; q != myQuadList.end() ; ++q ) + for ( size_t i = 0; i < (*q)->side.size(); ++i ) + if ( ! (*q)->side[i].grid->Edge(0).IsNull() && + //(*q)->nbNodeOut( i ) == 0 && + sidesOnEdge.insert( (*q)->side[i].grid.get() ).second ) + { + const vector& uvVec = (*q)->side[i].grid->GetUVPtStruct(); + for ( unsigned j = 0; j < uvVec.size(); ++j ) + { + TSmoothNode & sNode = smooNoMap[ uvVec[j].node ]; + sNode._uv = uvVec[j].UV(); + sNode._xyz = SMESH_TNodeXYZ( uvVec[j].node ); + } + } - // define refernce orientation in 2D - TNo2SmooNoMap::iterator n2sn = smooNoMap.begin(); - for ( ; n2sn != smooNoMap.end(); ++n2sn ) - if ( !n2sn->second._triangles.empty() ) - break; - if ( n2sn == smooNoMap.end() ) return; - const TSmoothNode & sampleNode = n2sn->second; - const bool refForward = ( sampleNode._triangles[0].IsForward( sampleNode._uv )); + // define reference orientation in 2D + TNo2SmooNoMap::iterator n2sn = smooNoMap.begin(); + for ( ; n2sn != smooNoMap.end(); ++n2sn ) + if ( !n2sn->second._triangles.empty() ) + break; + if ( n2sn == smooNoMap.end() ) return; + const TSmoothNode & sampleNode = n2sn->second; + const bool refForward = ( sampleNode._triangles[0].IsForward( sampleNode._uv )); - // Smoothing + // Smoothing - for ( int iLoop = 0; iLoop < 5; ++iLoop ) - { - for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn ) + for ( int iLoop = 0; iLoop < 5; ++iLoop ) { - TSmoothNode& sNode = n2sn->second; - if ( sNode._triangles.empty() ) - continue; // not movable node + for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn ) + { + TSmoothNode& sNode = n2sn->second; + if ( sNode._triangles.empty() ) + continue; // not movable node - gp_XY newUV; - bool isValid = false; - bool use3D = ( iLoop > 2 ); // 3 loops in 2D and 2, in 3D + gp_XY newUV; + bool isValid = false; + bool use3D = ( iLoop > 2 ); // 3 loops in 2D and 2, in 3D - if ( use3D ) - { - // compute a new XYZ - gp_XYZ newXYZ (0,0,0); - for ( unsigned i = 0; i < sNode._triangles.size(); ++i ) - newXYZ += sNode._triangles[i]._n1->_xyz; - newXYZ /= sNode._triangles.size(); - - // compute a new UV by projection - proj.Perform( newXYZ ); - isValid = ( proj.IsDone() && proj.NbPoints() > 0 ); - if ( isValid ) + if ( use3D ) { + // compute a new XYZ + gp_XYZ newXYZ (0,0,0); + for ( size_t i = 0; i < sNode._triangles.size(); ++i ) + newXYZ += sNode._triangles[i]._n1->_xyz; + newXYZ /= sNode._triangles.size(); + + // compute a new UV by projection + newUV = surface->NextValueOfUV( sNode._uv, newXYZ, 10*tol ).XY(); + // check validity of the newUV - Quantity_Parameter u,v; - proj.LowerDistanceParameters( u, v ); - newUV.SetCoord( u, v ); + for ( size_t i = 0; i < sNode._triangles.size() && isValid; ++i ) + isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward ); + } + if ( !isValid ) + { + // compute a new UV by averaging + newUV.SetCoord(0.,0.); + for ( unsigned i = 0; i < sNode._triangles.size(); ++i ) + newUV += sNode._triangles[i]._n1->_uv; + newUV /= sNode._triangles.size(); + + // check validity of the newUV + isValid = true; for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i ) isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward ); } - } - if ( !isValid ) - { - // compute a new UV by averaging - newUV.SetCoord(0.,0.); - for ( unsigned i = 0; i < sNode._triangles.size(); ++i ) - newUV += sNode._triangles[i]._n1->_uv; - newUV /= sNode._triangles.size(); - - // check validity of the newUV - isValid = true; - for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i ) - isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward ); - } - if ( isValid ) - { - sNode._uv = newUV; - sNode._xyz = surface->Value( newUV.X(), newUV.Y() ).XYZ(); + if ( isValid ) + { + sNode._uv = newUV; + sNode._xyz = surface->Value( newUV ).XYZ(); + } } } - } - // Set new XYZ to the smoothed nodes + // Set new XYZ to the smoothed nodes - for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn ) - { - TSmoothNode& sNode = n2sn->second; - if ( sNode._triangles.empty() ) - continue; // not movable node + for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn ) + { + TSmoothNode& sNode = n2sn->second; + if ( sNode._triangles.empty() ) + continue; // not movable node - SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( n2sn->first ); - gp_Pnt xyz = surface->Value( sNode._uv.X(), sNode._uv.Y() ); - meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() ); + SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( n2sn->first ); + gp_Pnt xyz = surface->Value( sNode._uv ); + meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() ); - // store the new UV - node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( sNode._uv.X(), sNode._uv.Y() ))); + // store the new UV + node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( sNode._uv.X(), sNode._uv.Y() ))); + } } // Move medium nodes in quadratic mesh @@ -4023,16 +4509,17 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad) if ( node->getshapeId() != myHelper->GetSubShapeID() ) continue; // medium node is on EDGE or VERTEX - gp_XY uv1 = myHelper->GetNodeUV( geomFace, link.node1(), node ); - gp_XY uv2 = myHelper->GetNodeUV( geomFace, link.node2(), node ); + gp_XYZ pm = 0.5 * ( SMESH_TNodeXYZ( link.node1() ) + SMESH_TNodeXYZ( link.node2() )); + gp_XY uvm = myHelper->GetNodeUV( quad->face, node ); + + gp_Pnt2d uv = surface->NextValueOfUV( uvm, pm, 10*tol ); + gp_Pnt xyz = surface->Value( uv ); - gp_XY uv = myHelper->GetMiddleUV( surface, uv1, uv2 ); node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() ))); - - gp_Pnt xyz = surface->Value( uv.X(), uv.Y() ); meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() ); } } + return; } //================================================================================ @@ -4059,11 +4546,11 @@ bool StdMeshers_Quadrangle_2D::check() { TError err; TSideVector wireVec = - StdMeshers_FaceSide::GetFaceWires( geomFace, *myHelper->GetMesh(), true, err ); + StdMeshers_FaceSide::GetFaceWires( geomFace, *myHelper->GetMesh(), true, err, myHelper ); StdMeshers_FaceSidePtr wire = wireVec[0]; // find a right angle VERTEX - int iVertex; + int iVertex = 0; double maxAngle = -1e100; for ( int i = 0; i < wire->NbEdges(); ++i ) { @@ -4128,11 +4615,18 @@ bool StdMeshers_Quadrangle_2D::check() if ( myHelper->HasSeam() ) for ( int i = 0; i < nbN && !nInFace; ++i ) if ( !myHelper->IsSeamShape( nn[i]->getshapeId() )) + { nInFace = nn[i]; + gp_XY uv = myHelper->GetNodeUV( geomFace, nInFace ); + if ( myHelper->IsOnSeam( uv )) + nInFace = NULL; + } + toCheckUV = true; for ( int i = 0; i < nbN; ++i ) uv[ i ] = myHelper->GetNodeUV( geomFace, nn[i], nInFace, &toCheckUV ); + bool isBad = false; switch ( nbN ) { case 4: { @@ -4145,28 +4639,45 @@ bool StdMeshers_Quadrangle_2D::check() if ( sign1 * sign2 < 0 ) continue; // this should not happen } - if ( sign1 * okSign < 0 ) - badFaces.push_back ( f ); + isBad = ( sign1 * okSign < 0 ); break; } case 3: { double sign = getArea( uv[0], uv[1], uv[2] ); - if ( sign * okSign < 0 ) - badFaces.push_back ( f ); + isBad = ( sign * okSign < 0 ); break; } default:; } + + // if ( isBad && myHelper->HasRealSeam() ) + // { + // // detect a case where a face intersects the seam + // for ( int iPar = 1; iPar < 3; ++iPar ) + // if ( iPar & myHelper->GetPeriodicIndex() ) + // { + // double min = uv[0].Coord( iPar ), max = uv[0].Coord( iPar ); + // for ( int i = 1; i < nbN; ++i ) + // { + // min = Min( min, uv[i].Coord( iPar )); + // max = Max( max, uv[i].Coord( iPar )); + // } + // } + // } + if ( isBad ) + badFaces.push_back ( f ); } if ( !badFaces.empty() ) { SMESH_subMesh* fSM = myHelper->GetMesh()->GetSubMesh( geomFace ); SMESH_ComputeErrorPtr& err = fSM->GetComputeError(); - err.reset ( new SMESH_ComputeError( COMPERR_ALGO_FAILED, - "Inverted elements generated")); - err->myBadElements.swap( badFaces ); + SMESH_BadInputElements* badElems = + new SMESH_BadInputElements( meshDS, COMPERR_ALGO_FAILED, + "Inverted elements generated"); + badElems->myBadElements.swap( badFaces ); + err.reset( badElems ); return !isOK; } @@ -4174,372 +4685,6 @@ bool StdMeshers_Quadrangle_2D::check() return isOK; } -/*//================================================================================ -/*! - * \brief Finds vertices at the most sharp face corners - * \param [in] theFace - the FACE - * \param [in,out] theWire - the ordered edges of the face. It can be modified to - * have the first VERTEX of the first EDGE in \a vertices - * \param [out] theVertices - the found corner vertices in the order corresponding to - * the order of EDGEs in \a theWire - * \param [out] theNbDegenEdges - nb of degenerated EDGEs in theFace - * \param [in] theConsiderMesh - if \c true, only meshed VERTEXes are considered - * as possible corners - * \return int - number of quad sides found: 0, 3 or 4 - */ -//================================================================================ - -int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, - SMESH_Mesh & theMesh, - std::list& theWire, - std::vector& theVertices, - int & theNbDegenEdges, - const bool theConsiderMesh) -{ - theNbDegenEdges = 0; - - SMESH_MesherHelper helper( theMesh ); - - // sort theVertices by angle - multimap vertexByAngle; - TopTools_DataMapOfShapeReal angleByVertex; - TopoDS_Edge prevE = theWire.back(); - if ( SMESH_Algo::isDegenerated( prevE )) - { - list::reverse_iterator edge = ++theWire.rbegin(); - while ( SMESH_Algo::isDegenerated( *edge )) - ++edge; - if ( edge == theWire.rend() ) - return false; - prevE = *edge; - } - list::iterator edge = theWire.begin(); - for ( ; edge != theWire.end(); ++edge ) - { - if ( SMESH_Algo::isDegenerated( *edge )) - { - ++theNbDegenEdges; - continue; - } - TopoDS_Vertex v = helper.IthVertex( 0, *edge ); - if ( !theConsiderMesh || SMESH_Algo::VertexNode( v, helper.GetMeshDS() )) - { - double angle = SMESH_MesherHelper::GetAngle( prevE, *edge, theFace, v ); - vertexByAngle.insert( make_pair( angle, v )); - angleByVertex.Bind( v, angle ); - } - prevE = *edge; - } - - // find out required nb of corners (3 or 4) - int nbCorners = 4; - TopoDS_Shape triaVertex = helper.GetMeshDS()->IndexToShape( myTriaVertexID ); - if ( !triaVertex.IsNull() && - triaVertex.ShapeType() == TopAbs_VERTEX && - helper.IsSubShape( triaVertex, theFace ) && - ( vertexByAngle.size() != 4 || vertexByAngle.begin()->first < 5 * M_PI/180. )) - nbCorners = 3; - else - triaVertex.Nullify(); - - // check nb of available corners - if ( nbCorners == 3 ) - { - if ( vertexByAngle.size() < 3 ) - return error(COMPERR_BAD_SHAPE, - TComm("Face must have 3 sides but not ") << vertexByAngle.size() ); - } - else - { - if ( vertexByAngle.size() == 3 && theNbDegenEdges == 0 ) - { - if ( myTriaVertexID < 1 ) - return error(COMPERR_BAD_PARMETERS, - "No Base vertex provided for a trilateral geometrical face"); - - TComm comment("Invalid Base vertex: "); - comment << myTriaVertexID << " its ID is not among [ "; - multimap::iterator a2v = vertexByAngle.begin(); - comment << helper.GetMeshDS()->ShapeToIndex( a2v->second ) << ", "; a2v++; - comment << helper.GetMeshDS()->ShapeToIndex( a2v->second ) << ", "; a2v++; - comment << helper.GetMeshDS()->ShapeToIndex( a2v->second ) << " ]"; - return error(COMPERR_BAD_PARMETERS, comment ); - } - if ( vertexByAngle.size() + ( theNbDegenEdges > 0 ) < 4 && - vertexByAngle.size() + theNbDegenEdges != 4 ) - return error(COMPERR_BAD_SHAPE, - TComm("Face must have 4 sides but not ") << vertexByAngle.size() ); - } - - // put all corner vertices in a map - TopTools_MapOfShape vMap; - if ( nbCorners == 3 ) - vMap.Add( triaVertex ); - multimap::reverse_iterator a2v = vertexByAngle.rbegin(); - for ( int iC = 0; a2v != vertexByAngle.rend() && iC < nbCorners; ++a2v, ++iC ) - vMap.Add( (*a2v).second ); - - // check if there are possible variations in choosing corners - bool haveVariants = false; - if ( vertexByAngle.size() > nbCorners ) - { - double lostAngle = a2v->first; - double lastAngle = ( --a2v, a2v->first ); - haveVariants = ( lostAngle * 1.1 >= lastAngle ); - } - - const double angleTol = 5.* M_PI/180; - myCheckOri = ( vertexByAngle.size() > nbCorners || - vertexByAngle.begin()->first < angleTol ); - - // make theWire begin from a corner vertex or triaVertex - if ( nbCorners == 3 ) - while ( !triaVertex.IsSame( ( helper.IthVertex( 0, theWire.front() ))) || - SMESH_Algo::isDegenerated( theWire.front() )) - theWire.splice( theWire.end(), theWire, theWire.begin() ); - else - while ( !vMap.Contains( helper.IthVertex( 0, theWire.front() )) || - SMESH_Algo::isDegenerated( theWire.front() )) - theWire.splice( theWire.end(), theWire, theWire.begin() ); - - // fill the result vector and prepare for its refinement - theVertices.clear(); - vector< double > angles; - vector< TopoDS_Edge > edgeVec; - vector< int > cornerInd, nbSeg; - int nbSegTot = 0; - angles .reserve( vertexByAngle.size() ); - edgeVec.reserve( vertexByAngle.size() ); - nbSeg .reserve( vertexByAngle.size() ); - cornerInd.reserve( nbCorners ); - for ( edge = theWire.begin(); edge != theWire.end(); ++edge ) - { - if ( SMESH_Algo::isDegenerated( *edge )) - continue; - TopoDS_Vertex v = helper.IthVertex( 0, *edge ); - bool isCorner = vMap.Contains( v ); - if ( isCorner ) - { - theVertices.push_back( v ); - cornerInd.push_back( angles.size() ); - } - angles .push_back( angleByVertex.IsBound( v ) ? angleByVertex( v ) : -M_PI ); - edgeVec.push_back( *edge ); - if ( theConsiderMesh && haveVariants ) - { - if ( SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( *edge )) - nbSeg.push_back( sm->NbNodes() + 1 ); - else - nbSeg.push_back( 0 ); - nbSegTot += nbSeg.back(); - } - } - - // refine the result vector - make sides equal by length if - // there are several equal angles - if ( haveVariants ) - { - if ( nbCorners == 3 ) - angles[0] = 2 * M_PI; // not to move the base triangle VERTEX - - // here we refer to VERTEX'es and EDGEs by indices in angles and edgeVec vectors - typedef int TGeoIndex; - - // for each vertex find a vertex till which there are nbSegHalf segments - const int nbSegHalf = ( nbSegTot % 2 || nbCorners == 3 ) ? 0 : nbSegTot / 2; - vector< TGeoIndex > halfDivider( angles.size(), -1 ); - int nbHalfDividers = 0; - if ( nbSegHalf ) - { - // get min angle of corners - double minAngle = 10.; - for ( size_t iC = 0; iC < cornerInd.size(); ++iC ) - minAngle = Min( minAngle, angles[ cornerInd[ iC ]]); - - // find halfDivider's - for ( TGeoIndex iV1 = 0; iV1 < TGeoIndex( angles.size() ); ++iV1 ) - { - int nbSegs = 0; - TGeoIndex iV2 = iV1; - do { - nbSegs += nbSeg[ iV2 ]; - iV2 = helper.WrapIndex( iV2 + 1, nbSeg.size() ); - } while ( nbSegs < nbSegHalf ); - - if ( nbSegs == nbSegHalf && - angles[ iV1 ] + angleTol >= minAngle && - angles[ iV2 ] + angleTol >= minAngle ) - { - halfDivider[ iV1 ] = iV2; - ++nbHalfDividers; - } - } - } - - set< TGeoIndex > refinedCorners, treatedCorners; - for ( size_t iC = 0; iC < cornerInd.size(); ++iC ) - { - TGeoIndex iV = cornerInd[iC]; - if ( !treatedCorners.insert( iV ).second ) - continue; - list< TGeoIndex > equVerts; // inds of vertices that can become corners - equVerts.push_back( iV ); - int nbC[2] = { 0, 0 }; - // find equal angles backward and forward from the iV-th corner vertex - for ( int isFwd = 0; isFwd < 2; ++isFwd ) - { - int dV = isFwd ? +1 : -1; - int iCNext = helper.WrapIndex( iC + dV, cornerInd.size() ); - TGeoIndex iVNext = helper.WrapIndex( iV + dV, angles.size() ); - while ( iVNext != iV ) - { - bool equal = Abs( angles[iV] - angles[iVNext] ) < angleTol; - if ( equal ) - equVerts.insert( isFwd ? equVerts.end() : equVerts.begin(), iVNext ); - if ( iVNext == cornerInd[ iCNext ]) - { - if ( !equal ) - { - if ( angles[iV] < angles[iVNext] ) - refinedCorners.insert( iVNext ); - break; - } - nbC[ isFwd ]++; - treatedCorners.insert( cornerInd[ iCNext ] ); - iCNext = helper.WrapIndex( iCNext + dV, cornerInd.size() ); - } - iVNext = helper.WrapIndex( iVNext + dV, angles.size() ); - } - if ( iVNext == iV ) - break; // all angles equal - } - - const bool allCornersSame = ( nbC[0] == 3 ); - if ( allCornersSame && nbHalfDividers > 0 ) - { - // select two halfDivider's as corners - TGeoIndex hd1, hd2 = -1; - int iC2; - for ( iC2 = 0; iC2 < cornerInd.size() && hd2 < 0; ++iC2 ) - { - hd1 = cornerInd[ iC2 ]; - hd2 = halfDivider[ hd1 ]; - if ( std::find( equVerts.begin(), equVerts.end(), hd2 ) == equVerts.end() ) - hd2 = -1; // hd2-th vertex can't become a corner - else - break; - } - if ( hd2 >= 0 ) - { - angles[ hd1 ] = 2 * M_PI; // make hd1-th vertex no more "equal" - angles[ hd2 ] = 2 * M_PI; - refinedCorners.insert( hd1 ); - refinedCorners.insert( hd2 ); - treatedCorners = refinedCorners; - // update cornerInd - equVerts.push_front( equVerts.back() ); - equVerts.push_back( equVerts.front() ); - list< TGeoIndex >::iterator hdPos = - std::find( equVerts.begin(), equVerts.end(), hd2 ); - if ( hdPos == equVerts.end() ) break; - cornerInd[ helper.WrapIndex( iC2 + 0, cornerInd.size()) ] = hd1; - cornerInd[ helper.WrapIndex( iC2 + 1, cornerInd.size()) ] = *( --hdPos ); - cornerInd[ helper.WrapIndex( iC2 + 2, cornerInd.size()) ] = hd2; - cornerInd[ helper.WrapIndex( iC2 + 3, cornerInd.size()) ] = *( ++hdPos, ++hdPos ); - - theVertices[ 0 ] = helper.IthVertex( 0, edgeVec[ cornerInd[0] ]); - theVertices[ 1 ] = helper.IthVertex( 0, edgeVec[ cornerInd[1] ]); - theVertices[ 2 ] = helper.IthVertex( 0, edgeVec[ cornerInd[2] ]); - theVertices[ 3 ] = helper.IthVertex( 0, edgeVec[ cornerInd[3] ]); - iC = -1; - continue; - } - } - - // move corners to make sides equal by length - int nbEqualV = equVerts.size(); - int nbExcessV = nbEqualV - ( 1 + nbC[0] + nbC[1] ); - if ( nbExcessV > 0 ) // there is nbExcessV vertices that can become corners - { - // calculate normalized length of each "side" enclosed between neighbor equVerts - vector< double > accuLength; - double totalLen = 0; - vector< TGeoIndex > evVec( equVerts.begin(), equVerts.end() ); - int iEV = 0; - TGeoIndex iE = cornerInd[ helper.WrapIndex( iC - nbC[0] - 1, cornerInd.size() )]; - TGeoIndex iEEnd = cornerInd[ helper.WrapIndex( iC + nbC[1] + 1, cornerInd.size() )]; - while ( accuLength.size() < nbEqualV + int( !allCornersSame ) ) - { - // accumulate length of edges before iEV-th equal vertex - accuLength.push_back( totalLen ); - do { - accuLength.back() += SMESH_Algo::EdgeLength( edgeVec[ iE ]); - iE = helper.WrapIndex( iE + 1, edgeVec.size()); - if ( iEV < evVec.size() && iE == evVec[ iEV ] ) { - iEV++; - break; // equal vertex reached - } - } - while( iE != iEEnd ); - totalLen = accuLength.back(); - } - accuLength.resize( equVerts.size() ); - for ( size_t iS = 0; iS < accuLength.size(); ++iS ) - accuLength[ iS ] /= totalLen; - - // find equVerts most close to the ideal sub-division of all sides - int iBestEV = 0; - int iCorner = helper.WrapIndex( iC - nbC[0], cornerInd.size() ); - int nbSides = Min( nbCorners, 2 + nbC[0] + nbC[1] ); - for ( int iS = 1; iS < nbSides; ++iS, ++iBestEV ) - { - double idealLen = iS / double( nbSides ); - double d, bestDist = 2.; - for ( iEV = iBestEV; iEV < accuLength.size(); ++iEV ) - { - d = Abs( idealLen - accuLength[ iEV ]); - - // take into account presence of a coresponding halfDivider - const double cornerWgt = 0.5 / nbSides; - const double vertexWgt = 0.25 / nbSides; - TGeoIndex hd = halfDivider[ evVec[ iEV ]]; - if ( hd < 0 ) - d += vertexWgt; - else if( refinedCorners.count( hd )) - d -= cornerWgt; - else - d -= vertexWgt; - - // choose vertex with the best d - if ( d < bestDist ) - { - bestDist = d; - iBestEV = iEV; - } - } - if ( iBestEV > iS-1 + nbExcessV ) - iBestEV = iS-1 + nbExcessV; - theVertices[ iCorner ] = helper.IthVertex( 0, edgeVec[ evVec[ iBestEV ]]); - refinedCorners.insert( evVec[ iBestEV ]); - iCorner = helper.WrapIndex( iCorner + 1, cornerInd.size() ); - } - - } // if ( nbExcessV > 0 ) - else - { - refinedCorners.insert( cornerInd[ iC ]); - } - } // loop on cornerInd - - // make theWire begin from the cornerInd[0]-th EDGE - while ( !theWire.front().IsSame( edgeVec[ cornerInd[0] ])) - theWire.splice( theWire.begin(), theWire, --theWire.end() ); - - } // if ( haveVariants ) - - return nbCorners; -} - //================================================================================ /*! * \brief Constructor of a side of quad @@ -4547,7 +4692,7 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, //================================================================================ FaceQuadStruct::Side::Side(StdMeshers_FaceSidePtr theGrid) - : grid(theGrid), nbNodeOut(0), from(0), to(theGrid ? theGrid->NbPoints() : 0 ), di(1) + : grid(theGrid), from(0), to(theGrid ? theGrid->NbPoints() : 0 ), di(1), nbNodeOut(0) { } @@ -4605,49 +4750,66 @@ bool StdMeshers_Quadrangle_2D::getEnforcedUV() Standard_Real u1,u2,v1,v2; const TopoDS_Face& face = TopoDS::Face( myHelper->GetSubShape() ); const double tol = BRep_Tool::Tolerance( face ); - Handle(Geom_Surface) surf = BRep_Tool::Surface( face ); - surf->Bounds( u1,u2,v1,v2 ); - GeomAPI_ProjectPointOnSurf project; - project.Init(surf, u1,u2, v1,v2, tol ); + Handle(ShapeAnalysis_Surface) project = myHelper->GetSurface( face ); + project->Bounds( u1,u2,v1,v2 ); Bnd_Box bbox; BRepBndLib::Add( face, bbox ); double farTol = 0.01 * sqrt( bbox.SquareExtent() ); + // get internal VERTEXes of the FACE to use them instead of equal points + typedef map< pair< double, double >, TopoDS_Vertex > TUV2VMap; + TUV2VMap uv2intV; + for ( TopExp_Explorer vExp( face, TopAbs_VERTEX, TopAbs_EDGE ); vExp.More(); vExp.Next() ) + { + TopoDS_Vertex v = TopoDS::Vertex( vExp.Current() ); + gp_Pnt2d uv = project->ValueOfUV( BRep_Tool::Pnt( v ), tol ); + uv2intV.insert( make_pair( make_pair( uv.X(), uv.Y() ), v )); + } + for ( size_t iP = 0; iP < points.size(); ++iP ) { - project.Perform( points[ iP ]); - if ( !project.IsDone() ) - { - if ( isStrictCheck && iP < nbPoints ) - return error - (TComm("Projection of an enforced point to the face failed - (") - << points[ iP ].X() << ", "<< points[ iP ].Y() << ", "<< points[ iP ].Z() << " )"); - continue; - } - if ( project.LowerDistance() > farTol ) + gp_Pnt2d uv = project->ValueOfUV( points[ iP ], tol ); + if ( project->Gap() > farTol ) { if ( isStrictCheck && iP < nbPoints ) return error (COMPERR_BAD_PARMETERS, TComm("An enforced point is too far from the face, dist = ") - << project.LowerDistance() << " - (" + << points[ iP ].Distance( project->Value( uv )) << " - (" << points[ iP ].X() << ", "<< points[ iP ].Y() << ", "<< points[ iP ].Z() << " )"); continue; } - Quantity_Parameter u, v; - project.LowerDistanceParameters(u, v); - gp_Pnt2d uv( u, v ); BRepClass_FaceClassifier clsf ( face, uv, tol ); switch ( clsf.State() ) { case TopAbs_IN: { - double edgeDist = ( Min( Abs( u - u1 ), Abs( u - u2 )) + - Min( Abs( v - v1 ), Abs( v - v2 ))); + double edgeDist = ( Min( Abs( uv.X() - u1 ), Abs( uv.X() - u2 )) + + Min( Abs( uv.Y() - v1 ), Abs( uv.Y() - v2 ))); ForcedPoint fp; fp.uv = uv.XY(); fp.xyz = points[ iP ].XYZ(); if ( iP >= nbPoints ) fp.vertex = TopoDS::Vertex( vMap( iP - nbPoints + 1 )); + TUV2VMap::iterator uv2v = uv2intV.lower_bound( make_pair( uv.X()-tol, uv.Y()-tol )); + for ( ; uv2v != uv2intV.end() && uv2v->first.first <= uv.X()+tol; ++uv2v ) + if ( uv.SquareDistance( gp_Pnt2d( uv2v->first.first, uv2v->first.second )) < tol*tol ) + { + fp.vertex = uv2v->second; + break; + } + + fp.node = 0; + if ( myHelper->IsSubShape( fp.vertex, myHelper->GetMesh() )) + { + SMESH_subMesh* sm = myHelper->GetMesh()->GetSubMesh( fp.vertex ); + sm->ComputeStateEngine( SMESH_subMesh::COMPUTE ); + fp.node = SMESH_Algo::VertexNode( fp.vertex, myHelper->GetMeshDS() ); + } + else + { + fp.node = myHelper->AddNode( fp.xyz.X(), fp.xyz.Y(), fp.xyz.Z(), + 0, fp.uv.X(), fp.uv.Y() ); + } sortedFP.insert( make_pair( edgeDist, fp )); break; } @@ -4707,8 +4869,6 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() quadsBySide[ (*quadIt)->side[iSide] ].push_back( *quadIt ); } - SMESH_Mesh* mesh = myHelper->GetMesh(); - SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); const TopoDS_Face& face = TopoDS::Face( myHelper->GetSubShape() ); Handle(Geom_Surface) surf = BRep_Tool::Surface( face ); @@ -4716,7 +4876,7 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() { bool isNodeEnforced = false; - // look for a quad enclosing a enforced point + // look for a quad enclosing an enforced point for ( quadIt = myQuadList.begin(); quadIt != myQuadList.end(); ++quadIt ) { FaceQuadStruct::Ptr quad = *quadIt; @@ -4771,8 +4931,9 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() } // make a node of a side forced vector& points = (vector&) side.GetUVPtStruct(); - points[ sideNodeIndex ].u = myForcedPnts[ iFP ].U(); - points[ sideNodeIndex ].v = myForcedPnts[ iFP ].V(); + points[ sideNodeIndex ].u = myForcedPnts[ iFP ].U(); + points[ sideNodeIndex ].v = myForcedPnts[ iFP ].V(); + points[ sideNodeIndex ].node = myForcedPnts[ iFP ].node; updateSideUV( side, sideNodeIndex, quadsBySide ); @@ -4819,6 +4980,9 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() FaceQuadStruct::Ptr newQuad = myQuadList.back(); FaceQuadStruct::Side& newSide = newQuad->side[ iNewSide ]; + vector& points = (vector&) newSide.GetUVPtStruct(); + points[ indForced ].node = myForcedPnts[ iFP ].node; + newSide.forced_nodes.insert( indForced ); quad->side[( iNewSide+2 ) % 4 ].forced_nodes.insert( indForced ); @@ -4855,6 +5019,7 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() << myForcedPnts[iFP].xyz.Y() << ", " << myForcedPnts[iFP].xyz.Z() << " )"); } + myNeedSmooth = true; } // loop on enforced points @@ -4872,25 +5037,31 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() if ( quadVec.size() <= 1 ) continue; // outer side - bool missedNodesOnSide = false; const vector& points = side.grid->GetUVPtStruct(); for ( size_t iC = 0; iC < side.contacts.size(); ++iC ) { + if ( side.contacts[iC].point < side.from || + side.contacts[iC].point >= side.to ) + continue; + if ( side.contacts[iC].other_point < side.contacts[iC].other_side->from || + side.contacts[iC].other_point >= side.contacts[iC].other_side->to ) + continue; const vector& oGrid = side.contacts[iC].other_side->grid->GetUVPtStruct(); const UVPtStruct& uvPt = points[ side.contacts[iC].point ]; - if ( side.contacts[iC].other_point >= oGrid.size() || - side.contacts[iC].point >= points.size() ) + if ( side.contacts[iC].other_point >= (int) oGrid .size() || + side.contacts[iC].point >= (int) points.size() ) throw SALOME_Exception( "StdMeshers_Quadrangle_2D::addEnforcedNodes(): wrong contact" ); if ( oGrid[ side.contacts[iC].other_point ].node ) (( UVPtStruct& ) uvPt).node = oGrid[ side.contacts[iC].other_point ].node; } + + bool missedNodesOnSide = false; for ( size_t iP = 0; iP < points.size(); ++iP ) if ( !points[ iP ].node ) { UVPtStruct& uvPnt = ( UVPtStruct& ) points[ iP ]; - gp_Pnt P = surf->Value( uvPnt.u, uvPnt.v ); - uvPnt.node = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace( uvPnt.node, myHelper->GetSubShapeID(), uvPnt.u, uvPnt.v ); + gp_Pnt P = surf->Value( uvPnt.u, uvPnt.v ); + uvPnt.node = myHelper->AddNode(P.X(), P.Y(), P.Z(), 0, uvPnt.u, uvPnt.v ); missedNodesOnSide = true; } if ( missedNodesOnSide ) @@ -4950,8 +5121,10 @@ int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J) newQuad->side[ QUAD_TOP_SIDE ].from = iTop; newQuad->name = ( TComm("Right of I=") << I ); - quad->side[ QUAD_BOTTOM_SIDE ].to = iBot + 1; - quad->side[ QUAD_TOP_SIDE ].to = iTop + 1; + bool bRev = quad->side[ QUAD_BOTTOM_SIDE ].IsReversed(); + bool tRev = quad->side[ QUAD_TOP_SIDE ].IsReversed(); + quad->side[ QUAD_BOTTOM_SIDE ].to = iBot + ( bRev ? -1 : +1 ); + quad->side[ QUAD_TOP_SIDE ].to = iTop + ( tRev ? -1 : +1 ); quad->uv_grid.clear(); return QUAD_LEFT_SIDE; @@ -4985,8 +5158,10 @@ int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J) // << " L " << &quad->side[ QUAD_LEFT_SIDE ] << " "<< quad->side[ QUAD_LEFT_SIDE].NbPoints() // << " R " << &quad->side[ QUAD_RIGHT_SIDE ] << " "<< quad->side[ QUAD_RIGHT_SIDE].NbPoints()<< endl; - newQuad->side[ QUAD_RIGHT_SIDE ].to = iRgt+1; - newQuad->side[ QUAD_LEFT_SIDE ].to = iLft+1; + bool rRev = newQuad->side[ QUAD_RIGHT_SIDE ].IsReversed(); + bool lRev = newQuad->side[ QUAD_LEFT_SIDE ].IsReversed(); + newQuad->side[ QUAD_RIGHT_SIDE ].to = iRgt + ( rRev ? -1 : +1 ); + newQuad->side[ QUAD_LEFT_SIDE ].to = iLft + ( lRev ? -1 : +1 ); newQuad->name = ( TComm("Below J=") << J ); quad->side[ QUAD_RIGHT_SIDE ].from = iRgt; @@ -5030,9 +5205,9 @@ void StdMeshers_Quadrangle_2D::updateSideUV( FaceQuadStruct::Side& side, return; } - const int iFrom = Min ( iForced, *iNext ); - const int iTo = Max ( iForced, *iNext ) + 1; - const int sideSize = iTo - iFrom; + const int iFrom = Min ( iForced, *iNext ); + const int iTo = Max ( iForced, *iNext ) + 1; + const size_t sideSize = iTo - iFrom; vector points[4]; // side points of a temporary quad @@ -5042,7 +5217,7 @@ void StdMeshers_Quadrangle_2D::updateSideUV( FaceQuadStruct::Side& side, for ( int is2nd = 0; is2nd < 2; ++is2nd ) { points[ is2nd ].reserve( sideSize ); - int nbLoops = 0; + size_t nbLoops = 0; while ( points[is2nd].size() < sideSize ) { int iCur = iFrom + points[is2nd].size() - int( !points[is2nd].empty() ); @@ -5057,6 +5232,8 @@ void StdMeshers_Quadrangle_2D::updateSideUV( FaceQuadStruct::Side& side, for ( iS = 0; iS < q->side.size(); ++iS ) if ( side.grid == q->side[ iS ].grid ) break; + if ( iS == q->side.size() ) + continue; bool isOut; if ( !q->side[ iS ].IsReversed() ) isOut = ( q->side[ iS ].from > iCur || q->side[ iS ].to-1 <= iCur ); @@ -5536,9 +5713,11 @@ bool FaceQuadStruct::Side::IsForced( int nodeIndex ) const void FaceQuadStruct::Side::AddContact( int ip, Side* side, int iop ) { - if ( ip >= GetUVPtStruct().size() || - iop >= side->GetUVPtStruct().size() ) + if ( ip >= (int) GetUVPtStruct().size() || + iop >= (int) side->GetUVPtStruct().size() ) throw SALOME_Exception( "FaceQuadStruct::Side::AddContact(): wrong point" ); + if ( ip < from || ip >= to ) + return; { contacts.resize( contacts.size() + 1 ); Contact& c = contacts.back();