X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=e9e6ad82744158deec6e4d168e76324b3eb57343;hp=7adf2294f674da4706867b3d61664c14bf8891ce;hb=373690534aa482a4daea68f2a5a14fa9de8a9a1b;hpb=076deefb2bbc440cf0c15f6059182f4722d7259c diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index 7adf2294f..e9e6ad827 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 @@ -67,9 +67,7 @@ #ifndef StdMeshers_Array2OfNode_HeaderFile #define StdMeshers_Array2OfNode_HeaderFile typedef const SMDS_MeshNode* SMDS_MeshNodePtr; -DEFINE_BASECOLLECTION (StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr) -DEFINE_ARRAY2(StdMeshers_Array2OfNode, - StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr) +typedef NCollection_Array2 StdMeshers_Array2OfNode; #endif using namespace std; @@ -90,6 +88,7 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, myTrianglePreference(false), myTriaVertexID(-1), myNeedSmooth(false), + myCheckOri(false), myParams( NULL ), myQuadType(QUAD_STANDARD), myHelper( NULL ) @@ -133,8 +132,7 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis myParams = NULL; myQuadList.clear(); - bool isOk = true; - aStatus = SMESH_Hypothesis::HYP_OK; + aStatus = SMESH_Hypothesis::HYP_OK; const list & hyps = GetUsedHypothesis(aMesh, aShape, false); @@ -200,7 +198,9 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis } } - return isOk; + error( StdMeshers_ViscousLayers2D::CheckHypothesis( aMesh, aShape, aStatus )); + + return aStatus == HYP_OK; } //============================================================================= @@ -227,7 +227,9 @@ 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 ); if (!quad) @@ -247,7 +249,7 @@ 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))) @@ -256,7 +258,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, res = computeQuadPref( aMesh, F, quad ); } } - else if (myQuadType == QUAD_REDUCED) + else if ( myQuadType == QUAD_REDUCED ) { int n13 = n1 - n3; int n24 = n2 - n4; @@ -275,7 +277,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, "two opposite sides should have same number of segments, " "but actual number of segments is different on all sides. " "'Standard' transion has been used."); - else + else if ( ! ( n1 == n3 && n2 == n4 )) error( COMPERR_WARNING, "To use 'Reduced' transition, " "two opposite sides should have an even difference in number of segments. " @@ -294,6 +296,9 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, if ( res == COMPUTE_OK && myNeedSmooth ) smooth( quad ); + if ( res == COMPUTE_OK ) + res = check(); + return ( res == COMPUTE_OK ); } @@ -334,7 +339,7 @@ bool StdMeshers_Quadrangle_2D::computeTriangles(SMESH_Mesh& aMesh, FaceQuadStruct::Ptr newQuad = myQuadList.back(); if ( quad != newQuad ) // split done { - { + { // update left side limit till where to make triangles FaceQuadStruct::Ptr botQuad = // a bottom part ( quad->side[ QUAD_LEFT_SIDE ].from == 0 ) ? quad : newQuad; if ( botQuad->nbNodeOut( QUAD_LEFT_SIDE ) > 0 ) @@ -354,12 +359,33 @@ bool StdMeshers_Quadrangle_2D::computeTriangles(SMESH_Mesh& aMesh, { splitQuad( quad, quad->iSize-2, 0 ); } - if ( quad->nbNodeOut( QUAD_LEFT_SIDE )) + if ( quad->nbNodeOut( QUAD_LEFT_SIDE )) { splitQuad( quad, 1, 0 ); + + if ( quad->nbNodeOut( QUAD_TOP_SIDE )) + { + newQuad = myQuadList.back(); + if ( newQuad == quad ) // too narrow to split + { + // update left side limit till where to make triangles + quad->side[ QUAD_LEFT_SIDE ].to--; + } + else + { + FaceQuadStruct::Ptr leftQuad = + ( quad->side[ QUAD_BOTTOM_SIDE ].from == 0 ) ? quad : newQuad; + leftQuad->nbNodeOut( QUAD_TOP_SIDE ) = 0; + } + } } - return computeQuadDominant( aMesh, aFace ); + if ( ! computeQuadDominant( aMesh, aFace )) + return false; + + // try to fix zero-area triangles near straight-angle corners + + return true; } //================================================================================ @@ -445,10 +471,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); } } @@ -496,7 +519,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()); @@ -510,6 +533,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++) { @@ -532,8 +556,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) @@ -543,8 +566,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); @@ -558,8 +580,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; @@ -583,23 +604,63 @@ 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 left edge is out, we will stop at a second node - //if (quad->nbNodeOut(3)) stop++; - if ( quad->nbNodeOut( QUAD_RIGHT_SIDE )) - quad->UVPt( nbhoriz-1, 0 ).node = uv_e1[ nbright-2 ].node; - if ( quad->nbNodeOut( QUAD_LEFT_SIDE )) - quad->UVPt( 0, 0 ).node = uv_e3[ nbleft-2 ].node; + 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 + } + else + { + if ( quad->nbNodeOut( QUAD_RIGHT_SIDE )) + quad->UVPt( nbhoriz-1, 0 ).node = uv_e1[ nbright-2 ].node; + if ( quad->nbNodeOut( QUAD_LEFT_SIDE )) + quad->UVPt( 0, 0 ).node = uv_e3[ nbleft-2 ].node; + if ( nbright > 2 ) // there was a split at J + quad->nbNodeOut( QUAD_LEFT_SIDE ) = 0; + } + const SMDS_MeshNode *a, *b, *c, *d; + i = nbup - 1; + // avoid creating zero-area triangles near a straight-angle corner + { + a = uv_e2[i].node; + b = uv_e2[i-1].node; + c = uv_e1[nbright-2].node; + SMESH_TNodeXYZ pa( a ), pb( b ), pc( c ); + double area = 0.5 * (( pb - pa ) ^ ( pc - pa )).Modulus(); + if ( Abs( area ) < 1e-20 ) + { + --g; + d = quad->UVPt( g, nbvertic-2 ).node; + if ( myTrianglePreference ) + { + myHelper->AddFace(a, d, c); + } + else + { + if ( SMDS_MeshFace* face = myHelper->AddFace(a, b, d, c)) + { + 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 ); + } + } + --i; + } + } + } // for each node of the up edge find nearest node // in the first row of the regular grid and link them - for (i = nbup - 1; i > stop; i--) { - const SMDS_MeshNode *a, *b, *c, *d; + 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; @@ -615,7 +676,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; @@ -628,8 +689,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) @@ -638,8 +698,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); @@ -652,8 +711,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; @@ -677,6 +735,7 @@ 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 c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node; @@ -702,8 +761,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) @@ -713,8 +771,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); @@ -727,8 +784,7 @@ 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; @@ -740,8 +796,41 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, int g = nbvertic - 1; // last processed node in the grid int stop = 0; i = quad->side[ QUAD_LEFT_SIDE ].to-1; // nbleft - 1; - for (; i > stop; i--) { - const SMDS_MeshNode *a, *b, *c, *d; + + const SMDS_MeshNode *a, *b, *c, *d; + // avoid creating zero-area triangles near a straight-angle corner + { + a = uv_e3[i].node; + b = uv_e3[i-1].node; + c = quad->UVPt( 1, g ).node; + SMESH_TNodeXYZ pa( a ), pb( b ), pc( c ); + double area = 0.5 * (( pb - pa ) ^ ( pc - pa )).Modulus(); + if ( Abs( area ) < 1e-20 ) + { + --g; + d = quad->UVPt( 1, g ).node; + if ( myTrianglePreference ) + { + myHelper->AddFace(a, d, c); + } + else + { + if ( SMDS_MeshFace* face = myHelper->AddFace(a, b, d, c)) + { + 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 ); + } + } + --i; + } + } + } + for (; i > stop; i--) // loop on nodes on the left side + { a = uv_e3[i].node; b = uv_e3[i - 1].node; gp_Pnt pb (b->X(), b->Y(), b->Z()); @@ -751,12 +840,13 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, if (i == stop + 1) { // down bondary reached c = quad->uv_grid[nbhoriz*jlow + 1].node; near = jlow; - } else { + } + else { double mind = RealLast(); for (int k = g; k >= jlow; k--) { const SMDS_MeshNode *nk; if (k > jup) - nk = uv_e2[1].node; + nk = quad->uv_grid[nbhoriz*jup + 1].node; //uv_e2[1].node; else nk = quad->uv_grid[nbhoriz*k + 1].node; gp_Pnt pnk (nk->X(), nk->Y(), nk->Z()); @@ -772,18 +862,15 @@ 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) - d = uv_e2[1].node; + d = quad->uv_grid[nbhoriz*jup + 1].node; //uv_e2[1].node; else d = quad->uv_grid[nbhoriz*(near + 1) + 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); + if (!myTrianglePreference) { + myHelper->AddFace(a, b, c, d); } else { splitQuadFace(meshDS, geomFaceID, a, b, c, d); @@ -793,11 +880,10 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, for (int k = near + 1; k < g; k++) { c = quad->uv_grid[nbhoriz*k + 1].node; if (k + 1 > jup) - d = uv_e2[1].node; + 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; @@ -918,21 +1004,17 @@ bool StdMeshers_Quadrangle_2D::IsApplicable( const TopoDS_Shape & aShape, bool t int nbFoundFaces = 0; for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces ) { - TopoDS_Face aFace = TopoDS::Face(exp.Current()); - if ( aFace.Orientation() >= TopAbs_INTERNAL ) aFace.Orientation( TopAbs_FORWARD ); - - list< TopoDS_Edge > aWire; - list< int > nbEdgesInWire; - int nbWire = SMESH_Block::GetOrderedEdges (aFace, aWire, nbEdgesInWire); + const TopoDS_Shape& aFace = exp.Current(); + int nbWire = SMESH_MesherHelper::Count( aFace, TopAbs_WIRE, false ); if ( nbWire != 1 ) { if ( toCheckAll ) return false; continue; } int nbNoDegenEdges = 0; - list::iterator edge = aWire.begin(); - for ( ; edge != aWire.end(); ++edge ) { - if ( !SMESH_Algo::isDegenerated( *edge )) + TopExp_Explorer eExp( aFace, TopAbs_EDGE ); + for ( ; eExp.More() && nbNoDegenEdges < 3; eExp.Next() ) { + if ( !SMESH_Algo::isDegenerated( TopoDS::Edge( eExp.Current() ))) ++nbNoDegenEdges; } if ( toCheckAll && nbNoDegenEdges < 3 ) return false; @@ -979,7 +1061,7 @@ FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD ); const bool ignoreMediumNodes = _quadraticMesh; - // verify 1 wire only, with 4 edges + // verify 1 wire only list< TopoDS_Edge > edges; list< int > nbEdgesInWire; int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire); @@ -1035,45 +1117,40 @@ FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & { list< TopoDS_Edge > sideEdges; TopoDS_Vertex nextSideV = corners[( iSide + 1 - nbUsedDegen ) % corners.size() ]; - while ( edgeIt != edges.end() && - !nextSideV.IsSame( myHelper->IthVertex( 0, *edgeIt ))) + bool nextSideVReached = false; + do { - if ( SMESH_Algo::isDegenerated( *edgeIt ) ) + const TopoDS_Edge& edge = *edgeIt; + nextSideVReached = nextSideV.IsSame( myHelper->IthVertex( 1, edge )); + if ( SMESH_Algo::isDegenerated( edge )) { - if ( myNeedSmooth ) - { - ++edgeIt; // no side on the degenerated EDGE - } - else + if ( !myNeedSmooth ) // need to make a side on a degen edge { if ( sideEdges.empty() ) { + sideEdges.push_back( edge ); ++nbUsedDegen; - sideEdges.push_back( *edgeIt++ ); // a degenerated side - break; + nextSideVReached = true; } else { - break; // do not append a degenerated EDGE to a regular side + break; } } } else { - sideEdges.push_back( *edgeIt++ ); + sideEdges.push_back( edge ); } + ++edgeIt; } + while ( edgeIt != edges.end() && !nextSideVReached ); + if ( !sideEdges.empty() ) { - quad->side.push_back( StdMeshers_FaceSide::New( F, sideEdges, &aMesh, iSide < QUAD_TOP_SIDE, - ignoreMediumNodes, myProxyMesh )); - ++iSide; - } - else if ( !SMESH_Algo::isDegenerated( *edgeIt ) && // closed EDGE - myHelper->IthVertex( 0, *edgeIt ).IsSame( myHelper->IthVertex( 1, *edgeIt ))) - { - quad->side.push_back( StdMeshers_FaceSide::New( F, *edgeIt++, &aMesh, iSide < QUAD_TOP_SIDE, - ignoreMediumNodes, myProxyMesh)); + quad->side.push_back + ( StdMeshers_FaceSide::New( F, sideEdges, &aMesh, iSide < QUAD_TOP_SIDE, + ignoreMediumNodes, myProxyMesh )); ++iSide; } if ( quad->side.size() == 4 ) @@ -1359,6 +1436,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 ) { @@ -1505,7 +1584,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 @@ -1517,6 +1596,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) @@ -1546,7 +1627,33 @@ 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(); + } } //======================================================================= @@ -1676,10 +1783,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); @@ -1731,6 +1838,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; @@ -1816,10 +1927,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 // @@ -1882,6 +1993,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; @@ -2019,7 +2140,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 @@ -2085,10 +2206,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)); } } } @@ -2142,10 +2261,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)); } } } @@ -2215,10 +2332,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)); } } } @@ -2247,10 +2362,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)); } } } @@ -2269,7 +2382,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); @@ -2353,10 +2466,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) @@ -2542,21 +2651,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); } } @@ -2868,7 +2972,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 @@ -2968,10 +3073,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)); } } } @@ -3023,10 +3126,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)); } } } @@ -3087,10 +3188,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 @@ -3187,11 +3286,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 ); @@ -3200,7 +3298,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; @@ -3687,6 +3789,18 @@ namespace // data for smoothing double d = v1 ^ v2; return d > 1e-100; } + //================================================================================ + /*! + * \brief Returns area of a triangle + */ + //================================================================================ + + double getArea( const gp_UV uv1, const gp_UV uv2, const gp_UV uv3 ) + { + gp_XY v1 = uv1 - uv2, v2 = uv3 - uv2; + double a = v2 ^ v1; + return a; + } } //================================================================================ @@ -3737,11 +3851,11 @@ void StdMeshers_Quadrangle_2D::updateDegenUV(FaceQuadStruct::Ptr quad) uv1.v = uv2.v = 0.5 * ( uv1.v + uv2.v ); } - else if ( quad->side.size() == 4 && myQuadType == QUAD_STANDARD) + else if ( quad->side.size() == 4 /*&& myQuadType == QUAD_STANDARD*/) // 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) )) @@ -3772,29 +3886,72 @@ 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 ); + + if ( myHelper->HasDegeneratedEdges() && myForcedPnts.empty() ) + { + // "smooth" by computing node positions using 3D TFI and further projection - // TODO: do not smooth fixed nodes + int nbhoriz = quad->iSize; + int nbvertic = quad->jSize; + + 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 ); + + 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(); + } + } + return; + } + + // Get nodes to smooth typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap; TNo2SmooNoMap smooNoMap; - 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 ) ); - - SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); - SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace ); - SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes(); + // fixed nodes + set< const SMDS_MeshNode* > fixedNodes; + for ( size_t i = 0; i < myForcedPnts.size(); ++i ) + { + 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( geomFace, 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 ); @@ -3812,16 +3969,22 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad) } } // set _uv of smooth nodes on FACE boundary - for ( unsigned i = 0; i < quad->side.size(); ++i ) - { - const vector& uvVec = quad->side[i].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 ); - } - } + 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(); @@ -3850,22 +4013,16 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad) { // compute a new XYZ gp_XYZ newXYZ (0,0,0); - for ( unsigned i = 0; i < sNode._triangles.size(); ++i ) + 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 - proj.Perform( newXYZ ); - isValid = ( proj.IsDone() && proj.NbPoints() > 0 ); - if ( isValid ) - { - // check validity of the newUV - Quantity_Parameter u,v; - proj.LowerDistanceParameters( u, v ); - newUV.SetCoord( u, v ); - for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i ) - isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward ); - } + newUV = surface->NextValueOfUV( sNode._uv, newXYZ, 10*tol ).XY(); + + // check validity of the newUV + for ( size_t i = 0; i < sNode._triangles.size() && isValid; ++i ) + isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward ); } if ( !isValid ) { @@ -3883,7 +4040,7 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad) if ( isValid ) { sNode._uv = newUV; - sNode._xyz = surface->Value( newUV.X(), newUV.Y() ).XYZ(); + sNode._xyz = surface->Value( newUV ).XYZ(); } } } @@ -3897,7 +4054,7 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad) continue; // not movable node SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( n2sn->first ); - gp_Pnt xyz = surface->Value( sNode._uv.X(), sNode._uv.Y() ); + gp_Pnt xyz = surface->Value( sNode._uv ); meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() ); // store the new UV @@ -3917,19 +4074,159 @@ 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() ); } } } -/*//================================================================================ +//================================================================================ +/*! + * \brief Checks validity of generated faces + */ +//================================================================================ + +bool StdMeshers_Quadrangle_2D::check() +{ + const bool isOK = true; + if ( !myCheckOri || myQuadList.empty() || !myQuadList.front() || !myHelper ) + return isOK; + + TopoDS_Face geomFace = TopoDS::Face( myHelper->GetSubShape() ); + SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); + SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace ); + bool toCheckUV; + if ( geomFace.Orientation() >= TopAbs_INTERNAL ) geomFace.Orientation( TopAbs_FORWARD ); + + // Get a reference orientation sign + + double okSign; + { + TError err; + TSideVector wireVec = + StdMeshers_FaceSide::GetFaceWires( geomFace, *myHelper->GetMesh(), true, err ); + StdMeshers_FaceSidePtr wire = wireVec[0]; + + // find a right angle VERTEX + int iVertex = 0; + double maxAngle = -1e100; + for ( int i = 0; i < wire->NbEdges(); ++i ) + { + int iPrev = myHelper->WrapIndex( i-1, wire->NbEdges() ); + const TopoDS_Edge& e1 = wire->Edge( iPrev ); + const TopoDS_Edge& e2 = wire->Edge( i ); + double angle = myHelper->GetAngle( e1, e2, geomFace, wire->FirstVertex( i )); + if (( maxAngle < angle ) && + ( 5.* M_PI/180 < angle && angle < 175.* M_PI/180 )) + { + maxAngle = angle; + iVertex = i; + } + } + if ( maxAngle < -2*M_PI ) return isOK; + + // get a sign of 2D area of a corner face + + int iPrev = myHelper->WrapIndex( iVertex-1, wire->NbEdges() ); + const TopoDS_Edge& e1 = wire->Edge( iPrev ); + const TopoDS_Edge& e2 = wire->Edge( iVertex ); + + gp_Vec2d v1, v2; gp_Pnt2d p; + double u[2]; + { + bool rev = ( e1.Orientation() == TopAbs_REVERSED ); + Handle(Geom2d_Curve) c = BRep_Tool::CurveOnSurface( e1, geomFace, u[0], u[1] ); + c->D1( u[ !rev ], p, v1 ); + if ( !rev ) + v1.Reverse(); + } + { + bool rev = ( e2.Orientation() == TopAbs_REVERSED ); + Handle(Geom2d_Curve) c = BRep_Tool::CurveOnSurface( e2, geomFace, u[0], u[1] ); + c->D1( u[ rev ], p, v2 ); + if ( rev ) + v2.Reverse(); + } + + okSign = v2 ^ v1; + + if ( maxAngle < 0 ) + okSign *= -1; + } + + // Look for incorrectly oriented faces + + std::list badFaces; + + const SMDS_MeshNode* nn [ 8 ]; // 8 is just for safety + gp_UV uv [ 8 ]; + SMDS_ElemIteratorPtr fIt = fSubMesh->GetElements(); + while ( fIt->more() ) // loop on faces bound to a FACE + { + const SMDS_MeshElement* f = fIt->next(); + + const int nbN = f->NbCornerNodes(); + for ( int i = 0; i < nbN; ++i ) + nn[ i ] = f->GetNode( i ); + + const SMDS_MeshNode* nInFace = 0; + if ( myHelper->HasSeam() ) + for ( int i = 0; i < nbN && !nInFace; ++i ) + if ( !myHelper->IsSeamShape( nn[i]->getshapeId() )) + nInFace = nn[i]; + + toCheckUV = true; + for ( int i = 0; i < nbN; ++i ) + uv[ i ] = myHelper->GetNodeUV( geomFace, nn[i], nInFace, &toCheckUV ); + + switch ( nbN ) { + case 4: + { + double sign1 = getArea( uv[0], uv[1], uv[2] ); + double sign2 = getArea( uv[0], uv[2], uv[3] ); + if ( sign1 * sign2 < 0 ) + { + sign2 = getArea( uv[1], uv[2], uv[3] ); + sign1 = getArea( uv[1], uv[3], uv[0] ); + if ( sign1 * sign2 < 0 ) + continue; // this should not happen + } + if ( sign1 * okSign < 0 ) + badFaces.push_back ( f ); + break; + } + case 3: + { + double sign = getArea( uv[0], uv[1], uv[2] ); + if ( sign * okSign < 0 ) + badFaces.push_back ( f ); + break; + } + default:; + } + } + + 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 ); + + return !isOK; + } + + return isOK; +} + +//================================================================================ /*! * \brief Finds vertices at the most sharp face corners * \param [in] theFace - the FACE @@ -3954,10 +4251,11 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, theNbDegenEdges = 0; SMESH_MesherHelper helper( theMesh ); + StdMeshers_FaceSide faceSide( theFace, theWire, &theMesh, /*isFwd=*/true, /*skipMedium=*/true); // sort theVertices by angle multimap vertexByAngle; - TopTools_DataMapOfShapeReal angleByVertex; + TopTools_DataMapOfShapeReal angleByVertex; TopoDS_Edge prevE = theWire.back(); if ( SMESH_Algo::isDegenerated( prevE )) { @@ -3969,17 +4267,17 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, prevE = *edge; } list::iterator edge = theWire.begin(); - for ( ; edge != theWire.end(); ++edge ) + for ( int iE = 0; edge != theWire.end(); ++edge, ++iE ) { if ( SMESH_Algo::isDegenerated( *edge )) { ++theNbDegenEdges; continue; } - TopoDS_Vertex v = helper.IthVertex( 0, *edge ); - if ( !theConsiderMesh || SMESH_Algo::VertexNode( v, helper.GetMeshDS() )) + if ( !theConsiderMesh || faceSide.VertexNode( iE )) { - double angle = SMESH_MesherHelper::GetAngle( prevE, *edge, theFace ); + TopoDS_Vertex v = helper.IthVertex( 0, *edge ); + double angle = helper.GetAngle( prevE, *edge, theFace, v ); vertexByAngle.insert( make_pair( angle, v )); angleByVertex.Bind( v, angle ); } @@ -3998,6 +4296,17 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, triaVertex.Nullify(); // check nb of available corners + 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 ( vertexByAngle.size() < 3 ) @@ -4031,18 +4340,22 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, if ( nbCorners == 3 ) vMap.Add( triaVertex ); multimap::reverse_iterator a2v = vertexByAngle.rbegin(); - for ( ; a2v != vertexByAngle.rend() && vMap.Extent() < nbCorners; ++a2v ) + 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 isThereVariants = false; - if ( vertexByAngle.size() > nbCorners ) + bool haveVariants = false; + if ((int) vertexByAngle.size() > nbCorners ) { double lostAngle = a2v->first; double lastAngle = ( --a2v, a2v->first ); - isThereVariants = ( lostAngle * 1.1 >= lastAngle ); + haveVariants = ( lostAngle * 1.1 >= lastAngle ); } + const double angleTol = 5.* M_PI/180; + myCheckOri = ( (int)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() ))) || @@ -4058,9 +4371,10 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, vector< double > angles; vector< TopoDS_Edge > edgeVec; vector< int > cornerInd, nbSeg; - angles.reserve( vertexByAngle.size() ); + int nbSegTot = 0; + angles .reserve( vertexByAngle.size() ); edgeVec.reserve( vertexByAngle.size() ); - nbSeg.reserve( vertexByAngle.size() ); + nbSeg .reserve( vertexByAngle.size() ); cornerInd.reserve( nbCorners ); for ( edge = theWire.begin(); edge != theWire.end(); ++edge ) { @@ -4073,105 +4387,219 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, theVertices.push_back( v ); cornerInd.push_back( angles.size() ); } - angles.push_back( angleByVertex.IsBound( v ) ? angleByVertex( v ) : -M_PI ); + angles .push_back( angleByVertex.IsBound( v ) ? angleByVertex( v ) : -M_PI ); edgeVec.push_back( *edge ); - if ( theConsiderMesh && isThereVariants ) + 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 elual by length if + // refine the result vector - make sides equal by length if // there are several equal angles - if ( isThereVariants ) + if ( haveVariants ) { if ( nbCorners == 3 ) angles[0] = 2 * M_PI; // not to move the base triangle VERTEX - set< int > refinedCorners; + // 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 ) { - int iV = cornerInd[iC]; - if ( !refinedCorners.insert( iV ).second ) + TGeoIndex iV = cornerInd[iC]; + if ( !treatedCorners.insert( iV ).second ) continue; - list< int > equalVertices; - equalVertices.push_back( iV ); + 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() ); - int iVNext = helper.WrapIndex( iV + dV, angles.size() ); + 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] ) < 0.1 * angles[iV]; + bool equal = Abs( angles[iV] - angles[iVNext] ) < angleTol; if ( equal ) - equalVertices.insert( isFwd ? equalVertices.end() : equalVertices.begin(), iVNext ); + 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 ]++; - refinedCorners.insert( cornerInd[ iCNext ] ); + 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; + size_t 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 = equalVertices.size(); + int nbEqualV = equVerts.size(); int nbExcessV = nbEqualV - ( 1 + nbC[0] + nbC[1] ); - if ( nbExcessV > 0 ) + if ( nbExcessV > 0 ) // there is nbExcessV vertices that can become corners { - // calculate normalized length of each side enclosed between neighbor equalVertices - vector< double > curLengths; + // calculate normalized length of each "side" enclosed between neighbor equVerts + vector< double > accuLength; double totalLen = 0; - vector< int > evVec( equalVertices.begin(), equalVertices.end() ); - int iEV = 0; - int iE = cornerInd[ helper.WrapIndex( iC - nbC[0] - 1, cornerInd.size() )]; - int iEEnd = cornerInd[ helper.WrapIndex( iC + nbC[1] + 1, cornerInd.size() )]; - while ( curLengths.size() < nbEqualV + 1 ) + vector< TGeoIndex > evVec( equVerts.begin(), equVerts.end() ); + size_t 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 ((int) accuLength.size() < nbEqualV + int( !allCornersSame ) ) { - curLengths.push_back( totalLen ); + // accumulate length of edges before iEV-th equal vertex + accuLength.push_back( totalLen ); do { - curLengths.back() += SMESH_Algo::EdgeLength( edgeVec[ iE ]); + accuLength.back() += SMESH_Algo::EdgeLength( edgeVec[ iE ]); iE = helper.WrapIndex( iE + 1, edgeVec.size()); - if ( iEV < evVec.size() && iE == evVec[ iEV++ ] ) - break; + if ( iEV < evVec.size() && iE == evVec[ iEV ] ) { + iEV++; + break; // equal vertex reached + } } while( iE != iEEnd ); - totalLen = curLengths.back(); + totalLen = accuLength.back(); } - curLengths.resize( equalVertices.size() ); - for ( size_t iS = 0; iS < curLengths.size(); ++iS ) - curLengths[ iS ] /= totalLen; + accuLength.resize( equVerts.size() ); + for ( size_t iS = 0; iS < accuLength.size(); ++iS ) + accuLength[ iS ] /= totalLen; - // find equalVertices most close to the ideal sub-division of all sides + // 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 = 2 + nbC[0] + nbC[1]; + 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 = 1.; - for ( iEV = iBestEV; iEV < curLengths.size(); ++iEV ) - if (( d = Abs( idealLen - curLengths[ iEV ])) < bestDist ) + 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; } @@ -4183,7 +4611,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) { } @@ -4241,49 +4669,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; } @@ -4343,8 +4788,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 ); @@ -4352,7 +4795,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; @@ -4407,8 +4850,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 ); @@ -4455,6 +4899,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 ); @@ -4491,6 +4938,7 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() << myForcedPnts[iFP].xyz.Y() << ", " << myForcedPnts[iFP].xyz.Z() << " )"); } + myNeedSmooth = true; } // loop on enforced points @@ -4508,25 +4956,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 ) @@ -4553,7 +5007,7 @@ int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J) myQuadList.push_back( FaceQuadStruct::Ptr( newQuad )); vector points; - if ( I > 0 ) + if ( I > 0 && I <= quad->iSize-2 ) { points.reserve( quad->jSize ); for ( int jP = 0; jP < quad->jSize; ++jP ) @@ -4586,13 +5040,15 @@ 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; } - else if ( J > 0 ) //// split horizontally, a new quad is below an old one + else if ( J > 0 && J <= quad->jSize-2 ) //// split horizontally, a new quad is below an old one { points.reserve( quad->iSize ); for ( int iP = 0; iP < quad->iSize; ++iP ) @@ -4621,8 +5077,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; @@ -4666,9 +5124,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 @@ -4678,7 +5136,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() ); @@ -4693,6 +5151,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 ); @@ -5090,6 +5550,7 @@ FaceQuadStruct::Side& FaceQuadStruct::Side::operator=(const Side& otherSide) oSide->contacts[iOC].other_side = this; } } + return *this; } //================================================================================ @@ -5136,6 +5597,7 @@ bool FaceQuadStruct::Side::Reverse(bool keepGrid) grid->Reverse(); } } + return (bool)grid; } //================================================================================ @@ -5170,9 +5632,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();