X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=3842b4733f56d8d937d25444bfb7de67257cb272;hb=2c8f4c513e0fed76a7f0331542e520e315740aa4;hp=174a2441128de1ef6d992212dd4f9a75eebd0246;hpb=8fa039a796957b302d86d90b22afc0a998573f83;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index 174a24411..3842b4733 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 @@ -86,7 +87,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; + myHelper = 0; } //============================================================================= @@ -202,9 +203,10 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, aMesh.GetSubMesh(aShape); SMESH_MesherHelper helper (aMesh); - myTool = &helper; + myHelper = &helper; - _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 +224,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 +241,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; } } @@ -305,7 +311,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, b = quad->uv_grid[j * nbhoriz + i + 1].node; c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node; d = quad->uv_grid[(j + 1) * nbhoriz + i].node; - SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + 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; } @@ -792,6 +801,7 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes FaceQuadStruct* 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(); @@ -843,13 +853,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, nbSides 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() @@ -867,10 +880,30 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes sideEdges.splice(sideEdges.begin(), edges, --edges.end()); } } + if ( sideEdges.size() == 1 && BRep_Tool::Degenerated( sideEdges.front() )) + degenSides.push_back( nbSides ); + quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, nbSidesside.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 @@ -913,7 +946,7 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes for (int i = 0; i < nbSides; ++i) { MESSAGE (" ("); for (int e = 0; e < quad->side[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; @@ -1229,6 +1262,9 @@ bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, //return error("Can't find nodes on sides"); return error(COMPERR_BAD_INPUT_MESH); + if ( myNeedSmooth ) + UpdateDegenUV( quad ); + // nodes Id on "in" edges if (! quad->isEdgeOut[0]) { int j = 0; @@ -1489,6 +1525,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); } @@ -1658,13 +1697,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); } @@ -1743,13 +1782,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); } @@ -1781,13 +1820,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); } @@ -1893,13 +1932,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); } @@ -1926,13 +1965,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); } @@ -2087,11 +2126,11 @@ bool StdMeshers_Quadrangle_2D::EvaluateQuadPref(SMESH_Mesh & aMesh, */ //============================================================================= 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) + 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()); @@ -2099,20 +2138,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) @@ -2132,7 +2348,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) // // .-----.-----.-----.-----. .-----.-----.-----.-----. // | / \ | / \ | | / \ | / \ | @@ -2153,7 +2369,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) { @@ -2162,7 +2377,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; @@ -2178,25 +2392,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; - int npair_top = ncol_top / 2; - // maximum number of bottom elements for "linear" simple reduce - //int max_lin = ncol_top + npair_top * 2 * nrows; - // maximum number of bottom elements for "tree" simple reduce - int max_tree = npair_top * pow(2.0, nrows + 1); - if (ncol_top > npair_top * 2) { - int delta = ncol_bot - max_tree; - for (int irow = 1; irow < nrows; irow++) { - int nfour = delta / 4; - delta -= nfour*2; - } - if (delta <= (ncol_top - npair_top * 2)) - max_tree = ncol_bot; - } - - if (ncol_bot > max_tree) + // 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; } @@ -2232,7 +2433,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; @@ -2254,6 +2455,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++) { @@ -2303,7 +2507,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; @@ -2352,8 +2556,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); } } @@ -2407,8 +2611,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); } } @@ -2471,7 +2675,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); } @@ -2505,19 +2709,66 @@ 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 - int max_lin = ncol_top + npair_top * 2 * nrows; - // maximum number of bottom elements for "tree" simple reduce - //int max_tree = npair_top * pow(2, nrows + 1); - - //if (ncol_bot > max_tree) - // MultipleReduce = true; + // 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 = 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; + } + } + // maximum number of bottom elements for "tree" simple reduce 3->1 + //int max_tree31 = ncol_top * pow(3.0, nrows); + bool is_lin_31 = false; + bool is_lin_42 = false; + bool is_tree_31 = false; + bool is_tree_42 = false; + int max_lin = max_lin42; + if (ncol_bot > max_lin42) { + if (ncol_bot <= max_lin31) { + is_lin_31 = true; + max_lin = max_lin31; + } + } + else { + // if ncol_bot is a 3*n or not 2*n + if ((ncol_bot/3)*3 == ncol_bot || (ncol_bot/2)*2 != ncol_bot) { + is_lin_31 = true; + max_lin = max_lin31; + } + else { + is_lin_42 = true; + } + } + if (ncol_bot > max_lin) { // not "linear" + is_tree_31 = (ncol_bot > max_tree42); + if (ncol_bot <= max_tree42) { + if ((ncol_bot/3)*3 == ncol_bot || (ncol_bot/2)*2 != ncol_bot) { + is_tree_31 = true; + } + else { + is_tree_42 = true; + } + } + } const vector& uv_eb = quad->side[0]->GetUVPtStruct(true,0); const vector& uv_er = quad->side[1]->GetUVPtStruct(false,1); @@ -2527,62 +2778,197 @@ 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 (ncol_bot > max_lin) { - // "tree" simple reduce 2->4->8->16->32->... + 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 + // .-------------------------------.-------------------------------. nr + // | \ | / | + // | \ .---------------.---------------. / | + // | | | | | + // .---------------.---------------.---------------.---------------. // | \ | / | \ | / | // | \ .-------.-------. / | \ .-------.-------. / | // | | | | | | | | | - // .-------.-------.-------.-------.-------.-------.-------.-------. + // .-------.-------.-------.-------.-------.-------.-------.-------. i // |\ | /|\ | /|\ | /|\ | /| - // | \.---.---./ | \.---.---./ | \.---.---./ | \.---.---./ | i + // | \.---.---./ | \.---.---./ | \.---.---./ | \.---.---./ | // | | | | | | | | | | | | | | | | | // .---.---.---.---.---.---.---.---.---.---.---.---.---.---.---.---. // |\ | /|\ | /|\ | /|\ | /|\ | /|\ | /|\ | /|\ | /| @@ -2591,237 +2977,102 @@ 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); - + // "tree" simple reduce "31": 1->3->9->27->... + // + // .-----------------------------------------------------. nr + // | \ / | + // | .-----------------. | + // | | | | + // .-----------------.-----------------.-----------------. + // | \ / | \ / | \ / | + // | .-----. | .-----. | .-----. | i + // | | | | | | | | | | + // .-----.-----.-----.-----.-----.-----.-----.-----.-----. + // |\ /|\ /|\ /|\ /|\ /|\ /|\ /|\ /|\ /| + // | .-. | .-. | .-. | .-. | .-. | .-. | .-. | .-. | .-. | + // | | | | | | | | | | | | | | | | | | | | | | | | | | | | + // .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1 + // 1 j nb + + PReduceFunction reduceFunction = & ( is_tree_42 ? reduce42 : reduce31 ); + + 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; - //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); + 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 + 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); + 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 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; + 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 - else { - // "linear" simple reduce 4->8->12->16 (3 steps) + + 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 // | \ | / | \ | / | @@ -2856,7 +3107,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) { @@ -2873,17 +3125,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) @@ -2891,46 +3138,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++) { @@ -2939,144 +3180,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) { @@ -3091,109 +3198,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 + + else { + return false; + } } // end Simple Reduce implementation bool isOk = true; return isOk; } + +//================================================================================ +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; + } +} + +//================================================================================ +/*! + * \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 + */ +//================================================================================ + +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" )); + + // 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 ); + } +} + +//================================================================================ +/*! + * \brief Perform smoothing of 2D elements on a FACE with ignored degenerated EDGE + */ +//================================================================================ + +void StdMeshers_Quadrangle_2D::Smooth (FaceQuadStruct* 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 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 ])); + } + } + // 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 ); + } + } + + // 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() ); + } + } +}