X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=ada3ab0ccc43d7e524949325cb589e543ffa13dc;hp=5146d08d215d295994281be14f2842b4053a3843;hb=8b6d98aa4acd8176f72ff5b96693bcfef0f9ebd5;hpb=9a54694a0ab1e5cbc558a35c4606ceea4f7af2ef diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index 5146d08d2..ada3ab0cc 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2012 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 @@ -6,7 +6,7 @@ // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either -// version 2.1 of the License. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -30,25 +30,33 @@ #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" +#include "SMESH_HypoFilter.hxx" #include "SMESH_Mesh.hxx" +#include "SMESH_MeshAlgos.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_subMesh.hxx" #include "StdMeshers_FaceSide.hxx" #include "StdMeshers_QuadrangleParams.hxx" #include "StdMeshers_ViscousLayers2D.hxx" +#include +#include #include +#include +#include #include #include #include -#include #include +#include #include #include #include +#include #include #include #include @@ -56,22 +64,19 @@ #include "utilities.h" #include "Utils_ExceptHandlers.hxx" -#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) -#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; + //============================================================================= /*! - * + * */ //============================================================================= @@ -82,10 +87,11 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, myTrianglePreference(false), myTriaVertexID(-1), myNeedSmooth(false), + myCheckOri(false), + myParams( NULL ), myQuadType(QUAD_STANDARD), - myHelper( 0 ) + myHelper( NULL ) { - MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D"); _name = "Quadrangle_2D"; _shapeType = (1 << TopAbs_FACE); _compatibleHypothesis.push_back("QuadrangleParams"); @@ -96,13 +102,12 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, //============================================================================= /*! - * + * */ //============================================================================= StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D() { - MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D"); } //============================================================================= @@ -116,29 +121,30 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis const TopoDS_Shape& aShape, SMESH_Hypothesis::Hypothesis_Status& aStatus) { - bool isOk = true; + myTriaVertexID = -1; + myQuadType = QUAD_STANDARD; + myQuadranglePreference = false; + myTrianglePreference = false; + myHelper = (SMESH_MesherHelper*)NULL; + myParams = NULL; + myQuadList.clear(); + aStatus = SMESH_Hypothesis::HYP_OK; const list & hyps = GetUsedHypothesis(aMesh, aShape, false); const SMESHDS_Hypothesis * aHyp = 0; - myTriaVertexID = -1; - myQuadType = QUAD_STANDARD; - myQuadranglePreference = false; - myTrianglePreference = false; - myQuadStruct.reset(); - bool isFirstParams = true; // First assigned hypothesis (if any) is processed now if (hyps.size() > 0) { aHyp = hyps.front(); - if (strcmp("QuadrangleParams", aHyp->GetName()) == 0) { - const StdMeshers_QuadrangleParams* aHyp1 = - (const StdMeshers_QuadrangleParams*)aHyp; - myTriaVertexID = aHyp1->GetTriaVertex(); - myQuadType = aHyp1->GetQuadType(); + if (strcmp("QuadrangleParams", aHyp->GetName()) == 0) + { + myParams = (const StdMeshers_QuadrangleParams*)aHyp; + myTriaVertexID = myParams->GetTriaVertex(); + myQuadType = myParams->GetQuadType(); if (myQuadType == QUAD_QUADRANGLE_PREF || myQuadType == QUAD_QUADRANGLE_PREF_REVERSED) myQuadranglePreference = true; @@ -151,7 +157,7 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis } else if (strcmp("TrianglePreference", aHyp->GetName()) == 0){ isFirstParams = false; - myTrianglePreference = true; + myTrianglePreference = true; } else { isFirstParams = false; @@ -164,18 +170,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 @@ -189,12 +195,14 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis } } - return isOk; + error( StdMeshers_ViscousLayers2D::CheckHypothesis( aMesh, aShape, aStatus )); + + return aStatus == HYP_OK; } //============================================================================= /*! - * + * */ //============================================================================= @@ -202,89 +210,233 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) { const TopoDS_Face& F = TopoDS::Face(aShape); - Handle(Geom_Surface) S = BRep_Tool::Surface(F); + aMesh.GetSubMesh( F ); - SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); - aMesh.GetSubMesh(aShape); + // do not initialize my fields before this as StdMeshers_ViscousLayers2D + // can call Compute() recursively + SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F ); + if ( !proxyMesh ) + return false; + + myProxyMesh = proxyMesh; SMESH_MesherHelper helper (aMesh); myHelper = &helper; - myProxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F ); - if ( !myProxyMesh ) - return false; - _quadraticMesh = myHelper->IsQuadraticSubMesh(aShape); + myHelper->SetElementsOnShape( true ); myNeedSmooth = false; + myCheckOri = false; - FaceQuadStruct::Ptr quad = CheckNbEdges(aMesh, aShape); + FaceQuadStruct::Ptr quad = CheckNbEdges( aMesh, F, /*considerMesh=*/true, myHelper ); if (!quad) return false; - myQuadStruct = quad; + myQuadList.clear(); + myQuadList.push_back( quad ); - if (myQuadranglePreference) { - int n1 = quad->side[0]->NbPoints(); - int n2 = quad->side[1]->NbPoints(); - int n3 = quad->side[2]->NbPoints(); - int n4 = quad->side[3]->NbPoints(); + if ( !getEnforcedUV() ) + return false; + + updateDegenUV( quad ); + + int n1 = quad->side[0].NbPoints(); + int n2 = quad->side[1].NbPoints(); + int n3 = quad->side[2].NbPoints(); + int n4 = quad->side[3].NbPoints(); + + enum { NOT_COMPUTED = -1, COMPUTE_FAILED = 0, COMPUTE_OK = 1 }; + int res = NOT_COMPUTED; + if ( myQuadranglePreference ) + { int nfull = n1+n2+n3+n4; - int ntmp = nfull/2; - ntmp = ntmp*2; - if (nfull == ntmp && ((n1 != n3) || (n2 != n4))) { - // special path for using only quandrangle faces - bool ok = ComputeQuadPref(aMesh, aShape, quad); - if ( ok && myNeedSmooth ) - Smooth( quad ); - return ok; - } - } - else if (myQuadType == QUAD_REDUCED) { - int n1 = quad->side[0]->NbPoints(); - int n2 = quad->side[1]->NbPoints(); - int n3 = quad->side[2]->NbPoints(); - int n4 = quad->side[3]->NbPoints(); - int n13 = n1 - n3; - int n24 = n2 - n4; + if ((nfull % 2) == 0 && ((n1 != n3) || (n2 != n4))) + { + // special path genarating only quandrangle faces + res = computeQuadPref( aMesh, F, quad ); + } + } + else if ( myQuadType == QUAD_REDUCED ) + { + int n13 = n1 - n3; + int n24 = n2 - n4; int n13tmp = n13/2; n13tmp = n13tmp*2; int n24tmp = n24/2; n24tmp = n24tmp*2; if ((n1 == n3 && n2 != n4 && n24tmp == n24) || - (n2 == n4 && n1 != n3 && n13tmp == n13)) { - bool ok = ComputeReduced(aMesh, aShape, quad); - if ( ok && myNeedSmooth ) - Smooth( quad ); - return ok; + (n2 == n4 && n1 != n3 && n13tmp == n13)) + { + res = computeReduced( aMesh, F, quad ); + } + else + { + if ( n1 != n3 && n2 != n4 ) + error( COMPERR_WARNING, + "To use 'Reduced' transition, " + "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 if ( ! ( n1 == n3 && n2 == n4 )) + error( COMPERR_WARNING, + "To use 'Reduced' transition, " + "two opposite sides should have an even difference in number of segments. " + "'Standard' transion has been used."); } } - // set normalized grid on unit square in parametric domain - - if (!SetNormalizedGrid(aMesh, aShape, quad)) + if ( res == NOT_COMPUTED ) + { + if ( n1 != n3 || n2 != n4 ) + res = computeTriangles( aMesh, F, quad ); + else + res = computeQuadDominant( aMesh, F ); + } + + if ( res == COMPUTE_OK && myNeedSmooth ) + smooth( quad ); + + if ( res == COMPUTE_OK ) + res = check(); + + return ( res == COMPUTE_OK ); +} + +//================================================================================ +/*! + * \brief Compute quadrangles and triangles on the quad + */ +//================================================================================ + +bool StdMeshers_Quadrangle_2D::computeTriangles(SMESH_Mesh& aMesh, + const TopoDS_Face& aFace, + FaceQuadStruct::Ptr quad) +{ + int nb = quad->side[0].grid->NbPoints(); + int nr = quad->side[1].grid->NbPoints(); + int nt = quad->side[2].grid->NbPoints(); + int nl = quad->side[3].grid->NbPoints(); + + // rotate the quad to have nbNodeOut sides on TOP [and LEFT] + if ( nb > nt ) + quad->shift( nl > nr ? 3 : 2, true ); + else if ( nr > nl ) + quad->shift( 1, true ); + else if ( nl > nr ) + quad->shift( nt > nb ? 0 : 3, true ); + + if ( !setNormalizedGrid( quad )) + return false; + + if ( quad->nbNodeOut( QUAD_TOP_SIDE )) + { + splitQuad( quad, 0, quad->jSize-2 ); + } + if ( quad->nbNodeOut( QUAD_BOTTOM_SIDE )) // this should not happen + { + splitQuad( quad, 0, 1 ); + } + 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 ) + botQuad->side[ QUAD_LEFT_SIDE ].to += botQuad->nbNodeOut( QUAD_LEFT_SIDE ); + else if ( botQuad->nbNodeOut( QUAD_RIGHT_SIDE ) > 0 ) + botQuad->side[ QUAD_RIGHT_SIDE ].to += botQuad->nbNodeOut( QUAD_RIGHT_SIDE ); + } + // make quad be a greatest one + if ( quad->side[ QUAD_LEFT_SIDE ].NbPoints() == 2 || + quad->side[ QUAD_RIGHT_SIDE ].NbPoints() == 2 ) + quad = newQuad; + if ( !setNormalizedGrid( quad )) + return false; + } + + if ( quad->nbNodeOut( QUAD_RIGHT_SIDE )) + { + splitQuad( quad, quad->iSize-2, 0 ); + } + 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; + } + } + } + + if ( ! computeQuadDominant( aMesh, aFace )) return false; - // --- compute 3D values on points, store points & quadrangles + // try to fix zero-area triangles near straight-angle corners - int nbdown = quad->side[0]->NbPoints(); - int nbup = quad->side[2]->NbPoints(); + return true; +} - int nbright = quad->side[1]->NbPoints(); - int nbleft = quad->side[3]->NbPoints(); +//================================================================================ +/*! + * \brief Compute quadrangles and possibly triangles on all quads of myQuadList + */ +//================================================================================ - int nbhoriz = Min(nbdown, nbup); - int nbvertic = Min(nbright, nbleft); +bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, + const TopoDS_Face& aFace) +{ + if ( !addEnforcedNodes() ) + return false; + + std::list< FaceQuadStruct::Ptr >::iterator quad = myQuadList.begin(); + for ( ; quad != myQuadList.end(); ++quad ) + if ( !computeQuadDominant( aMesh, aFace, *quad )) + return false; + + return true; +} + +//================================================================================ +/*! + * \brief Compute quadrangles and possibly triangles + */ +//================================================================================ + +bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, + const TopoDS_Face& aFace, + FaceQuadStruct::Ptr quad) +{ + // --- set normalized grid on unit square in parametric domain + + if ( !setNormalizedGrid( quad )) + return false; + + // --- create nodes on points, and create quadrangles + + int nbhoriz = quad->iSize; + int nbvertic = quad->jSize; // internal mesh nodes - int i, j, geomFaceID = meshDS->ShapeToIndex(F); - for (i = 1; i < nbhoriz - 1; i++) { - for (j = 1; j < nbvertic - 1; j++) { - int ij = j * nbhoriz + i; - double u = quad->uv_grid[ij].u; - double v = quad->uv_grid[ij].v; - gp_Pnt P = S->Value(u, v); - SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(node, geomFaceID, u, v); - quad->uv_grid[ij].node = node; + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + Handle(Geom_Surface) S = BRep_Tool::Surface(aFace); + int i,j, geomFaceID = meshDS->ShapeToIndex(aFace); + for (i = 1; i < nbhoriz - 1; i++) + for (j = 1; j < nbvertic - 1; j++) + { + UVPtStruct& uvPnt = quad->UVPt( i, j ); + gp_Pnt P = S->Value( uvPnt.u, uvPnt.v ); + uvPnt.node = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace( uvPnt.node, geomFaceID, uvPnt.u, uvPnt.v ); } - } // mesh faces @@ -300,43 +452,45 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, // i // [0] - i = 0; int ilow = 0; int iup = nbhoriz - 1; - if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; } + if (quad->nbNodeOut(3)) { ilow++; } else { if (quad->nbNodeOut(1)) iup--; } int jlow = 0; int jup = nbvertic - 1; - if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; } + if (quad->nbNodeOut(0)) { jlow++; } else { if (quad->nbNodeOut(2)) jup--; } // regular quadrangles for (i = ilow; i < iup; i++) { for (j = jlow; j < jup; j++) { const SMDS_MeshNode *a, *b, *c, *d; - a = quad->uv_grid[j * nbhoriz + i ].node; - b = quad->uv_grid[j * nbhoriz + i + 1].node; + a = quad->uv_grid[ j * nbhoriz + i ].node; + 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); } } - const vector& uv_e0 = quad->side[0]->GetUVPtStruct(true,0); - const vector& uv_e1 = quad->side[1]->GetUVPtStruct(false,1); - const vector& uv_e2 = quad->side[2]->GetUVPtStruct(true,1); - const vector& uv_e3 = quad->side[3]->GetUVPtStruct(false,0); + // Boundary elements (must always be on an outer boundary of the FACE) + + const vector& uv_e0 = quad->side[0].grid->GetUVPtStruct(); + const vector& uv_e1 = quad->side[1].grid->GetUVPtStruct(); + const vector& uv_e2 = quad->side[2].grid->GetUVPtStruct(); + const vector& uv_e3 = quad->side[3].grid->GetUVPtStruct(); if (uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty()) return error(COMPERR_BAD_INPUT_MESH); double eps = Precision::Confusion(); - // Boundary quadrangles - - if (quad->isEdgeOut[0]) { + int nbdown = (int) uv_e0.size(); + int nbup = (int) uv_e2.size(); + int nbright = (int) uv_e1.size(); + int nbleft = (int) uv_e3.size(); + + if (quad->nbNodeOut(0) && nbvertic == 2) // this should not occur + { // Down edge is out // // |___|___|___|___|___|___| @@ -353,12 +507,16 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, // number of last node of the down edge to be processed int stop = nbdown - 1; // if right edge is out, we will stop at a node, previous to the last one - if (quad->isEdgeOut[1]) stop--; - + //if (quad->nbNodeOut(1)) stop--; + if ( quad->nbNodeOut( QUAD_RIGHT_SIDE )) + quad->UVPt( nbhoriz-1, 1 ).node = uv_e1[1].node; + if ( quad->nbNodeOut( QUAD_LEFT_SIDE )) + quad->UVPt( 0, 1 ).node = uv_e3[1].node; + // 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()); @@ -372,6 +530,7 @@ bool StdMeshers_Quadrangle_2D::Compute (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++) { @@ -394,8 +553,7 @@ bool StdMeshers_Quadrangle_2D::Compute (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) @@ -405,11 +563,10 @@ bool StdMeshers_Quadrangle_2D::Compute (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 { - SplitQuad(meshDS, geomFaceID, a, b, c, d); + splitQuadFace(meshDS, geomFaceID, a, b, c, d); } // if node d is not at position g - make additional triangles @@ -420,15 +577,15 @@ bool StdMeshers_Quadrangle_2D::Compute (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; } } } else { - if (quad->isEdgeOut[2]) { + if (quad->nbNodeOut(2) && nbvertic == 2) + { // Up edge is out // // <-<-<-<-<-<-<-<-<-<-<-<-< -- direction of processing @@ -443,17 +600,65 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, int g = nbhoriz - 1; // last processed node in the regular grid + ilow = 0; + iup = nbhoriz - 1; + int stop = 0; - // if left edge is out, we will stop at a second node - if (quad->isEdgeOut[3]) stop++; + if ( quad->side[3].grid->Edge(0).IsNull() ) // left side is simulated one + { + if ( nbright == 2 ) // quad divided at I but not at J (2D_mesh_QuadranglePreference_01/B1) + 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; @@ -469,7 +674,7 @@ bool StdMeshers_Quadrangle_2D::Compute (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; @@ -482,8 +687,7 @@ bool StdMeshers_Quadrangle_2D::Compute (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) @@ -492,22 +696,20 @@ bool StdMeshers_Quadrangle_2D::Compute (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 { - SplitQuad(meshDS, geomFaceID, a, b, c, d); + splitQuadFace(meshDS, geomFaceID, a, b, c, d); } - if (near + 1 < g) { // if d not is at g - make additional triangles + if (near + 1 < g) { // if d is not at g - make additional triangles for (int k = near + 1; k < g; k++) { c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node; if (k + 1 > iup) 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; @@ -517,20 +719,23 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } // right or left boundary quadrangles - if (quad->isEdgeOut[1]) { -// MESSAGE("right edge is out"); + 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; - if (quad->isEdgeOut[2]) stop--; - for (i = 0; i < stop; i++) { + i = 0; + if (quad->side[ QUAD_RIGHT_SIDE ].from != i ) i++; + if (quad->side[ QUAD_RIGHT_SIDE ].to != stop ) stop--; + for ( ; i < stop; i++) { const SMDS_MeshNode *a, *b, *c, *d; a = uv_e1[i].node; b = uv_e1[i + 1].node; 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 { @@ -554,8 +759,7 @@ bool StdMeshers_Quadrangle_2D::Compute (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) @@ -565,11 +769,10 @@ bool StdMeshers_Quadrangle_2D::Compute (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 { - SplitQuad(meshDS, geomFaceID, a, b, c, d); + splitQuadFace(meshDS, geomFaceID, a, b, c, d); } if (near - 1 > g) { // if d not is at g - make additional triangles @@ -579,36 +782,69 @@ bool StdMeshers_Quadrangle_2D::Compute (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->isEdgeOut[3]) { -// MESSAGE("left edge is out"); + if (quad->nbNodeOut(3) && nbhoriz == 2) + { int g = nbvertic - 1; // last processed node in the grid int stop = 0; - if (quad->isEdgeOut[0]) stop++; - for (i = nbleft - 1; i > stop; i--) { - const SMDS_MeshNode *a, *b, *c, *d; + i = quad->side[ QUAD_LEFT_SIDE ].to-1; // nbleft - 1; + + 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()); // 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; - } 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()); @@ -624,32 +860,28 @@ bool StdMeshers_Quadrangle_2D::Compute (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 { - SplitQuad(meshDS, geomFaceID, a, b, c, d); + splitQuadFace(meshDS, geomFaceID, a, b, c, d); } if (near + 1 < g) { // if d not is at g - make additional triangles 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; @@ -658,9 +890,6 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } } - if ( myNeedSmooth ) - Smooth( quad ); - bool isOk = true; return isOk; } @@ -672,19 +901,19 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, */ //============================================================================= -bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape, - MapShapeNbElems& aResMap) +bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh, + const TopoDS_Shape& aFace, + MapShapeNbElems& aResMap) { - aMesh.GetSubMesh(aShape); + aMesh.GetSubMesh(aFace); std::vector aNbNodes(4); bool IsQuadratic = false; - if (!CheckNbEdgesForEvaluate(aMesh, aShape, aResMap, aNbNodes, IsQuadratic)) { + if (!checkNbEdgesForEvaluate(aMesh, aFace, aResMap, aNbNodes, IsQuadratic)) { std::vector aResVec(SMDSEntity_Last); for (int i=SMDSEntity_Node; iGetComputeError(); smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this)); @@ -701,7 +930,7 @@ bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh, ntmp = ntmp*2; if (nfull==ntmp && ((n1!=n3) || (n2!=n4))) { // special path for using only quandrangle faces - return EvaluateQuadPref(aMesh, aShape, aNbNodes, aResMap, IsQuadratic); + return evaluateQuadPref(aMesh, aFace, aNbNodes, aResMap, IsQuadratic); //return true; } } @@ -753,217 +982,615 @@ bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh, aVec[SMDSEntity_Quadrangle] = nbFaces4 - aNbNodes[3] + 1; } } - SMESH_subMesh * sm = aMesh.GetSubMesh(aShape); + SMESH_subMesh * sm = aMesh.GetSubMesh(aFace); aResMap.insert(std::make_pair(sm,aVec)); return true; } - //================================================================================ /*! - * \brief Return true if only two given edges meat at their common vertex + * \brief Return true if the algorithm can mesh this shape + * \param [in] aShape - shape to check + * \param [in] toCheckAll - if true, this check returns OK if all shapes are OK, + * else, returns OK if at least one shape is OK */ //================================================================================ -static bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1, - const TopoDS_Edge& e2, - SMESH_Mesh & mesh) +bool StdMeshers_Quadrangle_2D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll ) { - 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; -} + int nbFoundFaces = 0; + for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces ) + { + 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, totalNbEdges = 0; + TopExp_Explorer eExp( aFace, TopAbs_EDGE ); + for ( ; eExp.More() && nbNoDegenEdges < 3; eExp.Next(), ++totalNbEdges ) { + if ( !SMESH_Algo::isDegenerated( TopoDS::Edge( eExp.Current() ))) + ++nbNoDegenEdges; + } + if ( toCheckAll && ( totalNbEdges < 4 && nbNoDegenEdges < 3 )) return false; + if ( !toCheckAll && ( totalNbEdges >= 4 || nbNoDegenEdges >= 3 )) return true; + } + return ( toCheckAll && nbFoundFaces != 0 ); +} -FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMesh, - const TopoDS_Shape & aShape) +namespace { - if ( myQuadStruct && myQuadStruct->face.IsSame( aShape )) - return myQuadStruct; - - TopoDS_Face F = TopoDS::Face(aShape); - if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD ); - const bool ignoreMediumNodes = _quadraticMesh; + //================================================================================ + /*! + * \brief Return true if only two given edges meat at their common vertex + */ + //================================================================================ - // verify 1 wire only, with 4 edges - list< TopoDS_Edge > edges; - list< int > nbEdgesInWire; - int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire); - if (nbWire != 1) { - error(COMPERR_BAD_SHAPE, TComm("Wrong number of wires: ") << nbWire); - return FaceQuadStruct::Ptr(); + 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; } - FaceQuadStruct::Ptr quad( new FaceQuadStruct ); - quad->uv_grid = 0; - quad->side.reserve(nbEdgesInWire.front()); - quad->face = F; - int nbSides = 0; - list< TopoDS_Edge >::iterator edgeIt = edges.begin(); - if (nbEdgesInWire.front() == 3) // exactly 3 edges + //-------------------------------------------------------------------------------- + /*! + * \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 { - SMESH_Comment comment; - SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); - if (myTriaVertexID < 1) + 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 ) { - comment << "No Base vertex parameter provided for a trilateral geometrical face"; - } - else + 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 { - TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID)); - if (!V.IsNull()) { - TopoDS_Edge E1,E2,E3; - for (; edgeIt != edges.end(); ++edgeIt) { - TopoDS_Edge E = *edgeIt; - TopoDS_Vertex VF, VL; - TopExp::Vertices(E, VF, VL, true); - if (VF.IsSame(V)) - E1 = E; - else if (VL.IsSame(V)) - E3 = E; - else - E2 = E; - } - if (!E1.IsNull() && !E2.IsNull() && !E3.IsNull()) - { - quad->side.push_back(new StdMeshers_FaceSide(F, E1, &aMesh, true, - ignoreMediumNodes, myProxyMesh)); - quad->side.push_back(new StdMeshers_FaceSide(F, E2, &aMesh, true, - ignoreMediumNodes, myProxyMesh)); - quad->side.push_back(new StdMeshers_FaceSide(F, E3, &aMesh, false, - ignoreMediumNodes, myProxyMesh)); - const vector& UVPSleft = quad->side[0]->GetUVPtStruct(true,0); - /* vector& UVPStop = */quad->side[1]->GetUVPtStruct(false,1); - /* vector& UVPSright = */quad->side[2]->GetUVPtStruct(true,1); - const SMDS_MeshNode* aNode = UVPSleft[0].node; - gp_Pnt2d aPnt2d(UVPSleft[0].u, UVPSleft[0].v); - quad->side.push_back(new StdMeshers_FaceSide(aNode, aPnt2d, quad->side[1])); - return quad; - } - } - comment << "Invalid Base vertex parameter: " << myTriaVertexID << " is not among ["; - TopTools_MapOfShape vMap; - for (TopExp_Explorer v(aShape, TopAbs_VERTEX); v.More(); v.Next()) - if (vMap.Add(v.Current())) - comment << meshDS->ShapeToIndex(v.Current()) << (vMap.Extent()==3 ? "]" : ", "); - } - error(comment); - quad.reset(); - return quad; - } - else if (nbEdgesInWire.front() == 4) // exactly 4 edges + 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) { - for (; edgeIt != edges.end(); ++edgeIt, nbSides++) - quad->side.push_back(new StdMeshers_FaceSide(F, *edgeIt, &aMesh, nbSides < QUAD_TOP_SIDE, - ignoreMediumNodes, myProxyMesh)); - } - else if (nbEdgesInWire.front() > 4) // more than 4 edges - try to unite some - { - list< TopoDS_Edge > sideEdges; - vector< int > degenSides; - while (!edges.empty()) { - sideEdges.clear(); - sideEdges.splice(sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end() - bool sameSide = true; - while (!edges.empty() && sameSide) { - sameSide = SMESH_Algo::IsContinuous(sideEdges.back(), edges.front()); - if (sameSide) - sideEdges.splice(sideEdges.end(), edges, edges.begin()); - } - if (nbSides == 0) { // go backward from the first edge - sameSide = true; - while (!edges.empty() && sameSide) { - sameSide = SMESH_Algo::IsContinuous(sideEdges.front(), edges.back()); - if (sameSide) - sideEdges.splice(sideEdges.begin(), edges, --edges.end()); - } + // 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; } - if ( sideEdges.size() == 1 && BRep_Tool::Degenerated( sideEdges.front() )) - degenSides.push_back( nbSides ); - quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, nbSides < QUAD_TOP_SIDE, - ignoreMediumNodes, myProxyMesh)); - ++nbSides; + // 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 ); } - if ( !degenSides.empty() && nbSides - degenSides.size() == 4 ) + + theHaveConcaveVertices = ( nbConvexAngles < nbEdges ); + + if ((int) theVertices.size() == theNbCorners ) + return; + + theVertices.clear(); + + if ( !theConsiderMesh || theNbCorners < 4 || + nbConvexAngles <= theNbCorners || + nbSharpAngles == theNbCorners ) { - myNeedSmooth = true; - for ( unsigned i = QUAD_TOP_SIDE; i < quad->side.size(); ++i ) - quad->side[i]->Reverse(); + if ( nbEdges == theNbCorners ) // return all vertices + { + for ( e = edge0; (int) theVertices.size() < theNbCorners; e = e->myNext ) + theVertices.push_back( e->my1stVertex ); + return; + } - for ( int i = degenSides.size()-1; i > -1; --i ) + // 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 { - StdMeshers_FaceSide* degenSide = quad->side[ degenSides[ i ]]; - delete degenSide; - quad->side.erase( quad->side.begin() + degenSides[ i ] ); + e->myNbSegments += nodes.size() - 1; + totNbSeg += nodes.size() - 1; } - for ( unsigned i = QUAD_TOP_SIDE; i < quad->side.size(); ++i ) - quad->side[i]->Reverse(); - nbSides -= degenSides.size(); + // 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; + } } - // issue 20222. Try to unite only edges shared by two same faces - if (nbSides < 4) + + 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.reset( new FaceQuadStruct ); - quad->side.reserve(nbEdgesInWire.front()); - nbSides = 0; + quad.myCornerE[ 0 ] = edge0; - SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire); - while (!edges.empty()) { - sideEdges.clear(); - sideEdges.splice(sideEdges.end(), edges, edges.begin()); - bool sameSide = true; - while (!edges.empty() && sameSide) { - sameSide = - SMESH_Algo::IsContinuous(sideEdges.back(), edges.front()) && - twoEdgesMeatAtVertex(sideEdges.back(), edges.front(), aMesh); - if (sameSide) - sideEdges.splice(sideEdges.end(), edges, edges.begin()); + // 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 (nbSides == 0) { // go backward from the first edge - sameSide = true; - while (!edges.empty() && sameSide) { - sameSide = - SMESH_Algo::IsContinuous(sideEdges.front(), edges.back()) && - twoEdgesMeatAtVertex(sideEdges.front(), edges.back(), aMesh); - if (sameSide) - sideEdges.splice(sideEdges.begin(), edges, --edges.end()); + 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; } } - quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, - nbSides < QUAD_TOP_SIDE, - ignoreMediumNodes, myProxyMesh)); - ++nbSides; } } + + 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; } - if (nbSides != 4) { -#ifdef _DEBUG_ - MESSAGE ("StdMeshers_Quadrangle_2D. Edge IDs of " << nbSides << " sides:\n"); - for (int i = 0; i < nbSides; ++i) { - MESSAGE (" ("); - for (int e = 0; e < quad->side[i]->NbEdges(); ++e) - MESSAGE (myHelper->GetMeshDS()->ShapeToIndex(quad->side[i]->Edge(e)) << " "); - MESSAGE (")\n"); + +} // namespace + +//================================================================================ +/*! + * \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 ); + 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, + SMESH_MesherHelper* aFaceHelper) +{ + if ( !myQuadList.empty() && myQuadList.front()->face.IsSame( aShape )) + return myQuadList.front(); + + TopoDS_Face F = TopoDS::Face(aShape); + if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD ); + const bool ignoreMediumNodes = _quadraticMesh; + + // verify 1 wire only + list< TopoDS_Edge > edges; + list< int > nbEdgesInWire; + int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire); + if (nbWire != 1) { + error(COMPERR_BAD_SHAPE, TComm("Wrong number of wires: ") << nbWire); + return FaceQuadStruct::Ptr(); + } + + // 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 ) + { + return FaceQuadStruct::Ptr(); + } + FaceQuadStruct::Ptr quad( new FaceQuadStruct ); + quad->side.reserve(nbEdgesInWire.front()); + quad->face = F; + + list< TopoDS_Edge >::iterator edgeIt = edges.begin(); + if ( nbSides == 3 ) // 3 sides and corners[0] is a vertex with myTriaVertexID + { + for ( int iSide = 0; iSide < 3; ++iSide ) + { + list< TopoDS_Edge > sideEdges; + TopoDS_Vertex nextSideV = corners[( iSide + 1 ) % 3 ]; + while ( edgeIt != edges.end() && + !nextSideV.IsSame( SMESH_MesherHelper::IthVertex( 0, *edgeIt ))) + if ( SMESH_Algo::isDegenerated( *edgeIt )) + ++edgeIt; + else + sideEdges.push_back( *edgeIt++ ); + if ( !sideEdges.empty() ) + quad->side.push_back( StdMeshers_FaceSide::New(F, sideEdges, &aMesh, iSide < QUAD_TOP_SIDE, + ignoreMediumNodes, myHelper, myProxyMesh)); + else + --iSide; + } + const vector& UVPSleft = quad->side[0].GetUVPtStruct(true,0); + /* vector& UVPStop = */quad->side[1].GetUVPtStruct(false,1); + /* vector& UVPSright = */quad->side[2].GetUVPtStruct(true,1); + const SMDS_MeshNode* aNode = UVPSleft[0].node; + gp_Pnt2d aPnt2d = UVPSleft[0].UV(); + quad->side.push_back( StdMeshers_FaceSide::New( quad->side[1].grid.get(), aNode, &aPnt2d )); + myNeedSmooth = ( nbDegenEdges > 0 ); + return quad; + } + else // 4 sides + { + myNeedSmooth = ( corners.size() == 4 && nbDegenEdges > 0 ); + int iSide = 0, nbUsedDegen = 0, nbLoops = 0; + for ( ; edgeIt != edges.end(); ++nbLoops ) + { + list< TopoDS_Edge > sideEdges; + TopoDS_Vertex nextSideV = corners[( iSide + 1 - nbUsedDegen ) % corners.size() ]; + bool nextSideVReached = false; + do + { + const TopoDS_Edge& edge = *edgeIt; + nextSideVReached = nextSideV.IsSame( myHelper->IthVertex( 1, edge )); + if ( SMESH_Algo::isDegenerated( edge )) + { + if ( !myNeedSmooth ) // need to make a side on a degen edge + { + if ( sideEdges.empty() ) + { + sideEdges.push_back( edge ); + ++nbUsedDegen; + nextSideVReached = true; + } + else + { + break; + } + } + } + else //if ( !myHelper || !myHelper->IsRealSeam( edge )) + { + 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, myHelper, myProxyMesh )); + ++iSide; + } + if ( quad->side.size() == 4 ) + break; + if ( nbLoops > 8 ) + { + error(TComm("Bug: infinite loop in StdMeshers_Quadrangle_2D::CheckNbEdges()")); + quad.reset(); + break; + } + } + if ( quad && quad->side.size() != 4 ) + { + error(TComm("Bug: ") << quad->side.size() << " sides found instead of 4"); + quad.reset(); } -#endif - if (!nbSides) - nbSides = nbEdgesInWire.front(); - error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides); - quad.reset(); } return quad; @@ -976,11 +1603,11 @@ FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & */ //============================================================================= -bool StdMeshers_Quadrangle_2D::CheckNbEdgesForEvaluate(SMESH_Mesh& aMesh, +bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh& aMesh, const TopoDS_Shape & aShape, - MapShapeNbElems& aResMap, - std::vector& aNbNodes, - bool& IsQuadratic) + MapShapeNbElems& aResMap, + std::vector& aNbNodes, + bool& IsQuadratic) { const TopoDS_Face & F = TopoDS::Face(aShape); @@ -1172,36 +1799,12 @@ StdMeshers_Quadrangle_2D::CheckAnd2Dcompute (SMESH_Mesh & aMesh, if ( quad ) { // set normalized grid on unit square in parametric domain - if (!SetNormalizedGrid(aMesh, aShape, quad)) + if ( ! setNormalizedGrid( quad )) quad.reset(); } return quad; } -//============================================================================= -/*! - * - */ -//============================================================================= - -faceQuadStruct::~faceQuadStruct() -{ - for (size_t i = 0; i < side.size(); i++) { - if (side[i]) { - delete side[i]; - for (size_t j = i+1; j < side.size(); j++) - if ( side[i] == side[j] ) - side[j] = 0; - } - } - side.clear(); - - if (uv_grid) { - delete [] uv_grid; - uv_grid = 0; - } -} - namespace { inline const vector& getUVPtStructIn(FaceQuadStruct::Ptr& quad, int i, int nbSeg) @@ -1209,9 +1812,9 @@ namespace bool isXConst = (i == QUAD_BOTTOM_SIDE || i == QUAD_TOP_SIDE); double constValue = (i == QUAD_BOTTOM_SIDE || i == QUAD_LEFT_SIDE) ? 0 : 1; return - quad->isEdgeOut[i] ? - quad->side[i]->SimulateUVPtStruct(nbSeg,isXConst,constValue) : - quad->side[i]->GetUVPtStruct(isXConst,constValue); + quad->nbNodeOut(i) ? + quad->side[i].grid->SimulateUVPtStruct(nbSeg,isXConst,constValue) : + quad->side[i].grid->GetUVPtStruct (isXConst,constValue); } inline gp_UV calcUV(double x, double y, const gp_UV& a0,const gp_UV& a1,const gp_UV& a2,const gp_UV& a3, @@ -1229,20 +1832,16 @@ namespace */ //============================================================================= -bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, - const TopoDS_Shape& aShape, - FaceQuadStruct::Ptr & quad) +bool StdMeshers_Quadrangle_2D::setNormalizedGrid (FaceQuadStruct::Ptr quad) { + if ( !quad->uv_grid.empty() ) + return true; + // Algorithme décrit dans "Génération automatique de maillages" // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85 // traitement dans le domaine paramétrique 2d u,v // transport - projection sur le carré unité -// MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid"); -// const TopoDS_Face& F = TopoDS::Face(aShape); - - // 1 --- find orientation of the 4 edges, by test on extrema - // max min 0 x1 1 // |<----north-2-------^ a3 -------------> a2 // | | ^1 1^ @@ -1254,87 +1853,135 @@ bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, // min max 0 x0 1 // =down // + const FaceQuadStruct::Side & bSide = quad->side[0]; + const FaceQuadStruct::Side & rSide = quad->side[1]; + const FaceQuadStruct::Side & tSide = quad->side[2]; + const FaceQuadStruct::Side & lSide = quad->side[3]; - // 3 --- 2D normalized values on unit square [0..1][0..1] - - int nbhoriz = Min(quad->side[0]->NbPoints(), quad->side[2]->NbPoints()); - int nbvertic = Min(quad->side[1]->NbPoints(), quad->side[3]->NbPoints()); - - quad->isEdgeOut[0] = (quad->side[0]->NbPoints() > quad->side[2]->NbPoints()); - quad->isEdgeOut[1] = (quad->side[1]->NbPoints() > quad->side[3]->NbPoints()); - quad->isEdgeOut[2] = (quad->side[2]->NbPoints() > quad->side[0]->NbPoints()); - quad->isEdgeOut[3] = (quad->side[3]->NbPoints() > quad->side[1]->NbPoints()); - - UVPtStruct *uv_grid = quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz]; - - const vector& uv_e0 = getUVPtStructIn(quad, 0, nbhoriz - 1); - const vector& uv_e1 = getUVPtStructIn(quad, 1, nbvertic - 1); - const vector& uv_e2 = getUVPtStructIn(quad, 2, nbhoriz - 1); - const vector& uv_e3 = getUVPtStructIn(quad, 3, nbvertic - 1); + 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 ) + { + // all sub-quads must have NO sides with nbNodeOut > 0 + quad->nbNodeOut(0) = Max( 0, bSide.grid->NbPoints() - tSide.grid->NbPoints() ); + quad->nbNodeOut(1) = Max( 0, rSide.grid->NbPoints() - lSide.grid->NbPoints() ); + quad->nbNodeOut(2) = Max( 0, tSide.grid->NbPoints() - bSide.grid->NbPoints() ); + quad->nbNodeOut(3) = Max( 0, lSide.grid->NbPoints() - rSide.grid->NbPoints() ); + } + const vector& uv_e0 = bSide.GetUVPtStruct(); + const vector& uv_e1 = rSide.GetUVPtStruct(); + const vector& uv_e2 = tSide.GetUVPtStruct(); + const vector& uv_e3 = lSide.GetUVPtStruct(); if (uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty()) //return error("Can't find nodes on sides"); return error(COMPERR_BAD_INPUT_MESH); - if ( myNeedSmooth ) - UpdateDegenUV( quad ); + quad->uv_grid.resize( nbvertic * nbhoriz ); + quad->iSize = nbhoriz; + quad->jSize = nbvertic; + UVPtStruct *uv_grid = & quad->uv_grid[0]; + + quad->uv_box.Clear(); // copy data of face boundary - /*if (! quad->isEdgeOut[0])*/ { - const int j = 0; - for (int i = 0; i < nbhoriz; i++) // down - uv_grid[ j * nbhoriz + i ] = uv_e0[i]; + + FaceQuadStruct::SideIterator sideIter; + + { // BOTTOM + const int j = 0; + const double x0 = bSide.First().normParam; + const double dx = bSide.Last().normParam - bSide.First().normParam; + for ( sideIter.Init( bSide ); sideIter.More(); sideIter.Next() ) { + sideIter.UVPt().x = ( sideIter.UVPt().normParam - x0 ) / dx; + sideIter.UVPt().y = 0.; + uv_grid[ j * nbhoriz + sideIter.Count() ] = sideIter.UVPt(); + quad->uv_box.Add( sideIter.UVPt().UV() ); + } } - /*if (! quad->isEdgeOut[1])*/ { - const int i = nbhoriz - 1; - for (int j = 0; j < nbvertic; j++) // right - uv_grid[ j * nbhoriz + i ] = uv_e1[j]; + { // RIGHT + const int i = nbhoriz - 1; + const double y0 = rSide.First().normParam; + const double dy = rSide.Last().normParam - rSide.First().normParam; + sideIter.Init( rSide ); + if ( quad->UVPt( i, sideIter.Count() ).node ) + sideIter.Next(); // avoid copying from a split emulated side + for ( ; sideIter.More(); sideIter.Next() ) { + sideIter.UVPt().x = 1.; + sideIter.UVPt().y = ( sideIter.UVPt().normParam - y0 ) / dy; + uv_grid[ sideIter.Count() * nbhoriz + i ] = sideIter.UVPt(); + quad->uv_box.Add( sideIter.UVPt().UV() ); + } } - /*if (! quad->isEdgeOut[2])*/ { - const int j = nbvertic - 1; - for (int i = 0; i < nbhoriz; i++) // up - uv_grid[ j * nbhoriz + i ] = uv_e2[i]; + { // TOP + const int j = nbvertic - 1; + const double x0 = tSide.First().normParam; + const double dx = tSide.Last().normParam - tSide.First().normParam; + int i = 0, nb = nbhoriz; + sideIter.Init( tSide ); + if ( quad->UVPt( nb-1, j ).node ) --nb; // avoid copying from a split emulated side + for ( ; i < nb; i++, sideIter.Next()) { + sideIter.UVPt().x = ( sideIter.UVPt().normParam - x0 ) / dx; + sideIter.UVPt().y = 1.; + uv_grid[ j * nbhoriz + i ] = sideIter.UVPt(); + quad->uv_box.Add( sideIter.UVPt().UV() ); + } } - /*if (! quad->isEdgeOut[3])*/ { - int i = 0; - for (int j = 0; j < nbvertic; j++) // left - uv_grid[ j * nbhoriz + i ] = uv_e3[j]; + { // LEFT + const int i = 0; + const double y0 = lSide.First().normParam; + const double dy = lSide.Last().normParam - lSide.First().normParam; + int j = 0, nb = nbvertic; + sideIter.Init( lSide ); + if ( quad->UVPt( i, j ).node ) + ++j, sideIter.Next(); // avoid copying from a split emulated side + if ( quad->UVPt( i, nb-1 ).node ) + --nb; + for ( ; j < nb; j++, sideIter.Next()) { + sideIter.UVPt().x = 0.; + sideIter.UVPt().y = ( sideIter.UVPt().normParam - y0 ) / dy; + uv_grid[ j * nbhoriz + i ] = sideIter.UVPt(); + quad->uv_box.Add( sideIter.UVPt().UV() ); + } } // normalized 2d parameters on grid - for (int i = 0; i < nbhoriz; i++) { - for (int j = 0; j < nbvertic; j++) { - int ij = j * nbhoriz + i; - // --- droite i cste : x = x0 + y(x1-x0) - double x0 = uv_e0[i].normParam; // bas - sud - double x1 = uv_e2[i].normParam; // haut - nord - // --- droite j cste : y = y0 + x(y1-y0) - double y0 = uv_e3[j].normParam; // gauche-ouest - double y1 = uv_e1[j].normParam; // droite - est + + for (int i = 1; i < nbhoriz-1; i++) + { + const double x0 = quad->UVPt( i, 0 ).x; + const double x1 = quad->UVPt( i, nbvertic-1 ).x; + for (int j = 1; j < nbvertic-1; j++) + { + const double y0 = quad->UVPt( 0, j ).y; + const double y1 = quad->UVPt( nbhoriz-1, j ).y; // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0) double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); double y = y0 + x * (y1 - y0); + int ij = j * nbhoriz + i; uv_grid[ij].x = x; uv_grid[ij].y = y; - //MESSAGE("-xy-01 "<UVPt( 0, 0 ).UV(); + gp_UV a1 = quad->UVPt( nbhoriz-1, 0 ).UV(); + gp_UV a2 = quad->UVPt( nbhoriz-1, nbvertic-1 ).UV(); + gp_UV a3 = quad->UVPt( 0, nbvertic-1 ).UV(); + + for (int i = 1; i < nbhoriz-1; i++) { - gp_UV p0( uv_e0[i].u, uv_e0[i].v ); - gp_UV p2( uv_e2[i].u, uv_e2[i].v ); - for (int j = 0; j < nbvertic; j++) + gp_UV p0 = quad->UVPt( i, 0 ).UV(); + gp_UV p2 = quad->UVPt( i, nbvertic-1 ).UV(); + for (int j = 1; j < nbvertic-1; j++) { - gp_UV p1( uv_e1[j].u, uv_e1[j].v ); - gp_UV p3( uv_e3[j].u, uv_e3[j].v ); + gp_UV p1 = quad->UVPt( nbhoriz-1, j ).UV(); + gp_UV p3 = quad->UVPt( 0, j ).UV(); int ij = j * nbhoriz + i; double x = uv_grid[ij].x; @@ -1351,39 +1998,91 @@ bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, //======================================================================= //function : ShiftQuad -//purpose : auxilary function for ComputeQuadPref +//purpose : auxiliary function for computeQuadPref //======================================================================= -static void shiftQuad(FaceQuadStruct::Ptr& quad, const int num, bool) +void StdMeshers_Quadrangle_2D::shiftQuad(FaceQuadStruct::Ptr& quad, const int num ) { - quad->shift( num, /*ori=*/true ); + quad->shift( num, /*ori=*/true, /*keepGrid=*/myQuadList.size() > 1 ); } //================================================================================ /*! - * \brief Rotate sides of a quad by nb - * \param nb - number of rotation 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 + * are altered instead */ //================================================================================ -void FaceQuadStruct::shift( size_t nb, bool ori ) +void FaceQuadStruct::shift( size_t nb, bool ori, bool keepGrid ) { if ( nb == 0 ) return; - StdMeshers_FaceSide* sideArr[4] = { side[0], side[1], side[2], side[3] }; - for (int i = QUAD_BOTTOM_SIDE; i < NB_QUAD_SIDES; ++i) { + + 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) + { int id = (i + nb) % NB_QUAD_SIDES; - bool wasForward = (i < QUAD_TOP_SIDE); - bool newForward = (id < QUAD_TOP_SIDE); - if (ori && wasForward != newForward) - sideArr[ i ]->Reverse(); - side[ id ] = sideArr[ i ]; + if ( ori ) + { + bool wasForward = (i < QUAD_TOP_SIDE); + bool newForward = (id < QUAD_TOP_SIDE); + if ( wasForward != newForward ) + side[ i ].Reverse( keepGrid ); + } + newSides[ id ] = side[ i ]; + sidePtrs[ i ] = & side[ i ]; + } + // make newSides refer newSides via Side::Contact's + for ( size_t i = 0; i < newSides.size(); ++i ) + { + FaceQuadStruct::Side& ns = newSides[ i ]; + for ( size_t iC = 0; iC < ns.contacts.size(); ++iC ) + { + FaceQuadStruct::Side* oSide = ns.contacts[iC].other_side; + vector< Side* >::iterator sIt = std::find( sidePtrs.begin(), sidePtrs.end(), oSide ); + if ( sIt != sidePtrs.end() ) + ns.contacts[iC].other_side = & newSides[ *sIt - sidePtrs[0] ]; + } + } + newSides.swap( side ); + + 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, @@ -1394,10 +2093,10 @@ static gp_UV calcUV(double x0, double x1, double y0, double y1, double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); double y = y0 + x * (y1 - y0); - gp_UV p0 = quad->side[QUAD_BOTTOM_SIDE]->Value2d(x).XY(); - gp_UV p1 = quad->side[QUAD_RIGHT_SIDE ]->Value2d(y).XY(); - gp_UV p2 = quad->side[QUAD_TOP_SIDE ]->Value2d(x).XY(); - gp_UV p3 = quad->side[QUAD_LEFT_SIDE ]->Value2d(y).XY(); + gp_UV p0 = quad->side[QUAD_BOTTOM_SIDE].grid->Value2d(x).XY(); + gp_UV p1 = quad->side[QUAD_RIGHT_SIDE ].grid->Value2d(y).XY(); + gp_UV p2 = quad->side[QUAD_TOP_SIDE ].grid->Value2d(x).XY(); + gp_UV p3 = quad->side[QUAD_LEFT_SIDE ].grid->Value2d(y).XY(); gp_UV uv = calcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3); @@ -1406,7 +2105,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, @@ -1414,10 +2113,10 @@ static gp_UV calcUV2(double x, double y, const gp_UV& a0, const gp_UV& a1, const gp_UV& a2, const gp_UV& a3) { - gp_UV p0 = quad->side[QUAD_BOTTOM_SIDE]->Value2d(x).XY(); - gp_UV p1 = quad->side[QUAD_RIGHT_SIDE ]->Value2d(y).XY(); - gp_UV p2 = quad->side[QUAD_TOP_SIDE ]->Value2d(x).XY(); - gp_UV p3 = quad->side[QUAD_LEFT_SIDE ]->Value2d(y).XY(); + gp_UV p0 = quad->side[QUAD_BOTTOM_SIDE].grid->Value2d(x).XY(); + gp_UV p1 = quad->side[QUAD_RIGHT_SIDE ].grid->Value2d(y).XY(); + gp_UV p2 = quad->side[QUAD_TOP_SIDE ].grid->Value2d(x).XY(); + gp_UV p3 = quad->side[QUAD_LEFT_SIDE ].grid->Value2d(y).XY(); gp_UV uv = calcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3); @@ -1431,70 +2130,63 @@ static gp_UV calcUV2(double x, double y, */ //======================================================================= -bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, - const TopoDS_Shape& aShape, +bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, + const TopoDS_Face& aFace, FaceQuadStruct::Ptr quad) { - // Auxilary key in order to keep old variant - // of meshing after implementation new variant - // for bug 0016220 from Mantis. - bool OldVersion = false; - if (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED) - OldVersion = true; + const bool OldVersion = (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED); + const bool WisF = true; - SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); - const TopoDS_Face& F = TopoDS::Face(aShape); - Handle(Geom_Surface) S = BRep_Tool::Surface(F); - bool WisF = true; - int i,j,geomFaceID = meshDS->ShapeToIndex(F); + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + Handle(Geom_Surface) S = BRep_Tool::Surface(aFace); + int i,j, geomFaceID = meshDS->ShapeToIndex(aFace); - int nb = quad->side[0]->NbPoints(); - int nr = quad->side[1]->NbPoints(); - int nt = quad->side[2]->NbPoints(); - int nl = quad->side[3]->NbPoints(); + int nb = quad->side[0].NbPoints(); + int nr = quad->side[1].NbPoints(); + int nt = quad->side[2].NbPoints(); + int nl = quad->side[3].NbPoints(); int dh = abs(nb-nt); int dv = abs(nr-nl); - if (dh>=dv) { - if (nt>nb) { - // it is a base case => not shift quad but me be replacement is need - shiftQuad(quad,0,WisF); - } - else { - // we have to shift quad on 2 - shiftQuad(quad,2,WisF); - } - } - else { - if (nr>nl) { - // we have to shift quad on 1 - shiftQuad(quad,1,WisF); - } - else { - // we have to shift quad on 3 - shiftQuad(quad,3,WisF); - } + if ( myForcedPnts.empty() ) + { + // rotate sides to be as in the picture below and to have + // dh >= dv and nt > nb + if ( dh >= dv ) + shiftQuad( quad, ( nt > nb ) ? 0 : 2 ); + else + shiftQuad( quad, ( nr > nl ) ? 1 : 3 ); + } + else + { + // rotate the quad to have nt > nb [and nr > nl] + if ( nb > nt ) + shiftQuad ( quad, nr > nl ? 1 : 2 ); + else if ( nr > nl ) + shiftQuad( quad, nb == nt ? 1 : 0 ); + else if ( nl > nr ) + shiftQuad( quad, 3 ); } - nb = quad->side[0]->NbPoints(); - nr = quad->side[1]->NbPoints(); - nt = quad->side[2]->NbPoints(); - nl = quad->side[3]->NbPoints(); + nb = quad->side[0].NbPoints(); + nr = quad->side[1].NbPoints(); + nt = quad->side[2].NbPoints(); + nl = quad->side[3].NbPoints(); dh = abs(nb-nt); dv = abs(nr-nl); int nbh = Max(nb,nt); - int nbv = Max(nr,nl); + int nbv = Max(nr,nl); int addh = 0; int addv = 0; + // Orientation of face and 3 main domain for future faces // ----------- Old version --------------- - // orientation of face and 3 main domain for future faces // 0 top 1 // 1------------1 // | | | | - // | | | | + // | |C | | // | L | | R | - // left | | | | rigth + // left | |__| | right // | / \ | // | / C \ | // |/ \| @@ -1502,38 +2194,360 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, // 0 bottom 1 // ----------- New version --------------- - // orientation of face and 3 main domain for future faces // 0 top 1 // 1------------1 - // | |____| | + // | |__| | // | / \ | // | / C \ | - // left |/________\| rigth - // | | + // left |/________\| right // | | + // | C | // | | // 0------------0 // 0 bottom 1 - if (dh>dv) { - addv = (dh-dv)/2; - nbv = nbv + addv; - } - else { // dv>=dh - addh = (dv-dh)/2; - nbh = nbh + addh; + + //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 vector& uv_eb_vec = quad->side[0].GetUVPtStruct(true,0); + const vector& uv_er_vec = quad->side[1].GetUVPtStruct(false,1); + const vector& uv_et_vec = quad->side[2].GetUVPtStruct(true,1); + const vector& uv_el_vec = quad->side[3].GetUVPtStruct(false,0); + if (uv_eb_vec.empty() || + uv_er_vec.empty() || + uv_et_vec.empty() || + uv_el_vec.empty()) + return error(COMPERR_BAD_INPUT_MESH); } + FaceQuadStruct::SideIterator uv_eb, uv_er, uv_et, uv_el; + uv_eb.Init( quad->side[0] ); + uv_er.Init( quad->side[1] ); + uv_et.Init( quad->side[2] ); + uv_el.Init( quad->side[3] ); - const vector& uv_eb = quad->side[0]->GetUVPtStruct(true,0); - const vector& uv_er = quad->side[1]->GetUVPtStruct(false,1); - const vector& uv_et = quad->side[2]->GetUVPtStruct(true,1); - const vector& uv_el = quad->side[3]->GetUVPtStruct(false,0); + gp_UV a0,a1,a2,a3, p0,p1,p2,p3, uv; + double x,y; - if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl) - return error(COMPERR_BAD_INPUT_MESH); + a0 = uv_eb[ 0 ].UV(); + a1 = uv_er[ 0 ].UV(); + a2 = uv_er[ nr-1 ].UV(); + a3 = uv_et[ 0 ].UV(); - if ( myNeedSmooth ) - UpdateDegenUV( quad ); + if ( !myForcedPnts.empty() ) + { + if ( dv != 0 && dh != 0 ) // here myQuadList.size() == 1 + { + const int dmin = Min( dv, dh ); + + // Make a side separating domains L and Cb + StdMeshers_FaceSidePtr sideLCb; + UVPtStruct p3dom; // a point where 3 domains meat + { // dmin + vector pointsLCb( dmin+1 ); // 1--------1 + pointsLCb[0] = uv_eb[0]; // | | | + for ( int i = 1; i <= dmin; ++i ) // | |Ct| + { // | L | | + x = uv_et[ i ].normParam; // | |__| + y = uv_er[ i ].normParam; // | / | + p0 = quad->side[0].grid->Value2d( x ).XY(); // | / Cb |dmin + p1 = uv_er[ i ].UV(); // |/ | + p2 = uv_et[ i ].UV(); // 0--------0 + p3 = quad->side[3].grid->Value2d( y ).XY(); + uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); + pointsLCb[ i ].u = uv.X(); + pointsLCb[ i ].v = uv.Y(); + } + 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; + { + vector pointsLCt( nl ); + pointsLCt[0] = p3dom; + pointsLCt.back() = uv_et[ dmin ]; + x = uv_et[ dmin ].normParam; + p0 = quad->side[0].grid->Value2d( x ).XY(); + p2 = uv_et[ dmin ].UV(); + double y0 = uv_er[ dmin ].normParam; + for ( int i = 1; i < nl-1; ++i ) + { + y = y0 + i / ( nl-1. ) * ( 1. - y0 ); + p1 = quad->side[1].grid->Value2d( y ).XY(); + p3 = quad->side[3].grid->Value2d( y ).XY(); + uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); + pointsLCt[ i ].u = uv.X(); + pointsLCt[ i ].v = uv.Y(); + } + sideLCt = StdMeshers_FaceSide::New( pointsLCt, aFace ); + } + // Make a side separating domains Cb and Ct + StdMeshers_FaceSidePtr sideCbCt; + { + vector pointsCbCt( nb ); + pointsCbCt[0] = p3dom; + pointsCbCt.back() = uv_er[ dmin ]; + y = uv_er[ dmin ].normParam; + p1 = uv_er[ dmin ].UV(); + p3 = quad->side[3].grid->Value2d( y ).XY(); + double x0 = uv_et[ dmin ].normParam; + for ( int i = 1; i < nb-1; ++i ) + { + x = x0 + i / ( nb-1. ) * ( 1. - x0 ); + p2 = quad->side[2].grid->Value2d( x ).XY(); + p0 = quad->side[0].grid->Value2d( x ).XY(); + uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); + pointsCbCt[ i ].u = uv.X(); + pointsCbCt[ i ].v = uv.Y(); + } + sideCbCt = StdMeshers_FaceSide::New( pointsCbCt, aFace ); + } + // Make Cb quad + FaceQuadStruct* qCb = new FaceQuadStruct( quad->face, "Cb" ); + myQuadList.push_back( FaceQuadStruct::Ptr( qCb )); + qCb->side.resize(4); + qCb->side[0] = quad->side[0]; + qCb->side[1] = quad->side[1]; + qCb->side[2] = sideCbCt; + qCb->side[3] = sideLCb; + qCb->side[1].to = dmin+1; + // Make L quad + FaceQuadStruct* qL = new FaceQuadStruct( quad->face, "L" ); + myQuadList.push_back( FaceQuadStruct::Ptr( qL )); + qL->side.resize(4); + qL->side[0] = sideLCb; + qL->side[1] = sideLCt; + qL->side[2] = quad->side[2]; + qL->side[3] = quad->side[3]; + qL->side[2].to = dmin+1; + // Make Ct from the main quad + FaceQuadStruct::Ptr qCt = quad; + qCt->side[0] = sideCbCt; + qCt->side[3] = sideLCt; + qCt->side[1].from = dmin; + qCt->side[2].from = dmin; + qCt->uv_grid.clear(); + qCt->name = "Ct"; + + // Connect sides + qCb->side[3].AddContact( dmin, & qCb->side[2], 0 ); + qCb->side[3].AddContact( dmin, & qCt->side[3], 0 ); + qCt->side[3].AddContact( 0, & qCt->side[0], 0 ); + qCt->side[0].AddContact( 0, & qL ->side[0], dmin ); + qL ->side[0].AddContact( dmin, & qL ->side[1], 0 ); + qL ->side[0].AddContact( dmin, & qCb->side[2], 0 ); + + if ( dh == dv ) + return computeQuadDominant( aMesh, aFace ); + else + return computeQuadPref( aMesh, aFace, qCt ); + + } // 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 dt = quad->side[2].IsReversed() ? -1 : +1; + //const int dl = quad->side[3].IsReversed() ? -1 : +1; + + // Case dv == 0, here possibly myQuadList.size() > 1 + // + // lw nb lw = dh/2 + // +------------+ + // | | | | + // | | Ct | | + // | L | | R | + // | |____| | + // | / \ | + // | / Cb \ | + // |/ \| + // +------------+ + const int lw = dh/2; // lateral width + + double yCbL, yCbR; + { + double lL = quad->side[3].Length(); + double lLwL = quad->side[2].Length( tfrom, + tfrom + ( lw ) * dt ); + yCbL = lLwL / ( lLwL + lL ); + + double lR = quad->side[1].Length(); + double lLwR = quad->side[2].Length( tfrom + ( lw + nb-1 ) * dt, + tfrom + ( lw + nb-1 + lw ) * dt); + yCbR = lLwR / ( lLwR + lR ); + } + // Make sides separating domains Cb and L and R + StdMeshers_FaceSidePtr sideLCb, sideRCb; + UVPtStruct pTBL, pTBR; // points where 3 domains meat + { + vector pointsLCb( lw+1 ), pointsRCb( lw+1 ); + pointsLCb[0] = uv_eb[ 0 ]; + pointsRCb[0] = uv_eb[ nb-1 ]; + for ( int i = 1, i2 = nt-2; i <= lw; ++i, --i2 ) + { + x = quad->side[2].Param( i ); + y = yCbL * i / lw; + p0 = quad->side[0].Value2d( x ); + p1 = quad->side[1].Value2d( y ); + p2 = uv_et[ i ].UV(); + p3 = quad->side[3].Value2d( y ); + uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); + pointsLCb[ i ].u = uv.X(); + pointsLCb[ i ].v = uv.Y(); + pointsLCb[ i ].x = x; + + x = quad->side[2].Param( i2 ); + y = yCbR * i / lw; + p1 = quad->side[1].Value2d( y ); + p0 = quad->side[0].Value2d( x ); + p2 = uv_et[ i2 ].UV(); + p3 = quad->side[3].Value2d( y ); + uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); + pointsRCb[ i ].u = uv.X(); + pointsRCb[ i ].v = uv.Y(); + pointsRCb[ i ].x = x; + } + sideLCb = StdMeshers_FaceSide::New( pointsLCb, aFace ); + 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; + { + vector pointsLCt( nl ), pointsRCt( nl ); + pointsLCt[0] = pTBL; + pointsLCt.back() = uv_et[ lw ]; + pointsRCt[0] = pTBR; + pointsRCt.back() = uv_et[ lw + nb - 1 ]; + x = pTBL.x; + p0 = quad->side[0].Value2d( x ); + p2 = uv_et[ lw ].UV(); + int iR = lw + nb - 1; + double xR = pTBR.x; + gp_UV p0R = quad->side[0].Value2d( xR ); + gp_UV p2R = uv_et[ iR ].UV(); + for ( int i = 1; i < nl-1; ++i ) + { + y = yCbL + ( 1. - yCbL ) * i / (nl-1.); + p1 = quad->side[1].Value2d( y ); + p3 = quad->side[3].Value2d( y ); + uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); + pointsLCt[ i ].u = uv.X(); + pointsLCt[ i ].v = uv.Y(); + + y = yCbR + ( 1. - yCbR ) * i / (nl-1.); + p1 = quad->side[1].Value2d( y ); + p3 = quad->side[3].Value2d( y ); + uv = calcUV( xR,y, a0,a1,a2,a3, p0R,p1,p2R,p3 ); + pointsRCt[ i ].u = uv.X(); + pointsRCt[ i ].v = uv.Y(); + } + sideLCt = StdMeshers_FaceSide::New( pointsLCt, aFace ); + sideRCt = StdMeshers_FaceSide::New( pointsRCt, aFace ); + } + // Make a side separating domains Cb and Ct + StdMeshers_FaceSidePtr sideCbCt; + { + vector pointsCbCt( nb ); + pointsCbCt[0] = pTBL; + pointsCbCt.back() = pTBR; + p1 = quad->side[1].Value2d( yCbR ); + p3 = quad->side[3].Value2d( yCbL ); + for ( int i = 1; i < nb-1; ++i ) + { + x = quad->side[2].Param( i + lw ); + y = yCbL + ( yCbR - yCbL ) * i / (nb-1.); + p2 = uv_et[ i + lw ].UV(); + p0 = quad->side[0].Value2d( x ); + uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); + pointsCbCt[ i ].u = uv.X(); + pointsCbCt[ i ].v = uv.Y(); + } + sideCbCt = StdMeshers_FaceSide::New( pointsCbCt, aFace ); + } + // Make Cb quad + FaceQuadStruct* qCb = new FaceQuadStruct( quad->face, "Cb" ); + myQuadList.push_back( FaceQuadStruct::Ptr( qCb )); + qCb->side.resize(4); + qCb->side[0] = quad->side[0]; + qCb->side[1] = sideRCb; + qCb->side[2] = sideCbCt; + qCb->side[3] = sideLCb; + // Make L quad + FaceQuadStruct* qL = new FaceQuadStruct( quad->face, "L" ); + myQuadList.push_back( FaceQuadStruct::Ptr( qL )); + qL->side.resize(4); + qL->side[0] = sideLCb; + qL->side[1] = sideLCt; + qL->side[2] = quad->side[2]; + qL->side[3] = quad->side[3]; + qL->side[2].to = ( lw + 1 ) * dt + tfrom; + // Make R quad + FaceQuadStruct* qR = new FaceQuadStruct( quad->face, "R" ); + myQuadList.push_back( FaceQuadStruct::Ptr( qR )); + qR->side.resize(4); + qR->side[0] = sideRCb; + qR->side[0].from = lw; + qR->side[0].to = -1; + qR->side[0].di = -1; + qR->side[1] = quad->side[1]; + qR->side[2] = quad->side[2]; + qR->side[2].from = ( lw + nb-1 ) * dt + tfrom; + qR->side[3] = sideRCt; + // Make Ct from the main quad + FaceQuadStruct::Ptr qCt = quad; + qCt->side[0] = sideCbCt; + qCt->side[1] = sideRCt; + qCt->side[2].from = ( lw ) * dt + tfrom; + qCt->side[2].to = ( lw + nb ) * dt + tfrom; + qCt->side[3] = sideLCt; + qCt->uv_grid.clear(); + qCt->name = "Ct"; + + // Connect sides + qCb->side[3].AddContact( lw, & qCb->side[2], 0 ); + qCb->side[3].AddContact( lw, & qCt->side[3], 0 ); + qCt->side[3].AddContact( 0, & qCt->side[0], 0 ); + qCt->side[0].AddContact( 0, & qL ->side[0], lw ); + qL ->side[0].AddContact( lw, & qL ->side[1], 0 ); + qL ->side[0].AddContact( lw, & qCb->side[2], 0 ); + // + qCb->side[1].AddContact( lw, & qCb->side[2], nb-1 ); + qCb->side[1].AddContact( lw, & qCt->side[1], 0 ); + qCt->side[0].AddContact( nb-1, & qCt->side[1], 0 ); + qCt->side[0].AddContact( nb-1, & qR ->side[0], lw ); + qR ->side[3].AddContact( 0, & qR ->side[0], lw ); + qR ->side[3].AddContact( 0, & qCb->side[2], nb-1 ); + + return computeQuadDominant( aMesh, aFace ); + + } // if ( !myForcedPnts.empty() ) + + if ( dh > dv ) { + addv = (dh-dv)/2; + nbv = nbv + addv; + } + else { // dv >= dh + addh = (dv-dh)/2; + nbh = nbh + addh; + } // arrays for normalized params TColStd_SequenceOfReal npb, npr, npt, npl; @@ -1550,7 +2564,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 @@ -1566,14 +2580,9 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, npl.InsertAfter(1,npl.Value(2)-dpr); } } - - gp_XY a0(uv_eb.front().u, uv_eb.front().v); - gp_XY a1(uv_eb.back().u, uv_eb.back().v); - gp_XY a2(uv_et.back().u, uv_et.back().v); - gp_XY a3(uv_et.front().u, uv_et.front().v); 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; @@ -1586,7 +2595,7 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, NodesL.SetValue(1,j,uv_el[j-1].node); if (dl>0) { // add top nodes - for (i=1; i<=dl; i++) + for (i=1; i<=dl; i++) NodesL.SetValue(i+1,nl,uv_et[i].node); // create and add needed nodes TColgp_SequenceOfXY UVtmp; @@ -1621,16 +2630,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), + myHelper->AddFace(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); - } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1), - NodesL.Value(i+1,j+1), NodesL.Value(i+1,j)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } } } @@ -1641,15 +2642,15 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, UVL.Append(gp_UV (uv_el[i].u, uv_el[i].v)); } } - + // step2: create faces for right domain StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr); // add right nodes - for (j=1; j<=nr; j++) + for (j=1; j<=nr; j++) NodesR.SetValue(1,j,uv_er[nr-j].node); if (dr>0) { // add top nodes - for (i=1; i<=dr; i++) + for (i=1; i<=dr; i++) NodesR.SetValue(i+1,1,uv_et[nt-1-i].node); // create and add needed nodes TColgp_SequenceOfXY UVtmp; @@ -1684,16 +2685,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), + myHelper->AddFace(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); - } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1), - NodesR.Value(i+1,j+1), NodesR.Value(i+1,j)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } } } @@ -1704,7 +2697,7 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, UVR.Append(gp_UV(uv_er[i].u, uv_er[i].v)); } } - + // step3: create faces for central domain StdMeshers_Array2OfNode NodesC(1,nb,1,nbv); // add first line using NodesL @@ -1718,12 +2711,12 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, for (i=1; iAddFace(NodesC.Value(i,j), NodesC.Value(i+1,j), + myHelper->AddFace(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); - } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1), - NodesC.Value(i+1,j+1), NodesC.Value(i+1,j)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } } } @@ -1801,16 +2786,8 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, for (j=1; jAddFace(NodesBRD.Value(i,j), NodesBRD.Value(i+1,j), + myHelper->AddFace(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); - } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i,j+1), - NodesBRD.Value(i+1,j+1), NodesBRD.Value(i+1,j)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } } } @@ -1829,7 +2806,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); @@ -1855,7 +2832,7 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, double yy1 = y1 + dy1*i; double dyy = yy1 - yy0; for (j=1; j<=nb; j++) { - double x = npt.Value(i+1+drl) + + double x = npt.Value(i+1+drl) + npb.Value(j) * (npt.Value(nt-i) - npt.Value(i+1+drl)); double y = yy0 + dyy*x; gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3); @@ -1898,7 +2875,7 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, double yy1 = y1 + dy1*i; double dyy = yy1 - yy0; for (j=1; j<=nb; j++) { - double x = npt.Value(i+1) + + double x = npt.Value(i+1) + npb.Value(j) * (npt.Value(nt-i-drl) - npt.Value(i+1)); double y = yy0 + dyy*x; gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3); @@ -1913,16 +2890,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), + myHelper->AddFace(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); - } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1), - NodesC.Value(i+1,j+1), NodesC.Value(i+1,j)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } } } // end nrAddFace(NodesLast.Value(i,1), NodesLast.Value(i+1,1), + myHelper->AddFace(NodesLast.Value(i,1), NodesLast.Value(i+1,1), NodesLast.Value(i+1,2), NodesLast.Value(i,2)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); - } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesLast.Value(i,1), NodesLast.Value(i,2), - NodesLast.Value(i+1,2), NodesLast.Value(i+1,2)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } } } // if ((drl+addv) > 0) @@ -1973,13 +2934,13 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, */ //======================================================================= -bool StdMeshers_Quadrangle_2D::EvaluateQuadPref(SMESH_Mesh & aMesh, +bool StdMeshers_Quadrangle_2D::evaluateQuadPref(SMESH_Mesh & aMesh, const TopoDS_Shape& aShape, - std::vector& aNbNodes, - MapShapeNbElems& aResMap, - bool IsQuadratic) + std::vector& aNbNodes, + 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; @@ -2101,36 +3062,29 @@ bool StdMeshers_Quadrangle_2D::EvaluateQuadPref(SMESH_Mesh & aMesh, return true; } - //============================================================================= /*! Split quadrangle in to 2 triangles by smallest diagonal * */ //============================================================================= -void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS, - int theFaceID, - const SMDS_MeshNode* theNode1, - const SMDS_MeshNode* theNode2, - const SMDS_MeshNode* theNode3, - const SMDS_MeshNode* theNode4) + +void StdMeshers_Quadrangle_2D::splitQuadFace(SMESHDS_Mesh * theMeshDS, + int theFaceID, + const SMDS_MeshNode* theNode1, + const SMDS_MeshNode* theNode2, + const SMDS_MeshNode* theNode3, + const SMDS_MeshNode* theNode4) { - gp_Pnt a(theNode1->X(),theNode1->Y(),theNode1->Z()); - gp_Pnt b(theNode2->X(),theNode2->Y(),theNode2->Z()); - gp_Pnt c(theNode3->X(),theNode3->Y(),theNode3->Z()); - gp_Pnt d(theNode4->X(),theNode4->Y(),theNode4->Z()); - SMDS_MeshFace* face; - if (a.Distance(c) > b.Distance(d)){ - face = myHelper->AddFace(theNode2, theNode4 , theNode1); - if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID); - face = myHelper->AddFace(theNode2, theNode3, theNode4); - if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID); - - } - 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); + if ( SMESH_TNodeXYZ( theNode1 ).SquareDistance( theNode3 ) > + SMESH_TNodeXYZ( theNode2 ).SquareDistance( theNode4 ) ) + { + myHelper->AddFace(theNode2, theNode4 , theNode1); + myHelper->AddFace(theNode2, theNode3, theNode4); + } + else + { + myHelper->AddFace(theNode1, theNode2 ,theNode3); + myHelper->AddFace(theNode1, theNode3, theNode4); } } @@ -2145,8 +3099,8 @@ namespace SMESH_MesherHelper* helper, Handle(Geom_Surface) S) { - const vector& uv_eb = quad->side[QUAD_BOTTOM_SIDE]->GetUVPtStruct(); - const vector& uv_et = quad->side[QUAD_TOP_SIDE ]->GetUVPtStruct(); + const vector& uv_eb = quad->side[QUAD_BOTTOM_SIDE].GetUVPtStruct(); + const vector& uv_et = quad->side[QUAD_TOP_SIDE ].GetUVPtStruct(); double rBot = ( uv_eb.size() - 1 ) * uvPt.normParam; double rTop = ( uv_et.size() - 1 ) * uvPt.normParam; int iBot = int( rBot ); @@ -2157,9 +3111,9 @@ namespace gp_UV uv = calcUV(/*x,y=*/x, y, /*a0,...=*/UVs[UV_A0], UVs[UV_A1], UVs[UV_A2], UVs[UV_A3], - /*p0=*/quad->side[QUAD_BOTTOM_SIDE]->Value2d( x ).XY(), + /*p0=*/quad->side[QUAD_BOTTOM_SIDE].grid->Value2d( x ).XY(), /*p1=*/UVs[ UV_R ], - /*p2=*/quad->side[QUAD_TOP_SIDE ]->Value2d( x ).XY(), + /*p2=*/quad->side[QUAD_TOP_SIDE ].grid->Value2d( x ).XY(), /*p3=*/UVs[ UV_L ]); gp_Pnt P = S->Value( uv.X(), uv.Y() ); uvPt.u = uv.X(); @@ -2316,19 +3270,19 @@ namespace * Implementation of Reduced algorithm (meshing with quadrangles only) */ //======================================================================= -bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, - const TopoDS_Shape& aShape, + +bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, + const TopoDS_Face& aFace, FaceQuadStruct::Ptr quad) { - SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); - const TopoDS_Face& F = TopoDS::Face(aShape); - Handle(Geom_Surface) S = BRep_Tool::Surface(F); - int i,j,geomFaceID = meshDS->ShapeToIndex(F); + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + Handle(Geom_Surface) S = BRep_Tool::Surface(aFace); + int i,j,geomFaceID = meshDS->ShapeToIndex(aFace); - int nb = quad->side[0]->NbPoints(); - int nr = quad->side[1]->NbPoints(); - int nt = quad->side[2]->NbPoints(); - int nl = quad->side[3]->NbPoints(); + int nb = quad->side[0].NbPoints(); // bottom + int nr = quad->side[1].NbPoints(); // right + int nt = quad->side[2].NbPoints(); // top + int nl = quad->side[3].NbPoints(); // left // Simple Reduce 10->8->6->4 (3 steps) Multiple Reduce 10->4 (1 step) // @@ -2374,16 +3328,24 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, } // number of rows and columns - int nrows = nr1 - 1; + int nrows = nr1 - 1; int ncol_top = nt1 - 1; int ncol_bot = nb1 - 1; // number of rows needed to reduce ncol_bot to ncol_top using simple 3->1 "tree" (see below) - int nrows_tree31 = int( log( (double)(ncol_bot / ncol_top) ) / log((double) 3 )); // = log x base 3 + int nrows_tree31 = + int( ceil( log( double(ncol_bot) / ncol_top) / log( 3.))); // = log x base 3 if ( nrows < nrows_tree31 ) + { MultipleReduce = true; + error( COMPERR_WARNING, + SMESH_Comment("To use 'Reduced' transition, " + "number of face rows should be at least ") + << nrows_tree31 << ". Actual number of face rows is " << nrows << ". " + "'Quadrangle preference (reversed)' transion has been used."); + } } - if (MultipleReduce) { // == ComputeQuadPref QUAD_QUADRANGLE_PREF_REVERSED + if (MultipleReduce) { // == computeQuadPref QUAD_QUADRANGLE_PREF_REVERSED //================================================== int dh = abs(nb-nt); int dv = abs(nr-nl); @@ -2391,28 +3353,28 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, if (dh >= dv) { if (nt > nb) { // it is a base case => not shift quad but may be replacement is need - shiftQuad(quad,0,true); + shiftQuad(quad,0); } else { // we have to shift quad on 2 - shiftQuad(quad,2,true); + shiftQuad(quad,2); } } else { if (nr > nl) { // we have to shift quad on 1 - shiftQuad(quad,1,true); + shiftQuad(quad,1); } else { // we have to shift quad on 3 - shiftQuad(quad,3,true); + shiftQuad(quad,3); } } - nb = quad->side[0]->NbPoints(); - nr = quad->side[1]->NbPoints(); - nt = quad->side[2]->NbPoints(); - nl = quad->side[3]->NbPoints(); + nb = quad->side[0].NbPoints(); + nr = quad->side[1].NbPoints(); + nt = quad->side[2].NbPoints(); + nl = quad->side[3].NbPoints(); dh = abs(nb-nt); dv = abs(nr-nl); int nbh = Max(nb,nt); @@ -2429,17 +3391,15 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, nbh = nbh + addh; } - const vector& uv_eb = quad->side[0]->GetUVPtStruct(true,0); - const vector& uv_er = quad->side[1]->GetUVPtStruct(false,1); - const vector& uv_et = quad->side[2]->GetUVPtStruct(true,1); - const vector& uv_el = quad->side[3]->GetUVPtStruct(false,0); + const vector& uv_eb = quad->side[0].GetUVPtStruct(true,0); + const vector& uv_er = quad->side[1].GetUVPtStruct(false,1); + 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); - if ( myNeedSmooth ) - UpdateDegenUV( quad ); - // arrays for normalized params TColStd_SequenceOfReal npb, npr, npt, npl; for (j = 0; j < nb; j++) { @@ -2462,7 +3422,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, // | | | | // | | | | // | L | | R | - // left | | | | rigth + // left | | | | right // | / \ | // | / C \ | // |/ \| @@ -2489,7 +3449,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; @@ -2537,10 +3497,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)); } } } @@ -2592,10 +3550,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)); } } } @@ -2624,7 +3580,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, // add bottom nodes (first columns) for (i=2; 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)); } } - // TODO ??? } // end Multiple Reduce implementation else { // Simple Reduce (!MultipleReduce) //========================================================= @@ -2673,24 +3626,24 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, } else { // we have to shift quad on 2 - shiftQuad(quad,2,true); + shiftQuad(quad,2); } } else { if (nl > nr) { // we have to shift quad on 1 - shiftQuad(quad,1,true); + shiftQuad(quad,1); } else { // we have to shift quad on 3 - shiftQuad(quad,3,true); + shiftQuad(quad,3); } } - nb = quad->side[0]->NbPoints(); - nr = quad->side[1]->NbPoints(); - nt = quad->side[2]->NbPoints(); - nl = quad->side[3]->NbPoints(); + nb = quad->side[0].NbPoints(); + nr = quad->side[1].NbPoints(); + nt = quad->side[2].NbPoints(); + nl = quad->side[3].NbPoints(); // number of rows and columns int nrows = nr - 1; // and also == nl - 1 @@ -2752,16 +3705,15 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, } } - const vector& uv_eb = quad->side[0]->GetUVPtStruct(true,0); - const vector& uv_er = quad->side[1]->GetUVPtStruct(false,1); - const vector& uv_et = quad->side[2]->GetUVPtStruct(true,1); - const vector& uv_el = quad->side[3]->GetUVPtStruct(false,0); + const vector& uv_eb = quad->side[0].GetUVPtStruct(true,0); + const vector& uv_er = quad->side[1].GetUVPtStruct(false,1); + 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 ); @@ -2770,7 +3722,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; @@ -2803,10 +3759,10 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, vector nb_col_by_row; - int delta_all = nb - nt; + int delta_all = nb - nt; int delta_one_col = nrows * 2; - int nb_col = delta_all / delta_one_col; - int remainder = delta_all - nb_col * delta_one_col; + int nb_col = delta_all / delta_one_col; + int remainder = delta_all - nb_col * delta_one_col; if (remainder > 0) { nb_col++; } @@ -2817,13 +3773,13 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, nb_col = ( nt - 1 ) / col_top_size; nb_col_by_row.resize( nrows, nb_col ); int nbrows_not_full = nrows - 1; - int cur_top_size = nt - 1; + int cur_top_size = nt - 1; remainder = delta_all - nb_col * delta_one_col; while ( remainder > 0 ) { - delta_one_col = nbrows_not_full * 2; - int nb_col_add = remainder / delta_one_col; - cur_top_size += 2 * nb_col_by_row[ nbrows_not_full ]; + delta_one_col = nbrows_not_full * 2; + int nb_col_add = remainder / delta_one_col; + cur_top_size += 2 * nb_col_by_row[ nbrows_not_full ]; int nb_col_free = cur_top_size / col_top_size - nb_col_by_row[ nbrows_not_full-1 ]; if ( nb_col_add > nb_col_free ) nb_col_add = nb_col_free; @@ -3246,7 +4202,8 @@ namespace // data for smoothing */ struct TSmoothNode { - gp_XY _uv; + gp_XY _uv; + gp_XYZ _xyz; vector< TTriangle > _triangles; // if empty, then node is not movable }; // -------------------------------------------------------------------------------- @@ -3256,6 +4213,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; + } } //================================================================================ @@ -3266,43 +4235,69 @@ namespace // data for smoothing */ //================================================================================ -void StdMeshers_Quadrangle_2D::UpdateDegenUV(FaceQuadStruct::Ptr quad) +void StdMeshers_Quadrangle_2D::updateDegenUV(FaceQuadStruct::Ptr quad) { - for ( unsigned i = 0; i < quad->side.size(); ++i ) - { - StdMeshers_FaceSide* side = quad->side[i]; - const vector& uvVec = side->GetUVPtStruct(); + if ( myNeedSmooth ) - // find which end of the side is on degenerated shape - int degenInd = -1; - if ( myHelper->IsDegenShape( uvVec[0].node->getshapeId() )) - degenInd = 0; - else if ( myHelper->IsDegenShape( uvVec.back().node->getshapeId() )) - degenInd = uvVec.size() - 1; - else - continue; + // Set UV of nodes on degenerated VERTEXes in the middle of degenerated EDGE + // -------------------------------------------------------------------------- + for ( unsigned i = 0; i < quad->side.size(); ++i ) + { + const vector& uvVec = quad->side[i].GetUVPtStruct(); + + // find which end of the side is on degenerated shape + int degenInd = -1; + if ( myHelper->IsDegenShape( uvVec[0].node->getshapeId() )) + degenInd = 0; + else if ( myHelper->IsDegenShape( uvVec.back().node->getshapeId() )) + degenInd = uvVec.size() - 1; + else + continue; + + // find another side sharing the degenerated shape + bool isPrev = ( degenInd == 0 ); + if ( i >= QUAD_TOP_SIDE ) + isPrev = !isPrev; + int i2 = ( isPrev ? ( i + 3 ) : ( i + 1 )) % 4; + const vector& uvVec2 = quad->side[ i2 ].GetUVPtStruct(); + int degenInd2 = -1; + if ( uvVec[ degenInd ].node == uvVec2.front().node ) + degenInd2 = 0; + else if ( uvVec[ degenInd ].node == uvVec2.back().node ) + degenInd2 = uvVec2.size() - 1; + else + throw SALOME_Exception( LOCALIZED( "Logical error" )); - // find another side sharing the degenerated shape - bool isPrev = ( degenInd == 0 ); - if ( i >= QUAD_TOP_SIDE ) - isPrev = !isPrev; - int i2 = ( isPrev ? ( i + 3 ) : ( i + 1 )) % 4; - StdMeshers_FaceSide* side2 = quad->side[ i2 ]; - const vector& uvVec2 = side2->GetUVPtStruct(); - int degenInd2 = -1; - if ( uvVec[ degenInd ].node == uvVec2[0].node ) - degenInd2 = 0; - else if ( uvVec[ degenInd ].node == uvVec2.back().node ) - degenInd2 = uvVec2.size() - 1; - else - throw SALOME_Exception( LOCALIZED( "Logical error" )); + // move UV in the middle + uvPtStruct& uv1 = const_cast( uvVec [ degenInd ]); + uvPtStruct& uv2 = const_cast( uvVec2[ degenInd2 ]); + uv1.u = uv2.u = 0.5 * ( uv1.u + uv2.u ); + uv1.v = uv2.v = 0.5 * ( uv1.v + uv2.v ); + } - // move UV in the middle - uvPtStruct& uv1 = const_cast( uvVec [ degenInd ]); - uvPtStruct& uv2 = const_cast( uvVec2[ degenInd2 ]); - uv1.u = uv2.u = 0.5 * ( uv1.u + uv2.u ); - uv1.v = uv2.v = 0.5 * ( uv1.v + uv2.v ); - } + 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 ( size_t i = 0; i < quad->side.size(); ++i ) + { + StdMeshers_FaceSidePtr degSide = quad->side[i]; + if ( !myHelper->IsDegenShape( degSide->EdgeID(0) )) + continue; + StdMeshers_FaceSidePtr oppSide = quad->side[( i+2 ) % quad->side.size() ]; + if ( degSide->NbSegments() == oppSide->NbSegments() ) + continue; + + // make new side data + const vector& uvVecDegOld = degSide->GetUVPtStruct(); + const SMDS_MeshNode* n = uvVecDegOld[0].node; + Handle(Geom2d_Curve) c2d = degSide->Curve2d(0); + double f = degSide->FirstU(0), l = degSide->LastU(0); + gp_Pnt2d p1 = uvVecDegOld.front().UV(); + gp_Pnt2d p2 = uvVecDegOld.back().UV(); + + quad->side[i] = StdMeshers_FaceSide::New( oppSide.get(), n, &p1, &p2, c2d, f, l ); + } } //================================================================================ @@ -3311,102 +4306,189 @@ void StdMeshers_Quadrangle_2D::UpdateDegenUV(FaceQuadStruct::Ptr quad) */ //================================================================================ -void StdMeshers_Quadrangle_2D::Smooth (FaceQuadStruct::Ptr quad) +void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad) { if ( !myNeedSmooth ) return; - // Get nodes to smooth - - typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap; - TNo2SmooNoMap smooNoMap; - - const TopoDS_Face& geomFace = TopoDS::Face( myHelper->GetSubShape() ); 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 double tol = BRep_Tool::Tolerance( quad->face ); + Handle(ShapeAnalysis_Surface) surface = myHelper->GetSurface( quad->face ); + + if ( myHelper->HasDegeneratedEdges() && myForcedPnts.empty() ) { - const SMDS_MeshNode* node = nIt->next(); - TSmoothNode & sNode = smooNoMap[ node ]; - sNode._uv = myHelper->GetNodeUV( geomFace, node ); + // "smooth" by computing node positions using 3D TFI and further projection - // set sNode._triangles - SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face ); - while ( fIt->more() ) + list< FaceQuadStruct::Ptr >::iterator q = myQuadList.begin(); + for ( ; q != myQuadList.end() ; ++q ) { - 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 ])); + quad = *q; + 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(); + } + } } } - // 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.SetCoord( uvVec[j].u, uvVec[j].v ); + 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 - - // compute a new UV - gp_XY newUV (0,0); - for ( unsigned i = 0; i < sNode._triangles.size(); ++i ) - newUV += sNode._triangles[i]._n1->_uv; - newUV /= sNode._triangles.size(); + for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn ) + { + TSmoothNode& sNode = n2sn->second; + if ( sNode._triangles.empty() ) + continue; // not movable node - // check validity of the newUV - bool isValid = true; - for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i ) - isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward ); + gp_XY newUV; + bool isValid = false; + bool use3D = ( iLoop > 2 ); // 3 loops in 2D and 2, in 3D - if ( isValid ) - sNode._uv = newUV; + 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 + 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 ) + { + sNode._uv = newUV; + sNode._xyz = surface->Value( newUV ).XYZ(); + } + } } - } - // Set new XYZ to the smoothed nodes + // Set new XYZ to the smoothed nodes - Handle(Geom_Surface) surface = BRep_Tool::Surface( geomFace ); - - 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 @@ -3422,14 +4504,1277 @@ 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; +} + +//================================================================================ +/*! + * \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, myHelper ); + 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]; + 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: + { + 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 + } + isBad = ( sign1 * okSign < 0 ); + break; + } + case 3: + { + double sign = getArea( uv[0], uv[1], uv[2] ); + 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 ); + + return !isOK; + } + + return isOK; +} + +//================================================================================ +/*! + * \brief Constructor of a side of quad + */ +//================================================================================ + +FaceQuadStruct::Side::Side(StdMeshers_FaceSidePtr theGrid) + : grid(theGrid), from(0), to(theGrid ? theGrid->NbPoints() : 0 ), di(1), nbNodeOut(0) +{ +} + +//============================================================================= +/*! + * \brief Constructor of a quad + */ +//============================================================================= + +FaceQuadStruct::FaceQuadStruct(const TopoDS_Face& F, const std::string& theName) + : face( F ), name( theName ) +{ + side.reserve(4); +} + +//================================================================================ +/*! + * \brief Fills myForcedPnts + */ +//================================================================================ + +bool StdMeshers_Quadrangle_2D::getEnforcedUV() +{ + myForcedPnts.clear(); + if ( !myParams ) return true; // missing hypothesis + + std::vector< TopoDS_Shape > shapes; + std::vector< gp_Pnt > points; + myParams->GetEnforcedNodes( shapes, points ); + + TopTools_IndexedMapOfShape vMap; + for ( size_t i = 0; i < shapes.size(); ++i ) + if ( !shapes[i].IsNull() ) + TopExp::MapShapes( shapes[i], TopAbs_VERTEX, vMap ); + + size_t nbPoints = points.size(); + for ( int i = 1; i <= vMap.Extent(); ++i ) + points.push_back( BRep_Tool::Pnt( TopoDS::Vertex( vMap( i )))); + + // find out if all points must be in the FACE, which is so if + // myParams is a local hypothesis on the FACE being meshed + bool isStrictCheck = false; + { + SMESH_HypoFilter paramFilter( SMESH_HypoFilter::Is( myParams )); + TopoDS_Shape assignedTo; + if ( myHelper->GetMesh()->GetHypothesis( myHelper->GetSubShape(), + paramFilter, + /*ancestors=*/true, + &assignedTo )) + isStrictCheck = ( assignedTo.IsSame( myHelper->GetSubShape() )); + } + + multimap< double, ForcedPoint > sortedFP; // sort points by distance from EDGEs + + Standard_Real u1,u2,v1,v2; + const TopoDS_Face& face = TopoDS::Face( myHelper->GetSubShape() ); + const double tol = BRep_Tool::Tolerance( face ); + 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 ) + { + 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 = ") + << points[ iP ].Distance( project->Value( uv )) << " - (" + << points[ iP ].X() << ", "<< points[ iP ].Y() << ", "<< points[ iP ].Z() << " )"); + continue; + } + BRepClass_FaceClassifier clsf ( face, uv, tol ); + switch ( clsf.State() ) { + case TopAbs_IN: + { + 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; + } + case TopAbs_OUT: + { + if ( isStrictCheck && iP < nbPoints ) + return error + (COMPERR_BAD_PARMETERS, TComm("An enforced point is out of the face boundary - ") + << points[ iP ].X() << ", "<< points[ iP ].Y() << ", "<< points[ iP ].Z() << " )"); + break; + } + case TopAbs_ON: + { + if ( isStrictCheck && iP < nbPoints ) + return error + (COMPERR_BAD_PARMETERS, TComm("An enforced point is on the face boundary - ") + << points[ iP ].X() << ", "<< points[ iP ].Y() << ", "<< points[ iP ].Z() << " )"); + break; + } + default: + { + if ( isStrictCheck && iP < nbPoints ) + return error + (TComm("Classification of an enforced point ralative to the face boundary failed - ") + << points[ iP ].X() << ", "<< points[ iP ].Y() << ", "<< points[ iP ].Z() << " )"); + } + } + } + + multimap< double, ForcedPoint >::iterator d2uv = sortedFP.begin(); + for ( ; d2uv != sortedFP.end(); ++d2uv ) + myForcedPnts.push_back( (*d2uv).second ); + + return true; +} + +//================================================================================ +/*! + * \brief Splits quads by adding points of enforced nodes and create nodes on + * the sides shared by quads + */ +//================================================================================ + +bool StdMeshers_Quadrangle_2D::addEnforcedNodes() +{ + // if ( myForcedPnts.empty() ) + // return true; + + // make a map of quads sharing a side + map< StdMeshers_FaceSidePtr, vector< FaceQuadStruct::Ptr > > quadsBySide; + list< FaceQuadStruct::Ptr >::iterator quadIt = myQuadList.begin(); + for ( ; quadIt != myQuadList.end(); ++quadIt ) + for ( size_t iSide = 0; iSide < (*quadIt)->side.size(); ++iSide ) + { + if ( !setNormalizedGrid( *quadIt )) + return false; + quadsBySide[ (*quadIt)->side[iSide] ].push_back( *quadIt ); + } + + const TopoDS_Face& face = TopoDS::Face( myHelper->GetSubShape() ); + Handle(Geom_Surface) surf = BRep_Tool::Surface( face ); + + for ( size_t iFP = 0; iFP < myForcedPnts.size(); ++iFP ) + { + bool isNodeEnforced = false; + + // look for a quad enclosing an enforced point + for ( quadIt = myQuadList.begin(); quadIt != myQuadList.end(); ++quadIt ) + { + FaceQuadStruct::Ptr quad = *quadIt; + if ( !setNormalizedGrid( *quadIt )) + return false; + int i,j; + if ( !quad->findCell( myForcedPnts[ iFP ], i, j )) + continue; + + // a grid cell is found, select a node of the cell to move + // to the enforced point to and to split the quad at + multimap< double, pair< int, int > > ijByDist; + for ( int di = 0; di < 2; ++di ) + for ( int dj = 0; dj < 2; ++dj ) + { + double dist2 = ( myForcedPnts[ iFP ].uv - quad->UVPt( i+di,j+dj ).UV() ).SquareModulus(); + ijByDist.insert( make_pair( dist2, make_pair( di,dj ))); + } + // try all nodes starting from the closest one + set< FaceQuadStruct::Ptr > changedQuads; + multimap< double, pair< int, int > >::iterator d2ij = ijByDist.begin(); + for ( ; !isNodeEnforced && d2ij != ijByDist.end(); ++d2ij ) + { + int di = d2ij->second.first; + int dj = d2ij->second.second; + + // check if a node is at a side + int iSide = -1; + if ( dj== 0 && j == 0 ) + iSide = QUAD_BOTTOM_SIDE; + else if ( dj == 1 && j+2 == quad->jSize ) + iSide = QUAD_TOP_SIDE; + else if ( di == 0 && i == 0 ) + iSide = QUAD_LEFT_SIDE; + else if ( di == 1 && i+2 == quad->iSize ) + iSide = QUAD_RIGHT_SIDE; + + if ( iSide > -1 ) // ----- node is at a side + { + FaceQuadStruct::Side& side = quad->side[ iSide ]; + // check if this node can be moved + if ( quadsBySide[ side ].size() < 2 ) + continue; // its a face boundary -> can't move the node + + int quadNodeIndex = ( iSide % 2 ) ? j : i; + int sideNodeIndex = side.ToSideIndex( quadNodeIndex ); + if ( side.IsForced( sideNodeIndex )) + { + // the node is already moved to another enforced point + isNodeEnforced = quad->isEqual( myForcedPnts[ iFP ], i, j ); + continue; + } + // 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 ].node = myForcedPnts[ iFP ].node; + + updateSideUV( side, sideNodeIndex, quadsBySide ); + + // update adjacent sides + set< StdMeshers_FaceSidePtr > updatedSides; + updatedSides.insert( side ); + for ( size_t i = 0; i < side.contacts.size(); ++i ) + if ( side.contacts[i].point == sideNodeIndex ) + { + const vector< FaceQuadStruct::Ptr >& adjQuads = + quadsBySide[ *side.contacts[i].other_side ]; + if ( adjQuads.size() > 1 && + updatedSides.insert( * side.contacts[i].other_side ).second ) + { + updateSideUV( *side.contacts[i].other_side, + side.contacts[i].other_point, + quadsBySide ); + } + changedQuads.insert( adjQuads.begin(), adjQuads.end() ); + } + const vector< FaceQuadStruct::Ptr >& adjQuads = quadsBySide[ side ]; + changedQuads.insert( adjQuads.begin(), adjQuads.end() ); + + isNodeEnforced = true; + } + else // ------------------ node is inside the quad + { + i += di; + j += dj; + // make a new side passing through IJ node and split the quad + int indForced, iNewSide; + if ( quad->iSize < quad->jSize ) // split vertically + { + quad->updateUV( myForcedPnts[ iFP ].uv, i, j, /*isVert=*/true ); + indForced = j; + iNewSide = splitQuad( quad, i, 0 ); + } + else + { + quad->updateUV( myForcedPnts[ iFP ].uv, i, j, /*isVert=*/false ); + indForced = i; + iNewSide = splitQuad( quad, 0, j ); + } + 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 ); + + quadsBySide[ newSide ].push_back( quad ); + quadsBySide[ newQuad->side[0] ].push_back( newQuad ); + quadsBySide[ newQuad->side[1] ].push_back( newQuad ); + quadsBySide[ newQuad->side[2] ].push_back( newQuad ); + quadsBySide[ newQuad->side[3] ].push_back( newQuad ); + + isNodeEnforced = true; + + } // end of "node is inside the quad" + + } // loop on nodes of the cell + + // remove out-of-date uv grid of changedQuads + set< FaceQuadStruct::Ptr >::iterator qIt = changedQuads.begin(); + for ( ; qIt != changedQuads.end(); ++qIt ) + (*qIt)->uv_grid.clear(); + + if ( isNodeEnforced ) + break; + + } // loop on quads + + if ( !isNodeEnforced ) + { + if ( !myForcedPnts[ iFP ].vertex.IsNull() ) + return error(TComm("Unable to move any node to vertex #") + <GetMeshDS()->ShapeToIndex( myForcedPnts[ iFP ].vertex )); + else + return error(TComm("Unable to move any node to point ( ") + << myForcedPnts[iFP].xyz.X() << ", " + << myForcedPnts[iFP].xyz.Y() << ", " + << myForcedPnts[iFP].xyz.Z() << " )"); + } + myNeedSmooth = true; + + } // loop on enforced points + + // Compute nodes on all sides, where not yet present + + for ( quadIt = myQuadList.begin(); quadIt != myQuadList.end(); ++quadIt ) + { + FaceQuadStruct::Ptr quad = *quadIt; + for ( int iSide = 0; iSide < 4; ++iSide ) + { + FaceQuadStruct::Side & side = quad->side[ iSide ]; + if ( side.nbNodeOut > 0 ) + continue; // emulated side + vector< FaceQuadStruct::Ptr >& quadVec = quadsBySide[ side ]; + if ( quadVec.size() <= 1 ) + continue; // outer side + + 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 >= (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 = myHelper->AddNode(P.X(), P.Y(), P.Z(), 0, uvPnt.u, uvPnt.v ); + missedNodesOnSide = true; + } + if ( missedNodesOnSide ) + { + // clear uv_grid where nodes are missing + for ( size_t iQ = 0; iQ < quadVec.size(); ++iQ ) + quadVec[ iQ ]->uv_grid.clear(); + } + } + } + + return true; +} + +//================================================================================ +/*! + * \brief Splits a quad at I or J. Returns an index of a new side in the new quad + */ +//================================================================================ + +int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J) +{ + FaceQuadStruct* newQuad = new FaceQuadStruct( quad->face ); + myQuadList.push_back( FaceQuadStruct::Ptr( newQuad )); + + vector points; + if ( I > 0 && I <= quad->iSize-2 ) + { + points.reserve( quad->jSize ); + for ( int jP = 0; jP < quad->jSize; ++jP ) + points.push_back( quad->UVPt( I, jP )); + + newQuad->side.resize( 4 ); + newQuad->side[ QUAD_BOTTOM_SIDE ] = quad->side[ QUAD_BOTTOM_SIDE ]; + newQuad->side[ QUAD_RIGHT_SIDE ] = quad->side[ QUAD_RIGHT_SIDE ]; + newQuad->side[ QUAD_TOP_SIDE ] = quad->side[ QUAD_TOP_SIDE ]; + newQuad->side[ QUAD_LEFT_SIDE ] = StdMeshers_FaceSide::New( points, quad->face ); + + FaceQuadStruct::Side& newSide = newQuad->side[ QUAD_LEFT_SIDE ]; + FaceQuadStruct::Side& newSide2 = quad->side [ QUAD_RIGHT_SIDE ]; + + quad->side[ QUAD_RIGHT_SIDE ] = newSide; + + int iBot = quad->side[ QUAD_BOTTOM_SIDE ].ToSideIndex( I ); + int iTop = quad->side[ QUAD_TOP_SIDE ].ToSideIndex( I ); + + newSide.AddContact ( 0, & quad->side[ QUAD_BOTTOM_SIDE ], iBot ); + newSide2.AddContact( 0, & quad->side[ QUAD_BOTTOM_SIDE ], iBot ); + newSide.AddContact ( quad->jSize - 1, & quad->side[ QUAD_TOP_SIDE ], iTop ); + newSide2.AddContact( quad->jSize - 1, & quad->side[ QUAD_TOP_SIDE ], iTop ); + // cout << "Contact: L " << &newSide << " "<< newSide.NbPoints() + // << " R " << &newSide2 << " "<< newSide2.NbPoints() + // << " B " << &quad->side[ QUAD_BOTTOM_SIDE ] << " "<< quad->side[ QUAD_BOTTOM_SIDE].NbPoints() + // << " T " << &quad->side[ QUAD_TOP_SIDE ] << " "<< quad->side[ QUAD_TOP_SIDE].NbPoints()<< endl; + + newQuad->side[ QUAD_BOTTOM_SIDE ].from = iBot; + newQuad->side[ QUAD_TOP_SIDE ].from = iTop; + newQuad->name = ( TComm("Right of I=") << I ); + + 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 && 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 ) + points.push_back( quad->UVPt( iP, J )); + + newQuad->side.resize( 4 ); + newQuad->side[ QUAD_BOTTOM_SIDE ] = quad->side[ QUAD_BOTTOM_SIDE ]; + newQuad->side[ QUAD_RIGHT_SIDE ] = quad->side[ QUAD_RIGHT_SIDE ]; + newQuad->side[ QUAD_TOP_SIDE ] = StdMeshers_FaceSide::New( points, quad->face ); + newQuad->side[ QUAD_LEFT_SIDE ] = quad->side[ QUAD_LEFT_SIDE ]; + + FaceQuadStruct::Side& newSide = newQuad->side[ QUAD_TOP_SIDE ]; + FaceQuadStruct::Side& newSide2 = quad->side [ QUAD_BOTTOM_SIDE ]; + + quad->side[ QUAD_BOTTOM_SIDE ] = newSide; + + int iLft = quad->side[ QUAD_LEFT_SIDE ].ToSideIndex( J ); + int iRgt = quad->side[ QUAD_RIGHT_SIDE ].ToSideIndex( J ); + + newSide.AddContact ( 0, & quad->side[ QUAD_LEFT_SIDE ], iLft ); + newSide2.AddContact( 0, & quad->side[ QUAD_LEFT_SIDE ], iLft ); + newSide.AddContact ( quad->iSize - 1, & quad->side[ QUAD_RIGHT_SIDE ], iRgt ); + newSide2.AddContact( quad->iSize - 1, & quad->side[ QUAD_RIGHT_SIDE ], iRgt ); + // cout << "Contact: T " << &newSide << " "<< newSide.NbPoints() + // << " B " << &newSide2 << " "<< newSide2.NbPoints() + // << " 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; + + 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; + quad->side[ QUAD_LEFT_SIDE ].from = iLft; + quad->uv_grid.clear(); + + return QUAD_TOP_SIDE; + } + + myQuadList.pop_back(); + return -1; +} + +//================================================================================ +/*! + * \brief Updates UV of a side after moving its node + */ +//================================================================================ + +void StdMeshers_Quadrangle_2D::updateSideUV( FaceQuadStruct::Side& side, + int iForced, + const TQuadsBySide& quadsBySide, + int * iNext) +{ + if ( !iNext ) + { + side.forced_nodes.insert( iForced ); + + // update parts of the side before and after iForced + + set::iterator iIt = side.forced_nodes.upper_bound( iForced ); + int iEnd = Min( side.NbPoints()-1, ( iIt == side.forced_nodes.end() ) ? int(1e7) : *iIt ); + if ( iForced + 1 < iEnd ) + updateSideUV( side, iForced, quadsBySide, &iEnd ); + + iIt = side.forced_nodes.lower_bound( iForced ); + int iBeg = Max( 0, ( iIt == side.forced_nodes.begin() ) ? 0 : *--iIt ); + if ( iForced - 1 > iBeg ) + updateSideUV( side, iForced, quadsBySide, &iBeg ); + + return; + } + + 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 + + // from the quads get grid points adjacent to the side + // to make two sides of a temporary quad + vector< FaceQuadStruct::Ptr > quads = quadsBySide.find( side )->second; // copy! + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + points[ is2nd ].reserve( sideSize ); + size_t nbLoops = 0; + while ( points[is2nd].size() < sideSize ) + { + int iCur = iFrom + points[is2nd].size() - int( !points[is2nd].empty() ); + + // look for a quad adjacent to iCur-th point of the side + for ( size_t iQ = 0; iQ < quads.size(); ++iQ ) + { + FaceQuadStruct::Ptr q = quads[ iQ ]; + if ( !q ) + continue; + size_t iS; + 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 ); + else + isOut = ( q->side[ iS ].to >= iCur || q->side[ iS ].from <= iCur ); + if ( isOut ) + continue; + if ( !setNormalizedGrid( q )) + continue; + + // found - copy points + int i,j,di,dj,nb; + if ( iS % 2 ) // right or left + { + i = ( iS == QUAD_LEFT_SIDE ) ? 1 : q->iSize-2; + j = q->side[ iS ].ToQuadIndex( iCur ); + di = 0; + dj = ( q->side[ iS ].IsReversed() ) ? -1 : +1; + nb = ( q->side[ iS ].IsReversed() ) ? j+1 : q->jSize-j; + } + else // bottom or top + { + i = q->side[ iS ].ToQuadIndex( iCur ); + j = ( iS == QUAD_BOTTOM_SIDE ) ? 1 : q->jSize-2; + di = ( q->side[ iS ].IsReversed() ) ? -1 : +1; + dj = 0; + nb = ( q->side[ iS ].IsReversed() ) ? i+1 : q->iSize-i; + } + if ( !points[is2nd].empty() ) + { + gp_UV lastUV = points[is2nd].back().UV(); + gp_UV quadUV = q->UVPt( i, j ).UV(); + if ( ( lastUV - quadUV ).SquareModulus() > 1e-10 ) + continue; // quad is on the other side of the side + i += di; j += dj; --nb; + } + for ( ; nb > 0 ; --nb ) + { + points[ is2nd ].push_back( q->UVPt( i, j )); + if ( points[is2nd].size() >= sideSize ) + break; + i += di; j += dj; + } + quads[ iQ ].reset(); // not to use this quad anymore + + if ( points[is2nd].size() >= sideSize ) + break; + } // loop on quads + + if ( nbLoops++ > quads.size() ) + throw SALOME_Exception( "StdMeshers_Quadrangle_2D::updateSideUV() bug: infinite loop" ); + + } // while ( points[is2nd].size() < sideSize ) + } // two loops to fill points[0] and points[1] + + // points for other pair of opposite sides of the temporary quad + + enum { L,R,B,T }; // side index of points[] + + points[B].push_back( points[L].front() ); + points[B].push_back( side.GetUVPtStruct()[ iFrom ]); + points[B].push_back( points[R].front() ); + + points[T].push_back( points[L].back() ); + points[T].push_back( side.GetUVPtStruct()[ iTo-1 ]); + points[T].push_back( points[R].back() ); + + // make the temporary quad + FaceQuadStruct::Ptr tmpQuad + ( new FaceQuadStruct( TopoDS::Face( myHelper->GetSubShape() ), "tmpQuad")); + tmpQuad->side.push_back( StdMeshers_FaceSide::New( points[B] )); // bottom + tmpQuad->side.push_back( StdMeshers_FaceSide::New( points[R] )); // right + tmpQuad->side.push_back( StdMeshers_FaceSide::New( points[T] )); + tmpQuad->side.push_back( StdMeshers_FaceSide::New( points[L] )); + + // compute new UV of the side + setNormalizedGrid( tmpQuad ); + gp_UV uv = tmpQuad->UVPt(1,0).UV(); + tmpQuad->updateUV( uv, 1,0, /*isVertical=*/true ); + + // update UV of the side + vector& sidePoints = (vector&) side.GetUVPtStruct(); + for ( int i = iFrom; i < iTo; ++i ) + { + const uvPtStruct& uvPt = tmpQuad->UVPt( 1, i-iFrom ); + sidePoints[ i ].u = uvPt.u; + sidePoints[ i ].v = uvPt.v; + } +} + +//================================================================================ +/*! + * \brief Finds indices of a grid quad enclosing the given enforced UV + */ +//================================================================================ + +bool FaceQuadStruct::findCell( const gp_XY& UV, int & I, int & J ) +{ + // setNormalizedGrid() must be called before! + if ( uv_box.IsOut( UV )) + return false; + + // find an approximate position + double x = 0.5, y = 0.5; + gp_XY t0 = UVPt( iSize - 1, 0 ).UV(); + gp_XY t1 = UVPt( 0, jSize - 1 ).UV(); + gp_XY t2 = UVPt( 0, 0 ).UV(); + SMESH_MeshAlgos::GetBarycentricCoords( UV, t0, t1, t2, x, y ); + x = Min( 1., Max( 0., x )); + y = Min( 1., Max( 0., y )); + + // precise the position + normPa2IJ( x,y, I,J ); + if ( !isNear( UV, I,J )) + { + // look for the most close IJ by traversing uv_grid in the middle + double dist2, minDist2 = ( UV - UVPt( I,J ).UV() ).SquareModulus(); + for ( int isU = 0; isU < 2; ++isU ) + { + int ind1 = isU ? 0 : iSize / 2; + int ind2 = isU ? jSize / 2 : 0; + int di1 = isU ? Max( 2, iSize / 20 ) : 0; + int di2 = isU ? 0 : Max( 2, jSize / 20 ); + int i,nb = isU ? iSize / di1 : jSize / di2; + for ( i = 0; i < nb; ++i, ind1 += di1, ind2 += di2 ) + if (( dist2 = ( UV - UVPt( ind1,ind2 ).UV() ).SquareModulus() ) < minDist2 ) + { + I = ind1; + J = ind2; + if ( isNear( UV, I,J )) + return true; + minDist2 = ( UV - UVPt( I,J ).UV() ).SquareModulus(); + } + } + if ( !isNear( UV, I,J, Max( iSize, jSize ) /2 )) + return false; + } + return true; +} + +//================================================================================ +/*! + * \brief Find indices (i,j) of a point in uv_grid by normalized parameters (x,y) + */ +//================================================================================ + +void FaceQuadStruct::normPa2IJ(double X, double Y, int & I, int & J ) +{ + + I = Min( int ( iSize * X ), iSize - 2 ); + J = Min( int ( jSize * Y ), jSize - 2 ); + + int oldI, oldJ; + do + { + oldI = I, oldJ = J; + while ( X <= UVPt( I,J ).x && I != 0 ) + --I; + while ( X > UVPt( I+1,J ).x && I+2 < iSize ) + ++I; + while ( Y <= UVPt( I,J ).y && J != 0 ) + --J; + while ( Y > UVPt( I,J+1 ).y && J+2 < jSize ) + ++J; + } while ( oldI != I || oldJ != J ); +} + +//================================================================================ +/*! + * \brief Looks for UV in quads around a given (I,J) and precise (I,J) + */ +//================================================================================ + +bool FaceQuadStruct::isNear( const gp_XY& UV, int & I, int & J, int nbLoops ) +{ + if ( I+1 >= iSize ) I = iSize - 2; + if ( J+1 >= jSize ) J = jSize - 2; + + double bcI, bcJ; + gp_XY uvI, uvJ, uv0, uv1; + for ( int iLoop = 0; iLoop < nbLoops; ++iLoop ) + { + int oldI = I, oldJ = J; + + uvI = UVPt( I+1, J ).UV(); + uvJ = UVPt( I, J+1 ).UV(); + uv0 = UVPt( I, J ).UV(); + SMESH_MeshAlgos::GetBarycentricCoords( UV, uvI, uvJ, uv0, bcI, bcJ ); + if ( bcI >= 0. && bcJ >= 0. && bcI + bcJ <= 1.) + return true; + + if ( I > 0 && bcI < 0. ) --I; + if ( I+2 < iSize && bcI > 1. ) ++I; + if ( J > 0 && bcJ < 0. ) --J; + if ( J+2 < jSize && bcJ > 1. ) ++J; + + uv1 = UVPt( I+1,J+1).UV(); + if ( I != oldI || J != oldJ ) + { + uvI = UVPt( I+1, J ).UV(); + uvJ = UVPt( I, J+1 ).UV(); + } + SMESH_MeshAlgos::GetBarycentricCoords( UV, uvI, uvJ, uv1, bcI, bcJ ); + if ( bcI >= 0. && bcJ >= 0. && bcI + bcJ <= 1.) + return true; + + if ( I > 0 && bcI > 1. ) --I; + if ( I+2 < iSize && bcI < 0. ) ++I; + if ( J > 0 && bcJ > 1. ) --J; + if ( J+2 < jSize && bcJ < 0. ) ++J; + + if ( I == oldI && J == oldJ ) + return false; + + if ( iLoop+1 == nbLoops ) + { + uvI = UVPt( I+1, J ).UV(); + uvJ = UVPt( I, J+1 ).UV(); + uv0 = UVPt( I, J ).UV(); + SMESH_MeshAlgos::GetBarycentricCoords( UV, uvI, uvJ, uv0, bcI, bcJ ); + if ( bcI >= 0. && bcJ >= 0. && bcI + bcJ <= 1.) + return true; + + uv1 = UVPt( I+1,J+1).UV(); + SMESH_MeshAlgos::GetBarycentricCoords( UV, uvI, uvJ, uv1, bcI, bcJ ); + if ( bcI >= 0. && bcJ >= 0. && bcI + bcJ <= 1.) + return true; + } + } + return false; +} + +//================================================================================ +/*! + * \brief Checks if a given UV is equal to a given grid point + */ +//================================================================================ + +bool FaceQuadStruct::isEqual( const gp_XY& UV, int I, int J ) +{ + TopLoc_Location loc; + Handle(Geom_Surface) surf = BRep_Tool::Surface( face, loc ); + gp_Pnt p1 = surf->Value( UV.X(), UV.Y() ); + gp_Pnt p2 = surf->Value( UVPt( I,J ).u, UVPt( I,J ).v ); + + double dist2 = 1e100; + for ( int di = -1; di < 2; di += 2 ) + { + int i = I + di; + if ( i < 0 || i+1 >= iSize ) continue; + for ( int dj = -1; dj < 2; dj += 2 ) + { + int j = J + dj; + if ( j < 0 || j+1 >= jSize ) continue; + + dist2 = Min( dist2, + p2.SquareDistance( surf->Value( UVPt( i,j ).u, UVPt( i,j ).v ))); + } + } + double tol2 = dist2 / 1000.; + return p1.SquareDistance( p2 ) < tol2; +} + +//================================================================================ +/*! + * \brief Recompute UV of grid points around a moved point in one direction + */ +//================================================================================ + +void FaceQuadStruct::updateUV( const gp_XY& UV, int I, int J, bool isVertical ) +{ + UVPt( I, J ).u = UV.X(); + UVPt( I, J ).v = UV.Y(); + + if ( isVertical ) + { + // above J + if ( J+1 < jSize-1 ) + { + gp_UV a0 = UVPt( 0, J ).UV(); + gp_UV a1 = UVPt( iSize-1, J ).UV(); + gp_UV a2 = UVPt( iSize-1, jSize-1 ).UV(); + gp_UV a3 = UVPt( 0, jSize-1 ).UV(); + + gp_UV p0 = UVPt( I, J ).UV(); + gp_UV p2 = UVPt( I, jSize-1 ).UV(); + const double y0 = UVPt( I, J ).y, dy = 1. - y0; + for (int j = J+1; j < jSize-1; j++) + { + gp_UV p1 = UVPt( iSize-1, j ).UV(); + gp_UV p3 = UVPt( 0, j ).UV(); + + UVPtStruct& uvPt = UVPt( I, j ); + gp_UV uv = calcUV( uvPt.x, ( uvPt.y - y0 ) / dy, a0,a1,a2,a3, p0,p1,p2,p3); + uvPt.u = uv.X(); + uvPt.v = uv.Y(); + } + } + // under J + if ( J-1 > 0 ) + { + gp_UV a0 = UVPt( 0, 0 ).UV(); + gp_UV a1 = UVPt( iSize-1, 0 ).UV(); + gp_UV a2 = UVPt( iSize-1, J ).UV(); + gp_UV a3 = UVPt( 0, J ).UV(); + + gp_UV p0 = UVPt( I, 0 ).UV(); + gp_UV p2 = UVPt( I, J ).UV(); + const double y0 = 0., dy = UVPt( I, J ).y - y0; + for (int j = 1; j < J; j++) + { + gp_UV p1 = UVPt( iSize-1, j ).UV(); + gp_UV p3 = UVPt( 0, j ).UV(); + + UVPtStruct& uvPt = UVPt( I, j ); + gp_UV uv = calcUV( uvPt.x, ( uvPt.y - y0 ) / dy, a0,a1,a2,a3, p0,p1,p2,p3); + uvPt.u = uv.X(); + uvPt.v = uv.Y(); + } + } + } + else // horizontally + { + // before I + if ( I-1 > 0 ) + { + gp_UV a0 = UVPt( 0, 0 ).UV(); + gp_UV a1 = UVPt( I, 0 ).UV(); + gp_UV a2 = UVPt( I, jSize-1 ).UV(); + gp_UV a3 = UVPt( 0, jSize-1 ).UV(); + + gp_UV p1 = UVPt( I, J ).UV(); + gp_UV p3 = UVPt( 0, J ).UV(); + const double x0 = 0., dx = UVPt( I, J ).x - x0; + for (int i = 1; i < I; i++) + { + gp_UV p0 = UVPt( i, 0 ).UV(); + gp_UV p2 = UVPt( i, jSize-1 ).UV(); + + UVPtStruct& uvPt = UVPt( i, J ); + gp_UV uv = calcUV(( uvPt.x - x0 ) / dx , uvPt.y, a0,a1,a2,a3, p0,p1,p2,p3); + uvPt.u = uv.X(); + uvPt.v = uv.Y(); + } + } + // after I + if ( I+1 < iSize-1 ) + { + gp_UV a0 = UVPt( I, 0 ).UV(); + gp_UV a1 = UVPt( iSize-1, 0 ).UV(); + gp_UV a2 = UVPt( iSize-1, jSize-1 ).UV(); + gp_UV a3 = UVPt( I, jSize-1 ).UV(); + + gp_UV p1 = UVPt( iSize-1, J ).UV(); + gp_UV p3 = UVPt( I, J ).UV(); + const double x0 = UVPt( I, J ).x, dx = 1. - x0; + for (int i = I+1; i < iSize-1; i++) + { + gp_UV p0 = UVPt( i, 0 ).UV(); + gp_UV p2 = UVPt( i, jSize-1 ).UV(); + + UVPtStruct& uvPt = UVPt( i, J ); + gp_UV uv = calcUV(( uvPt.x - x0 ) / dx , uvPt.y, a0,a1,a2,a3, p0,p1,p2,p3); + uvPt.u = uv.X(); + uvPt.v = uv.Y(); + } + } + } +} + +//================================================================================ +/*! + * \brief Side copying + */ +//================================================================================ + +FaceQuadStruct::Side& FaceQuadStruct::Side::operator=(const Side& otherSide) +{ + grid = otherSide.grid; + from = otherSide.from; + to = otherSide.to; + di = otherSide.di; + forced_nodes = otherSide.forced_nodes; + contacts = otherSide.contacts; + nbNodeOut = otherSide.nbNodeOut; + + for ( size_t iC = 0; iC < contacts.size(); ++iC ) + { + FaceQuadStruct::Side* oSide = contacts[iC].other_side; + for ( size_t iOC = 0; iOC < oSide->contacts.size(); ++iOC ) + if ( oSide->contacts[iOC].other_side == & otherSide ) + { + // cout << "SHIFT old " << &otherSide << " " << otherSide.NbPoints() + // << " -> new " << this << " " << this->NbPoints() << endl; + oSide->contacts[iOC].other_side = this; + } + } + return *this; +} + +//================================================================================ +/*! + * \brief Converts node index of a quad to node index of this side + */ +//================================================================================ + +int FaceQuadStruct::Side::ToSideIndex( int quadNodeIndex ) const +{ + return from + di * quadNodeIndex; +} + +//================================================================================ +/*! + * \brief Converts node index of this side to node index of a quad + */ +//================================================================================ + +int FaceQuadStruct::Side::ToQuadIndex( int sideNodeIndex ) const +{ + return ( sideNodeIndex - from ) * di; +} + +//================================================================================ +/*! + * \brief Reverse the side + */ +//================================================================================ + +bool FaceQuadStruct::Side::Reverse(bool keepGrid) +{ + if ( grid ) + { + if ( keepGrid ) + { + from -= di; + to -= di; + std::swap( from, to ); + di *= -1; + } + else + { + grid->Reverse(); + } + } + return (bool)grid; +} + +//================================================================================ +/*! + * \brief Checks if a node is enforced + * \param [in] nodeIndex - an index of a node in a size + * \return bool - \c true if the node is forced + */ +//================================================================================ + +bool FaceQuadStruct::Side::IsForced( int nodeIndex ) const +{ + if ( nodeIndex < 0 || nodeIndex >= grid->NbPoints() ) + throw SALOME_Exception( " FaceQuadStruct::Side::IsForced(): wrong index" ); + + if ( forced_nodes.count( nodeIndex ) ) + return true; + + for ( size_t i = 0; i < this->contacts.size(); ++i ) + if ( contacts[ i ].point == nodeIndex && + contacts[ i ].other_side->forced_nodes.count( contacts[ i ].other_point )) + return true; + + return false; +} + +//================================================================================ +/*! + * \brief Sets up a contact between this and another side + */ +//================================================================================ + +void FaceQuadStruct::Side::AddContact( int ip, Side* side, int iop ) +{ + 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(); + c.point = ip; + c.other_side = side; + c.other_point = iop; + } + { + side->contacts.resize( side->contacts.size() + 1 ); + Contact& c = side->contacts.back(); + c.point = iop; + c.other_side = this; + c.other_point = ip; + } +} + +//================================================================================ +/*! + * \brief Returns a normalized parameter of a point indexed within a quadrangle + */ +//================================================================================ + +double FaceQuadStruct::Side::Param( int i ) const +{ + const vector& points = GetUVPtStruct(); + return (( points[ from + i * di ].normParam - points[ from ].normParam ) / + ( points[ to - 1 * di ].normParam - points[ from ].normParam )); +} + +//================================================================================ +/*! + * \brief Returns UV by a parameter normalized within a quadrangle + */ +//================================================================================ + +gp_XY FaceQuadStruct::Side::Value2d( double x ) const +{ + const vector& points = GetUVPtStruct(); + double u = ( points[ from ].normParam + + x * ( points[ to-di ].normParam - points[ from ].normParam )); + return grid->Value2d( u ).XY(); +} + +//================================================================================ +/*! + * \brief Returns side length + */ +//================================================================================ + +double FaceQuadStruct::Side::Length(int theFrom, int theTo) const +{ + if ( IsReversed() != ( theTo < theFrom )) + std::swap( theTo, theFrom ); + + const vector& points = GetUVPtStruct(); + double r; + if ( theFrom == theTo && theTo == -1 ) + r = Abs( First().normParam - + Last ().normParam ); + else if ( IsReversed() ) + r = Abs( points[ Max( to, theTo+1 ) ].normParam - + points[ Min( from, theFrom ) ].normParam ); + else + r = Abs( points[ Min( to, theTo-1 ) ].normParam - + points[ Max( from, theFrom ) ].normParam ); + return r * grid->Length(); }