X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=d5c4c9ad788c34b02eabf65666a4f478f16e3493;hb=61bbd63ee9b5f8b85237b8a45c8f454b766c25ca;hp=dbb505f9714dc5ab66f534b81e3f53c89c2ddcc8;hpb=3bc3b5cb5a6d85bcacffdc8334b87f4e56e46e6d;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index dbb505f97..d5c4c9ad7 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -1,23 +1,24 @@ -// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2012 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 +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS // -// 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. +// 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. // -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. // -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // File : StdMeshers_Quadrangle_2D.cxx // Author : Paul RASCLE, EDF @@ -25,21 +26,19 @@ #include "StdMeshers_Quadrangle_2D.hxx" -#include "StdMeshers_FaceSide.hxx" - -#include "StdMeshers_QuadrangleParams.hxx" - +#include "SMDS_EdgePosition.hxx" +#include "SMDS_FacePosition.hxx" +#include "SMDS_MeshElement.hxx" +#include "SMDS_MeshNode.hxx" +#include "SMESH_Block.hxx" +#include "SMESH_Comment.hxx" #include "SMESH_Gen.hxx" #include "SMESH_Mesh.hxx" -#include "SMESH_subMesh.hxx" #include "SMESH_MesherHelper.hxx" -#include "SMESH_Block.hxx" -#include "SMESH_Comment.hxx" - -#include "SMDS_MeshElement.hxx" -#include "SMDS_MeshNode.hxx" -#include "SMDS_EdgePosition.hxx" -#include "SMDS_FacePosition.hxx" +#include "SMESH_subMesh.hxx" +#include "StdMeshers_FaceSide.hxx" +#include "StdMeshers_QuadrangleParams.hxx" +#include "StdMeshers_ViscousLayers2D.hxx" #include #include @@ -78,7 +77,8 @@ typedef SMESH_Comment TComm; StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMESH_Gen* gen) - : SMESH_2D_Algo(hypId, studyId, gen) + : SMESH_2D_Algo(hypId, studyId, gen), + myHelper( 0 ) { MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D"); _name = "Quadrangle_2D"; @@ -86,7 +86,7 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, _compatibleHypothesis.push_back("QuadrangleParams"); _compatibleHypothesis.push_back("QuadranglePreference"); _compatibleHypothesis.push_back("TrianglePreference"); - myTool = 0; + _compatibleHypothesis.push_back("ViscousLayers2D"); } //============================================================================= @@ -192,19 +192,24 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis */ //============================================================================= -bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape)// throw (SALOME_Exception) +bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape) { - // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside - //Unexpect aCatchSalomeException); + const TopoDS_Face& F = TopoDS::Face(aShape); + Handle(Geom_Surface) S = BRep_Tool::Surface(F); SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); aMesh.GetSubMesh(aShape); SMESH_MesherHelper helper (aMesh); - myTool = &helper; + myHelper = &helper; + + myProxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F ); + if ( !myProxyMesh ) + return false; - _quadraticMesh = myTool->IsQuadraticSubMesh(aShape); + _quadraticMesh = myHelper->IsQuadraticSubMesh(aShape); + myNeedSmooth = false; FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape); std::auto_ptr quadDeleter (quad); // to delete quad at exit from Compute() @@ -222,6 +227,8 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, 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; } } @@ -237,6 +244,8 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, 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; } } @@ -257,9 +266,6 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, int nbhoriz = Min(nbdown, nbup); int nbvertic = Min(nbright, nbleft); - const TopoDS_Face& F = TopoDS::Face(aShape); - Handle(Geom_Surface) S = BRep_Tool::Surface(F); - // internal mesh nodes int i, j, geomFaceID = meshDS->ShapeToIndex(F); for (i = 1; i < nbhoriz - 1; i++) { @@ -301,11 +307,11 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, 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 = myTool->AddFace(a, b, c, d); + d = quad->uv_grid[(j + 1) * nbhoriz + i ].node; + SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d); if (face) { meshDS->SetMeshElementOnShape(face, geomFaceID); } @@ -382,7 +388,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = myTool->AddFace(a, b, c); + SMDS_MeshFace* face = myHelper->AddFace(a, b, c); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } else { // make quadrangle @@ -393,7 +399,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); if (!myTrianglePreference){ - SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } else { @@ -408,7 +414,7 @@ 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 = myTool->AddFace(a, c, d); + SMDS_MeshFace* face = myHelper->AddFace(a, c, d); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } } @@ -470,7 +476,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = myTool->AddFace(a, b, c); + SMDS_MeshFace* face = myHelper->AddFace(a, b, c); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } else { // make quadrangle @@ -480,7 +486,7 @@ 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 = myTool->AddFace(a, b, c, d); + SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } else { @@ -494,7 +500,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, d = uv_e1[nbright - 2].node; else d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node; - SMDS_MeshFace* face = myTool->AddFace(a, c, d); + SMDS_MeshFace* face = myHelper->AddFace(a, c, d); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } } @@ -542,7 +548,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = myTool->AddFace(a, b, c); + SMDS_MeshFace* face = myHelper->AddFace(a, b, c); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } else { // make quadrangle @@ -553,7 +559,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); if (!myTrianglePreference){ - SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } else { @@ -567,7 +573,7 @@ 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 = myTool->AddFace(a, c, d); + SMDS_MeshFace* face = myHelper->AddFace(a, c, d); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } } @@ -612,7 +618,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = myTool->AddFace(a, b, c); + SMDS_MeshFace* face = myHelper->AddFace(a, b, c); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } else { // make quadrangle @@ -622,7 +628,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, d = quad->uv_grid[nbhoriz*(near + 1) + 1].node; //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); if (!myTrianglePreference){ - SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } else { @@ -636,7 +642,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, d = uv_e2[1].node; else d = quad->uv_grid[nbhoriz*(k + 1) + 1].node; - SMDS_MeshFace* face = myTool->AddFace(a, c, d); + SMDS_MeshFace* face = myHelper->AddFace(a, c, d); if (face) meshDS->SetMeshElementOnShape(face, geomFaceID); } } @@ -646,6 +652,9 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } } + if ( myNeedSmooth ) + Smooth( quad ); + bool isOk = true; return isOk; } @@ -822,9 +831,16 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes } if (!E1.IsNull() && !E2.IsNull() && !E3.IsNull()) { - quad->side.push_back(new StdMeshers_FaceSide(F, E1, &aMesh, true, ignoreMediumNodes)); - quad->side.push_back(new StdMeshers_FaceSide(F, E2, &aMesh, true, ignoreMediumNodes)); - quad->side.push_back(new StdMeshers_FaceSide(F, E3, &aMesh, false,ignoreMediumNodes)); + if ( myProxyMesh->GetProxySubMesh( E1 ) || + myProxyMesh->GetProxySubMesh( E2 ) || + myProxyMesh->GetProxySubMesh( E3 ) ) + + 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); @@ -844,13 +860,16 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes delete quad; return quad = 0; } - else if (nbEdgesInWire.front() == 4) { // exactly 4 edges + else if (nbEdgesInWire.front() == 4) // exactly 4 edges + { for (; edgeIt != edges.end(); ++edgeIt, nbSides++) - quad->side.push_back(new StdMeshers_FaceSide(F, *edgeIt, &aMesh, - nbSidesside.push_back(new StdMeshers_FaceSide(F, *edgeIt, &aMesh, nbSides < TOP_SIDE, + ignoreMediumNodes, myProxyMesh)); } - else if (nbEdgesInWire.front() > 4) { // more than 4 edges - try to unite some + 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() @@ -868,10 +887,30 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes sideEdges.splice(sideEdges.begin(), edges, --edges.end()); } } - quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, - nbSidesside.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, nbSides < TOP_SIDE, + ignoreMediumNodes, myProxyMesh)); ++nbSides; } + if ( !degenSides.empty() && nbSides - degenSides.size() == 4 ) + { + myNeedSmooth = true; + for ( unsigned i = TOP_SIDE; i < quad->side.size(); ++i ) + quad->side[i]->Reverse(); + + for ( int i = degenSides.size()-1; i > -1; --i ) + { + StdMeshers_FaceSide* degenSide = quad->side[ degenSides[ i ]]; + delete degenSide; + quad->side.erase( quad->side.begin() + degenSides[ i ] ); + } + for ( unsigned i = TOP_SIDE; i < quad->side.size(); ++i ) + quad->side[i]->Reverse(); + + nbSides -= degenSides.size(); + } // issue 20222. Try to unite only edges shared by two same faces if (nbSides < 4) { // delete found sides @@ -903,7 +942,8 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes } } quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, - nbSidesside[i]->NbEdges(); ++e) - MESSAGE (myTool->GetMeshDS()->ShapeToIndex(quad->side[i]->Edge(e)) << " "); + MESSAGE (myHelper->GetMeshDS()->ShapeToIndex(quad->side[i]->Edge(e)) << " "); MESSAGE (")\n"); } //cout << endl; @@ -1230,37 +1270,32 @@ bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, //return error("Can't find nodes on sides"); return error(COMPERR_BAD_INPUT_MESH); - // nodes Id on "in" edges - if (! quad->isEdgeOut[0]) { - int j = 0; - for (int i = 0; i < nbhoriz; i++) { // down - int ij = j * nbhoriz + i; - uv_grid[ij].node = uv_e0[i].node; - } + if ( myNeedSmooth ) + UpdateDegenUV( quad ); + + // 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]; } - if (! quad->isEdgeOut[1]) { - int i = nbhoriz - 1; - for (int j = 0; j < nbvertic; j++) { // right - int ij = j * nbhoriz + i; - uv_grid[ij].node = uv_e1[j].node; - } + /*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]; } - if (! quad->isEdgeOut[2]) { - int j = nbvertic - 1; - for (int i = 0; i < nbhoriz; i++) { // up - int ij = j * nbhoriz + i; - uv_grid[ij].node = uv_e2[i].node; - } + /*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]; } - if (! quad->isEdgeOut[3]) { + /*if (! quad->isEdgeOut[3])*/ { int i = 0; - for (int j = 0; j < nbvertic; j++) { // left - int ij = j * nbhoriz + i; - uv_grid[ij].node = uv_e3[j].node; - } + for (int j = 0; j < nbvertic; j++) // left + uv_grid[ j * nbhoriz + i ] = uv_e3[j]; } - // normalized 2d values on grid + // normalized 2d parameters on grid for (int i = 0; i < nbhoriz; i++) { for (int j = 0; j < nbvertic; j++) { int ij = j * nbhoriz + i; @@ -1281,26 +1316,23 @@ bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, } // 4 --- projection on 2d domain (u,v) - gp_UV a0(uv_e0.front().u, uv_e0.front().v); - gp_UV a1(uv_e0.back().u, uv_e0.back().v); - gp_UV a2(uv_e2.back().u, uv_e2.back().v); - gp_UV a3(uv_e2.front().u, uv_e2.front().v); + gp_UV a0 (uv_e0.front().u, uv_e0.front().v); + gp_UV a1 (uv_e0.back().u, uv_e0.back().v ); + gp_UV a2 (uv_e2.back().u, uv_e2.back().v ); + gp_UV a3 (uv_e2.front().u, uv_e2.front().v); + + for (int i = 0; i < nbhoriz; 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 p1( uv_e1[j].u, uv_e1[j].v ); + gp_UV p3( uv_e3[j].u, uv_e3[j].v ); - for (int i = 0; i < nbhoriz; i++) { - for (int j = 0; j < nbvertic; j++) { int ij = j * nbhoriz + i; double x = uv_grid[ij].x; double y = uv_grid[ij].y; - double param_0 = uv_e0[0].normParam + x * (uv_e0.back().normParam - uv_e0[0].normParam); // sud - double param_2 = uv_e2[0].normParam + x * (uv_e2.back().normParam - uv_e2[0].normParam); // nord - double param_1 = uv_e1[0].normParam + y * (uv_e1.back().normParam - uv_e1[0].normParam); // est - double param_3 = uv_e3[0].normParam + y * (uv_e3.back().normParam - uv_e3[0].normParam); // ouest - - //MESSAGE("params "<side[0]->Value2d(param_0).XY(); - gp_UV p1 = quad->side[1]->Value2d(param_1).XY(); - gp_UV p2 = quad->side[2]->Value2d(param_2).XY(); - gp_UV p3 = quad->side[3]->Value2d(param_3).XY(); gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3); @@ -1318,7 +1350,8 @@ bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, static void ShiftQuad(FaceQuadStruct* quad, const int num, bool) { - StdMeshers_FaceSide* side[4] = { quad->side[0], quad->side[1], quad->side[2], quad->side[3] }; + StdMeshers_FaceSide* side[4] = { quad->side[0], quad->side[1], + quad->side[2], quad->side[3] }; for (int i = BOTTOM_SIDE; i < NB_SIDES; ++i) { int id = (i + num) % NB_SIDES; bool wasForward = (i < TOP_SIDE); @@ -1339,23 +1372,18 @@ static gp_UV CalcUV(double x0, double x1, double y0, double y1, const gp_UV& a0, const gp_UV& a1, const gp_UV& a2, const gp_UV& a3) { - 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); double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); double y = y0 + x * (y1 - y0); - double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam); - double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam); - double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam); - double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam); - - gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY(); - gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY(); - gp_UV p2 = quad->side[TOP_SIDE ]->Value2d(param_t).XY(); - gp_UV p3 = quad->side[LEFT_SIDE ]->Value2d(param_l).XY(); + gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(x).XY(); + gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(y).XY(); + gp_UV p2 = quad->side[TOP_SIDE ]->Value2d(x).XY(); + gp_UV p3 = quad->side[LEFT_SIDE ]->Value2d(y).XY(); gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3); @@ -1490,6 +1518,9 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl) return error(COMPERR_BAD_INPUT_MESH); + if ( myNeedSmooth ) + UpdateDegenUV( quad ); + // arrays for normalized params //cout<<"Dump B:"<AddFace(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 = - myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1), + 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); } @@ -1659,13 +1690,13 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, 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 = - myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1), + 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); } @@ -1744,13 +1775,13 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, for (j=1; jAddFace(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 = - myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1), + 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); } @@ -1782,13 +1813,13 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, for (i=1; iAddFace(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 = - myTool->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i,j+1), + 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); } @@ -1894,13 +1925,13 @@ 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 = - myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1), + 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); } @@ -1927,13 +1958,13 @@ bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, for (i=1; iAddFace(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 = - myTool->AddFace(NodesLast.Value(i,1), NodesLast.Value(i,2), + 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); } @@ -2100,20 +2131,197 @@ void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS, gp_Pnt d(theNode4->X(),theNode4->Y(),theNode4->Z()); SMDS_MeshFace* face; if (a.Distance(c) > b.Distance(d)){ - face = myTool->AddFace(theNode2, theNode4 , theNode1); + face = myHelper->AddFace(theNode2, theNode4 , theNode1); if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID); - face = myTool->AddFace(theNode2, theNode3, theNode4); + face = myHelper->AddFace(theNode2, theNode3, theNode4); if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID); } else{ - face = myTool->AddFace(theNode1, theNode2 ,theNode3); + face = myHelper->AddFace(theNode1, theNode2 ,theNode3); if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID); - face = myTool->AddFace(theNode1, theNode3, theNode4); + face = myHelper->AddFace(theNode1, theNode3, theNode4); if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID); } } +namespace +{ + enum uvPos { UV_A0, UV_A1, UV_A2, UV_A3, UV_B, UV_R, UV_T, UV_L, UV_SIZE }; + + inline SMDS_MeshNode* makeNode( UVPtStruct & uvPt, + const double y, + FaceQuadStruct* quad, + const gp_UV* UVs, + SMESH_MesherHelper* helper, + Handle(Geom_Surface) S) + { + const vector& uv_eb = quad->side[BOTTOM_SIDE]->GetUVPtStruct(); + const vector& uv_et = quad->side[TOP_SIDE ]->GetUVPtStruct(); + double rBot = ( uv_eb.size() - 1 ) * uvPt.normParam; + double rTop = ( uv_et.size() - 1 ) * uvPt.normParam; + int iBot = int( rBot ); + int iTop = int( rTop ); + double xBot = uv_eb[ iBot ].normParam + ( rBot - iBot ) * ( uv_eb[ iBot+1 ].normParam - uv_eb[ iBot ].normParam ); + double xTop = uv_et[ iTop ].normParam + ( rTop - iTop ) * ( uv_et[ iTop+1 ].normParam - uv_et[ iTop ].normParam ); + double x = xBot + y * ( xTop - xBot ); + + gp_UV uv = CalcUV(/*x,y=*/x, y, + /*a0,...=*/UVs[UV_A0], UVs[UV_A1], UVs[UV_A2], UVs[UV_A3], + /*p0=*/quad->side[BOTTOM_SIDE]->Value2d( x ).XY(), + /*p1=*/UVs[ UV_R ], + /*p2=*/quad->side[TOP_SIDE ]->Value2d( x ).XY(), + /*p3=*/UVs[ UV_L ]); + gp_Pnt P = S->Value( uv.X(), uv.Y() ); + uvPt.u = uv.X(); + uvPt.v = uv.Y(); + return helper->AddNode(P.X(), P.Y(), P.Z(), 0, uv.X(), uv.Y() ); + } + + void reduce42( const vector& curr_base, + vector& next_base, + const int j, + int & next_base_len, + FaceQuadStruct* quad, + gp_UV* UVs, + const double y, + SMESH_MesherHelper* helper, + Handle(Geom_Surface)& S) + { + // add one "HH": nodes a,b,c,d,e and faces 1,2,3,4,5,6 + // + // .-----a-----b i + 1 + // |\ 5 | 6 /| + // | \ | / | + // | c--d--e | + // |1 |2 |3 |4 | + // | | | | | + // .--.--.--.--. i + // + // j j+2 j+4 + + // a (i + 1, j + 2) + const SMDS_MeshNode*& Na = next_base[ ++next_base_len ].node; + if ( !Na ) + Na = makeNode( next_base[ next_base_len ], y, quad, UVs, helper, S ); + + // b (i + 1, j + 4) + const SMDS_MeshNode*& Nb = next_base[ ++next_base_len ].node; + if ( !Nb ) + Nb = makeNode( next_base[ next_base_len ], y, quad, UVs, helper, S ); + + // c + double u = (curr_base[j + 2].u + next_base[next_base_len - 2].u) / 2.0; + double v = (curr_base[j + 2].v + next_base[next_base_len - 2].v) / 2.0; + gp_Pnt P = S->Value(u,v); + SMDS_MeshNode* Nc = helper->AddNode(P.X(), P.Y(), P.Z(), 0, u, v); + + // d + u = (curr_base[j + 2].u + next_base[next_base_len - 1].u) / 2.0; + v = (curr_base[j + 2].v + next_base[next_base_len - 1].v) / 2.0; + P = S->Value(u,v); + SMDS_MeshNode* Nd = helper->AddNode(P.X(), P.Y(), P.Z(), 0, u, v); + + // e + u = (curr_base[j + 2].u + next_base[next_base_len].u) / 2.0; + v = (curr_base[j + 2].v + next_base[next_base_len].v) / 2.0; + P = S->Value(u,v); + SMDS_MeshNode* Ne = helper->AddNode(P.X(), P.Y(), P.Z(), 0, u, v); + + // Faces + helper->AddFace(curr_base[j + 0].node, + curr_base[j + 1].node, Nc, + next_base[next_base_len - 2].node); + + helper->AddFace(curr_base[j + 1].node, + curr_base[j + 2].node, Nd, Nc); + + helper->AddFace(curr_base[j + 2].node, + curr_base[j + 3].node, Ne, Nd); + + helper->AddFace(curr_base[j + 3].node, + curr_base[j + 4].node, Nb, Ne); + + helper->AddFace(Nc, Nd, Na, next_base[next_base_len - 2].node); + + helper->AddFace(Nd, Ne, Nb, Na); + } + + void reduce31( const vector& curr_base, + vector& next_base, + const int j, + int & next_base_len, + FaceQuadStruct* quad, + gp_UV* UVs, + const double y, + SMESH_MesherHelper* helper, + Handle(Geom_Surface)& S) + { + // add one "H": nodes b,c,e and faces 1,2,4,5 + // + // .---------b i + 1 + // |\ 5 /| + // | \ / | + // | c---e | + // |1 |2 |4 | + // | | | | + // .--.---.--. i + // + // j j+1 j+2 j+3 + + // b (i + 1, j + 3) + const SMDS_MeshNode*& Nb = next_base[ ++next_base_len ].node; + if ( !Nb ) + Nb = makeNode( next_base[ next_base_len ], y, quad, UVs, helper, S ); + + // c and e + double u1 = (curr_base[ j ].u + next_base[ next_base_len-1 ].u ) / 2.0; + double u2 = (curr_base[ j+3 ].u + next_base[ next_base_len ].u ) / 2.0; + double u3 = (u2 - u1) / 3.0; + // + double v1 = (curr_base[ j ].v + next_base[ next_base_len-1 ].v ) / 2.0; + double v2 = (curr_base[ j+3 ].v + next_base[ next_base_len ].v ) / 2.0; + double v3 = (v2 - v1) / 3.0; + // c + double u = u1 + u3; + double v = v1 + v3; + gp_Pnt P = S->Value(u,v); + SMDS_MeshNode* Nc = helper->AddNode( P.X(), P.Y(), P.Z(), 0, u, v ); + // e + u = u1 + u3 + u3; + v = v1 + v3 + v3; + P = S->Value(u,v); + SMDS_MeshNode* Ne = helper->AddNode( P.X(), P.Y(), P.Z(), 0, u, v ); + + // Faces + // 1 + helper->AddFace( curr_base[ j + 0 ].node, + curr_base[ j + 1 ].node, + Nc, + next_base[ next_base_len - 1 ].node); + // 2 + helper->AddFace( curr_base[ j + 1 ].node, + curr_base[ j + 2 ].node, Ne, Nc); + // 4 + helper->AddFace( curr_base[ j + 2 ].node, + curr_base[ j + 3 ].node, Nb, Ne); + // 5 + helper->AddFace(Nc, Ne, Nb, + next_base[ next_base_len - 1 ].node); + } + + typedef void (* PReduceFunction) ( const vector& curr_base, + vector& next_base, + const int j, + int & next_base_len, + FaceQuadStruct* quad, + gp_UV* UVs, + const double y, + SMESH_MesherHelper* helper, + Handle(Geom_Surface)& S); + +} // namespace + //======================================================================= /*! * Implementation of Reduced algorithm (meshing with quadrangles only) @@ -2133,7 +2341,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, int nt = quad->side[2]->NbPoints(); int nl = quad->side[3]->NbPoints(); - // Simple Reduce 8->6->4->2 (3 steps) Multiple Reduce 8->2 (1 step) + // Simple Reduce 10->8->6->4 (3 steps) Multiple Reduce 10->4 (1 step) // // .-----.-----.-----.-----. .-----.-----.-----.-----. // | / \ | / \ | | / \ | / \ | @@ -2154,7 +2362,6 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, int nb1 = nb; int nr1 = nr; int nt1 = nt; - int nl1 = nl; if (nr == nl) { if (nb < nt) { @@ -2163,7 +2370,6 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, } } else if (nb == nt) { - nl1 = nb; // and == nt nr1 = nb; // and == nt if (nl < nr) { nt1 = nl; @@ -2179,12 +2385,12 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, } // number of rows and columns - int nrows = nr1 - 1; // and also == nl1 - 1 + int nrows = nr1 - 1; int ncol_top = nt1 - 1; int ncol_bot = nb1 - 1; - // maximum number of bottom elements for "tree" simple reduce 3->1 - int max_tree31 = ncol_top * pow(3.0, nrows); - if (ncol_bot > max_tree31) + // 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 + if ( nrows < nrows_tree31 ) MultipleReduce = true; } @@ -2220,7 +2426,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, nl = quad->side[3]->NbPoints(); dh = abs(nb-nt); dv = abs(nr-nl); - int nbh = Max(nb,nt); + int nbh = Max(nb,nt); int nbv = Max(nr,nl); int addh = 0; int addv = 0; @@ -2242,6 +2448,9 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || 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++) { @@ -2291,7 +2500,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 nodes + // auxilary 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; @@ -2340,8 +2549,8 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, for (i=1; i<=dl; i++) { for (j=1; jAddFace(NodesL.Value(i,j), NodesL.Value(i+1,j), - NodesL.Value(i+1,j+1), NodesL.Value(i,j+1)); + 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); } } @@ -2395,8 +2604,8 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, for (i=1; i<=dr; i++) { for (j=1; jAddFace(NodesR.Value(i,j), NodesR.Value(i+1,j), - NodesR.Value(i+1,j+1), NodesR.Value(i,j+1)); + 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); } } @@ -2459,7 +2668,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (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); } @@ -2493,26 +2702,31 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, 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 int ncol_top = nt - 1; int ncol_bot = nb - 1; int npair_top = ncol_top / 2; // maximum number of bottom elements for "linear" simple reduce 4->2 - int max_lin = ncol_top + npair_top * 2 * nrows; - // maximum number of bottom elements for "linear" simple reduce 4->2 + int max_lin42 = ncol_top + npair_top * 2 * nrows; + // maximum number of bottom elements for "linear" simple reduce 3->1 int max_lin31 = ncol_top + ncol_top * 2 * nrows; // maximum number of bottom elements for "tree" simple reduce 4->2 - int max_tree42 = npair_top * pow(2.0, nrows + 1); - if (ncol_top > npair_top * 2) { - int delta = ncol_bot - max_tree42; - for (int irow = 1; irow < nrows; irow++) { - int nfour = delta / 4; - delta -= nfour * 2; + int max_tree42 = 0; + // number of rows needed to reduce ncol_bot to ncol_top using simple 4->2 "tree" + int nrows_tree42 = int( log( (double)(ncol_bot / ncol_top) )/log((double)2) ); // needed to avoid overflow at pow(2) while computing max_tree42 + if (nrows_tree42 < nrows) { + max_tree42 = npair_top * pow(2.0, nrows + 1); + if ( ncol_top > npair_top * 2 ) { + int delta = ncol_bot - max_tree42; + for (int irow = 1; irow < nrows; irow++) { + int nfour = delta / 4; + delta -= nfour * 2; + } + if (delta <= (ncol_top - npair_top * 2)) + max_tree42 = ncol_bot; } - if (delta <= (ncol_top - npair_top * 2)) - max_tree42 = ncol_bot; } // maximum number of bottom elements for "tree" simple reduce 3->1 //int max_tree31 = ncol_top * pow(3.0, nrows); @@ -2520,7 +2734,8 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, bool is_lin_42 = false; bool is_tree_31 = false; bool is_tree_42 = false; - if (ncol_bot > max_lin) { + int max_lin = max_lin42; + if (ncol_bot > max_lin42) { if (ncol_bot <= max_lin31) { is_lin_31 = true; max_lin = max_lin31; @@ -2556,53 +2771,184 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl) return error(COMPERR_BAD_INPUT_MESH); - // arrays for normalized params - TColStd_SequenceOfReal npb, npr, npt, npl; - for (j = 0; j < nb; j++) { - npb.Append(uv_eb[j].normParam); - } - for (i = 0; i < nr; i++) { - npr.Append(uv_er[i].normParam); - } - for (j = 0; j < nt; j++) { - npt.Append(uv_et[j].normParam); - } - for (i = 0; i < nl; i++) { - npl.Append(uv_el[i].normParam); - } + myHelper->SetElementsOnShape( true ); - // We will ajust new points to this grid - if (!SetNormalizedGrid(aMesh, aShape, quad)) - return false; + 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 ); + uv[ UV_A2 ].SetCoord( uv_et.back().u, uv_et.back().v ); + uv[ UV_A3 ].SetCoord( uv_et.front().u, uv_et.front().v); - // TODO ??? - 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); - //========================================================= + vector curr_base = uv_eb, next_base; - TColStd_SequenceOfInteger curr_base, next_base; - TColStd_SequenceOfReal curr_par_u, curr_par_v; - TColStd_SequenceOfReal next_par_u, next_par_v; - StdMeshers_Array2OfNode NodesBRD (1,nb, 1,nr); - for (j = 1; j <= nb; j++) { - NodesBRD.SetValue(j, 1, uv_eb[j - 1].node); // bottom - curr_base.Append(j); - next_base.Append(-1); - curr_par_u.Append(uv_eb[j-1].u); - curr_par_v.Append(uv_eb[j-1].v); - next_par_u.Append(0.); - next_par_v.Append(0.); - } - for (j = 1; j <= nt; j++) { - NodesBRD.SetValue(j, nr, uv_et[j - 1].node); // top - } + UVPtStruct nullUVPtStruct; nullUVPtStruct.node = 0; int curr_base_len = nb; int next_base_len = 0; - if (is_tree_42) { + if ( true ) + { // ------------------------------------------------------------------ + // New algorithm implemented by request of IPAL22856 + // "2D quadrangle mesher of reduced type works wrong" + // http://bugtracker.opencascade.com/show_bug.cgi?id=22856 + + // the algorithm is following: all reduces are centred in horizontal + // direction and are distributed among all rows + + if (ncol_bot > max_tree42) { + is_lin_31 = true; + } + else { + if ((ncol_top/3)*3 == ncol_top ) { + is_lin_31 = true; + } + else { + is_lin_42 = true; + } + } + + const int col_top_size = is_lin_42 ? 2 : 1; + const int col_base_size = is_lin_42 ? 4 : 3; + + // Compute nb of "columns" (like in "linear" simple reducing) in all rows + + vector nb_col_by_row; + + 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; + if (remainder > 0) { + nb_col++; + } + if ( nb_col * col_top_size >= nt ) // == "tree" reducing situation + { + // top row is full (all elements reduced), add "columns" one by one + // in rows below until all bottom elements are reduced + 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; + 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 ]; + 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; + for ( int irow = 0; irow < nbrows_not_full; ++irow ) + nb_col_by_row[ irow ] += nb_col_add; + nbrows_not_full --; + remainder -= nb_col_add * delta_one_col; + } + } + else // == "linear" reducing situation + { + nb_col_by_row.resize( nrows, nb_col ); + if (remainder > 0) + for ( int irow = remainder / 2; irow < nrows; ++irow ) + nb_col_by_row[ irow ]--; + } + + // Make elements + + PReduceFunction reduceFunction = & ( is_lin_42 ? reduce42 : reduce31 ); + + const int reduce_grp_size = is_lin_42 ? 4 : 3; + + for (i = 1; i < nr; i++) // layer by layer + { + nb_col = nb_col_by_row[ i-1 ]; + int nb_next = curr_base_len - nb_col * 2; + if (nb_next < nt) nb_next = nt; + + const double y = uv_el[ i ].normParam; + + if ( i + 1 == nr ) // top + { + next_base = uv_et; + } + else + { + next_base.clear(); + next_base.resize( nb_next, nullUVPtStruct ); + next_base.front() = uv_el[i]; + next_base.back() = uv_er[i]; + + // compute normalized param u + double du = 1. / ( nb_next - 1 ); + next_base[0].normParam = 0.; + for ( j = 1; j < nb_next; ++j ) + next_base[j].normParam = next_base[j-1].normParam + du; + } + uv[ UV_L ].SetCoord( next_base.front().u, next_base.front().v ); + uv[ UV_R ].SetCoord( next_base.back().u, next_base.back().v ); + + int free_left = ( curr_base_len - 1 - nb_col * col_base_size ) / 2; + int free_middle = curr_base_len - 1 - nb_col * col_base_size - 2 * free_left; + + // not reduced left elements + for (j = 0; j < free_left; j++) + { + // f (i + 1, j + 1) + const SMDS_MeshNode*& Nf = next_base[++next_base_len].node; + if ( !Nf ) + Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S ); + + myHelper->AddFace(curr_base[ j ].node, + curr_base[ j+1 ].node, + Nf, + next_base[ next_base_len-1 ].node); + } + + for (int icol = 1; icol <= nb_col; icol++) + { + // add "H" + reduceFunction( curr_base, next_base, j, next_base_len, quad, uv, y, myHelper, S ); + + j += reduce_grp_size; + + // elements in the middle of "columns" added for symmetry + if ( free_middle > 0 && ( nb_col % 2 == 0 ) && icol == nb_col / 2 ) + { + for (int imiddle = 1; imiddle <= free_middle; imiddle++) { + // f (i + 1, j + imiddle) + const SMDS_MeshNode*& Nf = next_base[++next_base_len].node; + if ( !Nf ) + Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S ); + + myHelper->AddFace(curr_base[ j-1+imiddle ].node, + curr_base[ j +imiddle ].node, + Nf, + next_base[ next_base_len-1 ].node); + } + j += free_middle; + } + } + + // not reduced right elements + for (; j < curr_base_len-1; j++) { + // f (i + 1, j + 1) + const SMDS_MeshNode*& Nf = next_base[++next_base_len].node; + if ( !Nf ) + Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S ); + + myHelper->AddFace(curr_base[ j ].node, + curr_base[ j+1 ].node, + Nf, + next_base[ next_base_len-1 ].node); + } + + curr_base_len = next_base_len + 1; + next_base_len = 0; + curr_base.swap( next_base ); + } + + } + else if ( is_tree_42 || is_tree_31 ) + { // "tree" simple reduce "42": 2->4->8->16->32->... // // .-------------------------------.-------------------------------. nr @@ -2624,236 +2970,6 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, // .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1 // 1 j nb - for (i = 1; i < nr; i++) { // layer by layer - // left - NodesBRD.SetValue(1, i+1, uv_el[i].node); - next_base.SetValue(++next_base_len, 1); - // right - NodesBRD.SetValue(nb, i+1, uv_er[i].node); - - next_par_u.SetValue(next_base_len, uv_el[i].u); - next_par_v.SetValue(next_base_len, uv_el[i].v); - - // to stop reducing, if number of nodes reaches nt - int delta = curr_base_len - nt; - - //double du = uv_er[i].u - uv_el[i].u; - //double dv = uv_er[i].v - uv_el[i].v; - - // to calculate normalized parameter, we must know number of points in next layer - int nb_four = (curr_base_len - 1) / 4; - int nb_next = nb_four*2 + (curr_base_len - nb_four*4); - if (nb_next < nt) nb_next = nt; - - for (j = 1; j + 4 <= curr_base_len && delta > 0; j += 4, delta -= 2) { - // add one "HH": nodes a,b,c,d,e and faces 1,2,3,4,5,6 - // - // .-----a-----b i + 1 - // |\ 5 | 6 /| - // | \ | / | - // | c--d--e | - // |1 |2 |3 |4 | - // | | | | | - // .--.--.--.--. i - // - // j j+2 j+4 - - double u,v; - - // a (i + 1, j + 2) - const SMDS_MeshNode* Na; - next_base_len++; - next_base.SetValue(next_base_len, curr_base.Value(j + 2)); - if (i + 1 == nr) { // top - Na = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else { - //double norm_par = double(next_base_len - 1)/double(nb_next - 1); - //u = uv_el[i].u + du * norm_par; - //v = uv_el[i].v + dv * norm_par; - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 2)); - //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 2)); - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Na1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Na1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na1); - Na = Na1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - - // b (i + 1, j + 4) - const SMDS_MeshNode* Nb; - next_base_len++; - next_base.SetValue(next_base_len, curr_base.Value(j + 4)); - if (i + 1 == nr) { // top - Nb = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 4 == curr_base_len) { // right - Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - //double norm_par = double(next_base_len - 1)/double(nb_next - 1); - //u = uv_el[i].u + du * norm_par; - //v = uv_el[i].v + dv * norm_par; - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 4)); - //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 4)); - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1); - Nb = Nb1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - - // c - u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 2)) / 2.0; - v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 2)) / 2.0; - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nc, geomFaceID, u, v); - - // d - u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 1)) / 2.0; - v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 1)) / 2.0; - P = S->Value(u,v); - SMDS_MeshNode* Nd = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nd, geomFaceID, u, v); - - // e - u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len)) / 2.0; - v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len)) / 2.0; - P = S->Value(u,v); - SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Ne, geomFaceID, u, v); - - // Faces - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i), - NodesBRD.Value(curr_base.Value(j + 1), i), - Nc, - NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); - - SMDS_MeshFace* F2 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i), - NodesBRD.Value(curr_base.Value(j + 2), i), - Nd, Nc); - if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID); - - SMDS_MeshFace* F3 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i), - NodesBRD.Value(curr_base.Value(j + 3), i), - Ne, Nd); - if (F3) meshDS->SetMeshElementOnShape(F3, geomFaceID); - - SMDS_MeshFace* F4 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 3), i), - NodesBRD.Value(curr_base.Value(j + 4), i), - Nb, Ne); - if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID); - - SMDS_MeshFace* F5 = myTool->AddFace(Nc, Nd, Na, - NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1)); - if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID); - - SMDS_MeshFace* F6 = myTool->AddFace(Nd, Ne, Nb, Na); - if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID); - } - - // not reduced side elements (if any) - for (; j < curr_base_len; j++) { - // f (i + 1, j + 1) - const SMDS_MeshNode* Nf; - double u,v; - next_base.SetValue(++next_base_len, curr_base.Value(j + 1)); - if (i + 1 == nr) { // top - Nf = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 1 == curr_base_len) { // right - Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - //double norm_par = double(next_base_len - 1)/double(nb_next - 1); - //u = uv_el[i].u + du * norm_par; - //v = uv_el[i].v + dv * norm_par; - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 1)); - //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 1)); - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1); - Nf = Nf1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i), - NodesBRD.Value(curr_base.Value(j + 1), i), - NodesBRD.Value(next_base.Value(next_base_len), i + 1), - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); - } - - curr_base_len = next_base_len; - curr_base = next_base; - curr_par_u = next_par_u; - curr_par_v = next_par_v; - next_base_len = 0; - } - } // end "tree" simple reduce "42" - else if (is_tree_31) { // "tree" simple reduce "31": 1->3->9->27->... // // .-----------------------------------------------------. nr @@ -2871,178 +2987,84 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, // .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1 // 1 j nb - for (i = 1; i < nr; i++) { // layer by layer - // left - NodesBRD.SetValue(1, i+1, uv_el[i].node); - next_base.SetValue(++next_base_len, 1); - // right - NodesBRD.SetValue(nb, i+1, uv_er[i].node); + PReduceFunction reduceFunction = & ( is_tree_42 ? reduce42 : reduce31 ); - next_par_u.SetValue(next_base_len, uv_el[i].u); - next_par_v.SetValue(next_base_len, uv_el[i].v); + const int reduce_grp_size = is_tree_42 ? 4 : 3; + for (i = 1; i < nr; i++) // layer by layer + { // to stop reducing, if number of nodes reaches nt int delta = curr_base_len - nt; // to calculate normalized parameter, we must know number of points in next layer - int nb_three = (curr_base_len - 1) / 3; - int nb_next = nb_three + (curr_base_len - nb_three*3); + int nb_reduce_groups = (curr_base_len - 1) / reduce_grp_size; + int nb_next = nb_reduce_groups * (reduce_grp_size-2) + (curr_base_len - nb_reduce_groups*reduce_grp_size); if (nb_next < nt) nb_next = nt; - for (j = 1; j + 3 <= curr_base_len && delta > 0; j += 3, delta -= 2) { - // add one "H": nodes b,c,e and faces 1,2,4,5 - // - // .---------b i + 1 - // |\ 5 /| - // | \ / | - // | c---e | - // |1 |2 |4 | - // | | | | - // .--.---.--. i - // - // j j+1 j+2 j+3 - - double u,v; - - // b (i + 1, j + 3) - const SMDS_MeshNode* Nb; - next_base_len++; - next_base.SetValue(next_base_len, curr_base.Value(j + 3)); - if (i + 1 == nr) { // top - Nb = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 3 == curr_base_len) { // right - Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1); - Nb = Nb1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - - // c and e - double u1 = (curr_par_u.Value(j) + next_par_u.Value(next_base_len - 1)) / 2.0; - double u2 = (curr_par_u.Value(j + 3) + next_par_u.Value(next_base_len)) / 2.0; - double u3 = (u2 - u1) / 3.0; - - double v1 = (curr_par_v.Value(j) + next_par_v.Value(next_base_len - 1)) / 2.0; - double v2 = (curr_par_v.Value(j + 3) + next_par_v.Value(next_base_len)) / 2.0; - double v3 = (v2 - v1) / 3.0; - - // c - u = u1 + u3; - v = v1 + v3; - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nc, geomFaceID, u, v); - - // e - u = u1 + u3 + u3; - v = v1 + v3 + v3; - P = S->Value(u,v); - SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Ne, geomFaceID, u, v); - - // Faces - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i), - NodesBRD.Value(curr_base.Value(j + 1), i), - Nc, - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); - - SMDS_MeshFace* F2 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i), - NodesBRD.Value(curr_base.Value(j + 2), i), - Ne, Nc); - if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID); - - SMDS_MeshFace* F4 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i), - NodesBRD.Value(curr_base.Value(j + 3), i), - Nb, Ne); - if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID); - - SMDS_MeshFace* F5 = myTool->AddFace(Nc, Ne, Nb, - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID); + const double y = uv_el[ i ].normParam; + + if ( i + 1 == nr ) // top + { + next_base = uv_et; + } + else + { + next_base.clear(); + next_base.resize( nb_next, nullUVPtStruct ); + next_base.front() = uv_el[i]; + next_base.back() = uv_er[i]; + + // compute normalized param u + double du = 1. / ( nb_next - 1 ); + next_base[0].normParam = 0.; + for ( j = 1; j < nb_next; ++j ) + next_base[j].normParam = next_base[j-1].normParam + du; + } + uv[ UV_L ].SetCoord( next_base.front().u, next_base.front().v ); + uv[ UV_R ].SetCoord( next_base.back().u, next_base.back().v ); + + for (j = 0; j+reduce_grp_size < curr_base_len && delta > 0; j+=reduce_grp_size, delta-=2) + { + reduceFunction( curr_base, next_base, j, next_base_len, quad, uv, y, myHelper, S ); } // not reduced side elements (if any) - for (; j < curr_base_len; j++) { + for (; j < curr_base_len-1; j++) + { // f (i + 1, j + 1) - const SMDS_MeshNode* Nf; - double u,v; - next_base.SetValue(++next_base_len, curr_base.Value(j + 1)); - if (i + 1 == nr) { // top - Nf = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 1 == curr_base_len) { // right - Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1); - Nf = Nf1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i), - NodesBRD.Value(curr_base.Value(j + 1), i), - NodesBRD.Value(next_base.Value(next_base_len), i + 1), - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); - } - - curr_base_len = next_base_len; - curr_base = next_base; - curr_par_u = next_par_u; - curr_par_v = next_par_v; + const SMDS_MeshNode*& Nf = next_base[++next_base_len].node; + if ( !Nf ) + Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S ); + + myHelper->AddFace(curr_base[ j ].node, + curr_base[ j+1 ].node, + Nf, + next_base[ next_base_len-1 ].node); + } + curr_base_len = next_base_len + 1; next_base_len = 0; + curr_base.swap( next_base ); } - } // end "tree" simple reduce "31" - else if (is_lin_42) { + } // end "tree" simple reduce + + else if ( is_lin_42 || is_lin_31 ) { + // "linear" simple reduce "31": 2->6->10->14 + // + // .-----------------------------.-----------------------------. nr + // | \ / | \ / | + // | .---------. | .---------. | + // | | | | | | | + // .---------.---------.---------.---------.---------.---------. + // | / \ / \ | / \ / \ | + // | / .-----. \ | / .-----. \ | i + // | / | | \ | / | | \ | + // .-----.-----.-----.-----.-----.-----.-----.-----.-----.-----. + // | / / \ / \ \ | / / \ / \ \ | + // | / / .-. \ \ | / / .-. \ \ | + // | / / / \ \ \ | / / / \ \ \ | + // .--.----.---.-----.---.-----.-.--.----.---.-----.---.-----.-. 1 + // 1 j nb + // "linear" simple reduce "42": 4->8->12->16 // // .---------------.---------------.---------------.---------------. nr @@ -3078,7 +3100,8 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, if (remainder > 0) { nb_col++; } - int free_left = ((nt - 1) - nb_col * 2) / 2; + const int col_top_size = is_lin_42 ? 2 : 1; + int free_left = ((nt - 1) - nb_col * col_top_size) / 2; free_left += nr - 2; int free_middle = (nr - 2) * 2; if (remainder > 0 && nb_col == 1) { @@ -3095,17 +3118,12 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, //int free_left = 2; //int free_middle = 4; - for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) { // layer by layer - // left - NodesBRD.SetValue(1, i+1, uv_el[i].node); - next_base.SetValue(++next_base_len, 1); - // right - NodesBRD.SetValue(nb, i+1, uv_er[i].node); + PReduceFunction reduceFunction = & ( is_lin_42 ? reduce42 : reduce31 ); - // left - next_par_u.SetValue(next_base_len, uv_el[i].u); - next_par_v.SetValue(next_base_len, uv_el[i].v); + const int reduce_grp_size = is_lin_42 ? 4 : 3; + for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) // layer by layer + { // to calculate normalized parameter, we must know number of points in next layer int nb_next = curr_base_len - nb_col * 2; if (remainder > 0 && i > remainder / 2) @@ -3113,46 +3131,40 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, nb_next += 2; if (nb_next < nt) nb_next = nt; + const double y = uv_el[ i ].normParam; + + if ( i + 1 == nr ) // top + { + next_base = uv_et; + } + else + { + next_base.clear(); + next_base.resize( nb_next, nullUVPtStruct ); + next_base.front() = uv_el[i]; + next_base.back() = uv_er[i]; + + // compute normalized param u + double du = 1. / ( nb_next - 1 ); + next_base[0].normParam = 0.; + for ( j = 1; j < nb_next; ++j ) + next_base[j].normParam = next_base[j-1].normParam + du; + } + uv[ UV_L ].SetCoord( next_base.front().u, next_base.front().v ); + uv[ UV_R ].SetCoord( next_base.back().u, next_base.back().v ); + // not reduced left elements - for (j = 1; j <= free_left; j++) { + for (j = 0; j < free_left; j++) + { // f (i + 1, j + 1) - const SMDS_MeshNode* Nf; - double u,v; - next_base.SetValue(++next_base_len, curr_base.Value(j + 1)); - if (i + 1 == nr) { // top - Nf = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1); - Nf = Nf1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i), - NodesBRD.Value(curr_base.Value(j + 1), i), - NodesBRD.Value(next_base.Value(next_base_len), i + 1), - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); + const SMDS_MeshNode*& Nf = next_base[++next_base_len].node; + if ( !Nf ) + Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S ); + + myHelper->AddFace(curr_base[ j ].node, + curr_base[ j+1 ].node, + Nf, + next_base[ next_base_len-1 ].node); } for (int icol = 1; icol <= nb_col; icol++) { @@ -3161,144 +3173,10 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, // stop short "column" break; - // add one "HH": nodes a,b,c,d,e and faces 1,2,3,4,5,6 - // - // .-----a-----b i + 1 - // |\ 5 | 6 /| - // | \ | / | - // | c--d--e | - // |1 |2 |3 |4 | - // | | | | | - // .--.--.--.--. i - // - // j j+2 j+4 - - double u,v; - - // a (i + 1, j + 2) - const SMDS_MeshNode* Na; - next_base_len++; - next_base.SetValue(next_base_len, curr_base.Value(j + 2)); - if (i + 1 == nr) { // top - Na = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Na1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Na1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na1); - Na = Na1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - - // b (i + 1, j + 4) - const SMDS_MeshNode* Nb; - next_base_len++; - next_base.SetValue(next_base_len, curr_base.Value(j + 4)); - if (i + 1 == nr) { // top - Nb = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 4 == curr_base_len) { // right - Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1); - Nb = Nb1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - - // c - u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 2)) / 2.0; - v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 2)) / 2.0; - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nc, geomFaceID, u, v); - - // d - u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 1)) / 2.0; - v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 1)) / 2.0; - P = S->Value(u,v); - SMDS_MeshNode* Nd = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nd, geomFaceID, u, v); - - // e - u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len)) / 2.0; - v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len)) / 2.0; - P = S->Value(u,v); - SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Ne, geomFaceID, u, v); - - // Faces - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i), - NodesBRD.Value(curr_base.Value(j + 1), i), - Nc, - NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); - - SMDS_MeshFace* F2 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i), - NodesBRD.Value(curr_base.Value(j + 2), i), - Nd, Nc); - if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID); - - SMDS_MeshFace* F3 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i), - NodesBRD.Value(curr_base.Value(j + 3), i), - Ne, Nd); - if (F3) meshDS->SetMeshElementOnShape(F3, geomFaceID); - - SMDS_MeshFace* F4 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 3), i), - NodesBRD.Value(curr_base.Value(j + 4), i), - Nb, Ne); - if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID); - - SMDS_MeshFace* F5 = myTool->AddFace(Nc, Nd, Na, - NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1)); - if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID); - - SMDS_MeshFace* F6 = myTool->AddFace(Nd, Ne, Nb, Na); - if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID); - - j += 4; + // add "H" + reduceFunction( curr_base, next_base, j, next_base_len, quad, uv, y, myHelper, S ); + + j += reduce_grp_size; // not reduced middle elements if (icol < nb_col) { @@ -3313,424 +3191,256 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, for (int imiddle = 1; imiddle <= free_add; imiddle++) { // f (i + 1, j + imiddle) - const SMDS_MeshNode* Nf; - double u,v; - next_base.SetValue(++next_base_len, curr_base.Value(j + imiddle)); - if (i + 1 == nr) { // top - Nf = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + imiddle == curr_base_len) { // right - Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1); - Nf = Nf1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j - 1 + imiddle), i), - NodesBRD.Value(curr_base.Value(j + imiddle), i), - NodesBRD.Value(next_base.Value(next_base_len), i + 1), - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); + const SMDS_MeshNode*& Nf = next_base[++next_base_len].node; + if ( !Nf ) + Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S ); + + myHelper->AddFace(curr_base[ j-1+imiddle ].node, + curr_base[ j +imiddle ].node, + Nf, + next_base[ next_base_len-1 ].node); } j += free_add; } } // not reduced right elements - for (; j < curr_base_len; j++) { + for (; j < curr_base_len-1; j++) { // f (i + 1, j + 1) - const SMDS_MeshNode* Nf; - double u,v; - next_base.SetValue(++next_base_len, curr_base.Value(j + 1)); - if (i + 1 == nr) { // top - Nf = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 1 == curr_base_len) { // right - Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1); - Nf = Nf1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i), - NodesBRD.Value(curr_base.Value(j + 1), i), - NodesBRD.Value(next_base.Value(next_base_len), i + 1), - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); - } - - curr_base_len = next_base_len; - curr_base = next_base; - curr_par_u = next_par_u; - curr_par_v = next_par_v; + const SMDS_MeshNode*& Nf = next_base[++next_base_len].node; + if ( !Nf ) + Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S ); + + myHelper->AddFace(curr_base[ j ].node, + curr_base[ j+1 ].node, + Nf, + next_base[ next_base_len-1 ].node); + } + + curr_base_len = next_base_len + 1; next_base_len = 0; + curr_base.swap( next_base ); } - } // end "linear" simple reduce "42" - else if (is_lin_31) { - // "linear" simple reduce "31": 2->6->10->14 - // - // .-----------------------------.-----------------------------. nr - // | \ / | \ / | - // | .---------. | .---------. | - // | | | | | | | - // .---------.---------.---------.---------.---------.---------. - // | / \ / \ | / \ / \ | - // | / .-----. \ | / .-----. \ | i - // | / | | \ | / | | \ | - // .-----.-----.-----.-----.-----.-----.-----.-----.-----.-----. - // | / / \ / \ \ | / / \ / \ \ | - // | / / .-. \ \ | / / .-. \ \ | - // | / / / \ \ \ | / / / \ \ \ | - // .--.----.---.-----.---.-----.-.--.----.---.-----.---.-----.-. 1 - // 1 j nb - int delta_all = nb - nt; - int delta_one_col = (nr - 1) * 2; - int nb_col = delta_all / delta_one_col; - int remainder = delta_all - nb_col * delta_one_col; - if (remainder > 0) { - nb_col++; - } - int free_left = ((nt - 1) - nb_col) / 2; - free_left += nr - 2; - int free_middle = (nr - 2) * 2; - if (remainder > 0 && nb_col == 1) { - int nb_rows_short_col = remainder / 2; - int nb_rows_thrown = (nr - 1) - nb_rows_short_col; - free_left -= nb_rows_thrown; - } + } // end "linear" simple reduce - for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) { // layer by layer - // left - NodesBRD.SetValue(1, i+1, uv_el[i].node); - next_base.SetValue(++next_base_len, 1); - // right - NodesBRD.SetValue(nb, i+1, uv_er[i].node); + else { + return false; + } + } // end Simple Reduce implementation - // left - next_par_u.SetValue(next_base_len, uv_el[i].u); - next_par_v.SetValue(next_base_len, uv_el[i].v); + bool isOk = true; + return isOk; +} - // to calculate normalized parameter, we must know number of points in next layer - int nb_next = curr_base_len - nb_col * 2; - if (remainder > 0 && i > remainder / 2) - // take into account short "column" - nb_next += 2; - if (nb_next < nt) nb_next = nt; +//================================================================================ +namespace // data for smoothing +{ + struct TSmoothNode; + // -------------------------------------------------------------------------------- + /*! + * \brief Structure used to check validity of node position after smoothing. + * It holds two nodes connected to a smoothed node and belonging to + * one mesh face + */ + struct TTriangle + { + TSmoothNode* _n1; + TSmoothNode* _n2; + TTriangle( TSmoothNode* n1=0, TSmoothNode* n2=0 ): _n1(n1), _n2(n2) {} + + inline bool IsForward( gp_UV uv ) const; + }; + // -------------------------------------------------------------------------------- + /*! + * \brief Data of a smoothed node + */ + struct TSmoothNode + { + gp_XY _uv; + vector< TTriangle > _triangles; // if empty, then node is not movable + }; + // -------------------------------------------------------------------------------- + inline bool TTriangle::IsForward( gp_UV uv ) const + { + gp_Vec2d v1( uv, _n1->_uv ), v2( uv, _n2->_uv ); + double d = v1 ^ v2; + return d > 1e-100; + } +} - // not reduced left elements - for (j = 1; j <= free_left; j++) { - // f (i + 1, j + 1) - const SMDS_MeshNode* Nf; - double u,v; - next_base.SetValue(++next_base_len, curr_base.Value(j + 1)); - if (i + 1 == nr) { // top - Nf = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1); - Nf = Nf1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i), - NodesBRD.Value(curr_base.Value(j + 1), i), - NodesBRD.Value(next_base.Value(next_base_len), i + 1), - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); - } +//================================================================================ +/*! + * \brief Set UV of nodes on degenerated VERTEXes in the middle of degenerated EDGE + * + * WARNING: this method must be called AFTER retrieving UVPtStruct's from quad + */ +//================================================================================ - for (int icol = 1; icol <= nb_col; icol++) { +void StdMeshers_Quadrangle_2D::UpdateDegenUV(FaceQuadStruct* quad) +{ + for ( unsigned i = 0; i < quad->side.size(); ++i ) + { + StdMeshers_FaceSide* side = quad->side[i]; + const vector& uvVec = side->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 >= 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" )); - if (remainder > 0 && icol == nb_col && i > remainder / 2) - // stop short "column" - break; + // 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 ); + } +} - // add one "H": nodes b,c,e and faces 1,2,4,5 - // - // .---------b i + 1 - // |\ 5 /| - // | \ / | - // | c---e | - // |1 |2 |4 | - // | | | | - // .--.---.--. i - // - // j j+1 j+2 j+3 - - double u,v; - - // b (i + 1, j + 3) - const SMDS_MeshNode* Nb; - next_base_len++; - next_base.SetValue(next_base_len, curr_base.Value(j + 3)); - if (i + 1 == nr) { // top - Nb = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 3 == curr_base_len) { // right - Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1); - Nb = Nb1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - - // c and d - double u1 = (curr_par_u.Value(j) + next_par_u.Value(next_base_len - 1)) / 2.0; - double u2 = (curr_par_u.Value(j + 3) + next_par_u.Value(next_base_len)) / 2.0; - double u3 = (u2 - u1) / 3.0; - - double v1 = (curr_par_v.Value(j) + next_par_v.Value(next_base_len - 1)) / 2.0; - double v2 = (curr_par_v.Value(j + 3) + next_par_v.Value(next_base_len)) / 2.0; - double v3 = (v2 - v1) / 3.0; - - // c - u = u1 + u3; - v = v1 + v3; - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nc, geomFaceID, u, v); - - // e - u = u1 + u3 + u3; - v = v1 + v3 + v3; - P = S->Value(u,v); - SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Ne, geomFaceID, u, v); - - // Faces - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i), - NodesBRD.Value(curr_base.Value(j + 1), i), - Nc, - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); - - SMDS_MeshFace* F2 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i), - NodesBRD.Value(curr_base.Value(j + 2), i), - Ne, Nc); - if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID); - - SMDS_MeshFace* F4 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i), - NodesBRD.Value(curr_base.Value(j + 3), i), - Nb, Ne); - if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID); - - SMDS_MeshFace* F5 = myTool->AddFace(Nc, Ne, Nb, - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID); - - j += 3; +//================================================================================ +/*! + * \brief Perform smoothing of 2D elements on a FACE with ignored degenerated EDGE + */ +//================================================================================ - // not reduced middle elements - if (icol < nb_col) { - if (remainder > 0 && icol == nb_col - 1 && i > remainder / 2) - // pass middle elements before stopped short "column" - break; +void StdMeshers_Quadrangle_2D::Smooth (FaceQuadStruct* quad) +{ + if ( !myNeedSmooth ) return; - int free_add = free_middle; - if (remainder > 0 && icol == nb_col - 1) - // next "column" is short - free_add -= (nr - 1) - (remainder / 2); + // Get nodes to smooth - for (int imiddle = 1; imiddle <= free_add; imiddle++) { - // f (i + 1, j + imiddle) - const SMDS_MeshNode* Nf; - double u,v; - next_base.SetValue(++next_base_len, curr_base.Value(j + imiddle)); - if (i + 1 == nr) { // top - Nf = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + imiddle == curr_base_len) { // right - Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1); - Nf = Nf1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j - 1 + imiddle), i), - NodesBRD.Value(curr_base.Value(j + imiddle), i), - NodesBRD.Value(next_base.Value(next_base_len), i + 1), - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); - } - j += free_add; - } - } + typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap; + TNo2SmooNoMap smooNoMap; - // not reduced right elements - for (; j < curr_base_len; j++) { - // f (i + 1, j + 1) - const SMDS_MeshNode* Nf; - double u,v; - next_base.SetValue(++next_base_len, curr_base.Value(j + 1)); - if (i + 1 == nr) { // top - Nf = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 1 == curr_base_len) { // right - Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1); - Nf = Nf1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i), - NodesBRD.Value(curr_base.Value(j + 1), i), - NodesBRD.Value(next_base.Value(next_base_len), i + 1), - NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); - if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID); - } - - curr_base_len = next_base_len; - curr_base = next_base; - curr_par_u = next_par_u; - curr_par_v = next_par_v; - next_base_len = 0; - } - } // end "linear" simple reduce "31" - else { + 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 SMDS_MeshNode* node = nIt->next(); + TSmoothNode & sNode = smooNoMap[ node ]; + sNode._uv = myHelper->GetNodeUV( geomFace, node ); + + // set sNode._triangles + SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face ); + while ( fIt->more() ) + { + const SMDS_MeshElement* face = fIt->next(); + const int nbN = face->NbCornerNodes(); + const int nInd = face->GetNodeIndex( node ); + const int prevInd = myHelper->WrapIndex( nInd - 1, nbN ); + const int nextInd = myHelper->WrapIndex( nInd + 1, nbN ); + const SMDS_MeshNode* prevNode = face->GetNode( prevInd ); + const SMDS_MeshNode* nextNode = face->GetNode( nextInd ); + sNode._triangles.push_back( TTriangle( & smooNoMap[ prevNode ], + & smooNoMap[ nextNode ])); } - } // end Simple Reduce implementation + } + // set _uv of smooth nodes on FACE boundary + for ( unsigned i = 0; i < quad->side.size(); ++i ) + { + const vector& uvVec = quad->side[i]->GetUVPtStruct(); + for ( unsigned j = 0; j < uvVec.size(); ++j ) + { + TSmoothNode & sNode = smooNoMap[ uvVec[j].node ]; + sNode._uv.SetCoord( uvVec[j].u, uvVec[j].v ); + } + } - bool isOk = true; - return isOk; + // 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 )); + + // Smoothing + + for ( int iLoop = 0; iLoop < 5; ++iLoop ) + { + for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn ) + { + 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(); + + // 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 ); + + if ( isValid ) + sNode._uv = newUV; + } + } + + // 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 + + 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() ); + + // store the new UV + node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( sNode._uv.X(), sNode._uv.Y() ))); + } + + // Move medium nodes in quadratic mesh + if ( _quadraticMesh ) + { + const TLinkNodeMap& links = myHelper->GetTLinkNodeMap(); + TLinkNodeMap::const_iterator linkIt = links.begin(); + for ( ; linkIt != links.end(); ++linkIt ) + { + const SMESH_TLink& link = linkIt->first; + SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( linkIt->second ); + + 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_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() ); + } + } }