X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=23be5d62faaab0306318ae709e11d59422abe35a;hb=f4f14a98e3d222a333b9f1b8cbb1a2a29824bc2b;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..23be5d62f 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-2011 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,16 +2138,16 @@ 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); } } @@ -2153,7 +2192,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 +2200,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 +2215,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; } @@ -2254,6 +2278,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++) { @@ -2352,7 +2379,7 @@ 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), + 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,7 +2434,7 @@ 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), + 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 +2498,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); } @@ -2511,13 +2538,63 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, 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 + // 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 "tree" simple reduce - //int max_tree = npair_top * pow(2, nrows + 1); + // maximum number of bottom elements for "linear" simple reduce 4->2 + 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" +#ifdef WIN32 + // of the MSVC doesn't contain log2 + int nrows_tree42 = int( log( (double)(ncol_bot / ncol_top) )/log((double)2) ); // needed to avoid overflow at pow(2) +#else + int nrows_tree42 = int( log2( ncol_bot / ncol_top )); // needed to avoid overflow at pow(2) +#endif - //if (ncol_bot > max_tree) - // MultipleReduce = true; + if (ncol_top > npair_top * 2 && nrows_tree42 < nrows) { + max_tree42 = npair_top * pow(2.0, nrows + 1); + int delta = ncol_bot - int( 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; + if (ncol_bot > max_lin) { + 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); @@ -2573,16 +2650,20 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, int curr_base_len = nb; int next_base_len = 0; - if (ncol_bot > max_lin) { - // "tree" simple reduce 2->4->8->16->32->... + if (is_tree_42) { + // "tree" simple reduce "42": 2->4->8->16->32->... // - // .---------------.---------------.---------------.---------------. nr + // .-------------------------------.-------------------------------. nr + // | \ | / | + // | \ .---------------.---------------. / | + // | | | | | + // .---------------.---------------.---------------.---------------. // | \ | / | \ | / | // | \ .-------.-------. / | \ .-------.-------. / | // | | | | | | | | | - // .-------.-------.-------.-------.-------.-------.-------.-------. + // .-------.-------.-------.-------.-------.-------.-------.-------. i // |\ | /|\ | /|\ | /|\ | /| - // | \.---.---./ | \.---.---./ | \.---.---./ | \.---.---./ | i + // | \.---.---./ | \.---.---./ | \.---.---./ | \.---.---./ | // | | | | | | | | | | | | | | | | | // .---.---.---.---.---.---.---.---.---.---.---.---.---.---.---.---. // |\ | /|\ | /|\ | /|\ | /|\ | /|\ | /|\ | /|\ | /| @@ -2732,32 +2813,32 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, meshDS->SetNodeOnFace(Ne, geomFaceID, u, v); // Faces - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i), + SMDS_MeshFace* F1 = myHelper->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), + SMDS_MeshFace* F2 = myHelper->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), + SMDS_MeshFace* F3 = myHelper->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), + SMDS_MeshFace* F4 = myHelper->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, + SMDS_MeshFace* F5 = myHelper->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); + SMDS_MeshFace* F6 = myHelper->AddFace(Nd, Ne, Nb, Na); if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID); } @@ -2806,7 +2887,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, } 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), + SMDS_MeshFace* F1 = myHelper->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)); @@ -2819,9 +2900,198 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, curr_par_v = next_par_v; next_base_len = 0; } - } // end "tree" simple reduce - else { - // "linear" simple reduce 4->8->12->16 (3 steps) + } // end "tree" simple reduce "42" + else if (is_tree_31) { + // "tree" simple reduce "31": 1->3->9->27->... + // + // .-----------------------------------------------------. nr + // | \ / | + // | .-----------------. | + // | | | | + // .-----------------.-----------------.-----------------. + // | \ / | \ / | \ / | + // | .-----. | .-----. | .-----. | i + // | | | | | | | | | | + // .-----.-----.-----.-----.-----.-----.-----.-----.-----. + // |\ /|\ /|\ /|\ /|\ /|\ /|\ /|\ /|\ /| + // | .-. | .-. | .-. | .-. | .-. | .-. | .-. | .-. | .-. | + // | | | | | | | | | | | | | | | | | | | | | | | | | | | | + // .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 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; + + // 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); + 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 = myHelper->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 = myHelper->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 = myHelper->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 = myHelper->AddFace(Nc, Ne, Nb, + NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); + if (F5) meshDS->SetMeshElementOnShape(F5, 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 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 = myHelper->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 "31" + else if (is_lin_42) { + // "linear" simple reduce "42": 4->8->12->16 // // .---------------.---------------.---------------.---------------. nr // | \ | / | \ | / | @@ -2926,7 +3196,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, } 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), + SMDS_MeshFace* F1 = myHelper->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)); @@ -3048,32 +3318,32 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, meshDS->SetNodeOnFace(Ne, geomFaceID, u, v); // Faces - SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i), + SMDS_MeshFace* F1 = myHelper->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), + SMDS_MeshFace* F2 = myHelper->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), + SMDS_MeshFace* F3 = myHelper->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), + SMDS_MeshFace* F4 = myHelper->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, + SMDS_MeshFace* F5 = myHelper->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); + SMDS_MeshFace* F6 = myHelper->AddFace(Nd, Ne, Nb, Na); if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID); j += 4; @@ -3128,7 +3398,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, } 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), + SMDS_MeshFace* F1 = myHelper->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)); @@ -3178,22 +3448,549 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, } 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), + SMDS_MeshFace* F1 = myHelper->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 "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; + } + + 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); + + // left + next_par_u.SetValue(next_base_len, uv_el[i].u); + next_par_v.SetValue(next_base_len, uv_el[i].v); + + // 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; + + // 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 = myHelper->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); } + for (int icol = 1; icol <= nb_col; icol++) { + + if (remainder > 0 && icol == nb_col && i > remainder / 2) + // stop short "column" + break; + + // 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 = myHelper->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 = myHelper->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 = myHelper->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 = myHelper->AddFace(Nc, Ne, Nb, + NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1)); + if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID); + + j += 3; + + // 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; + + int free_add = free_middle; + if (remainder > 0 && icol == nb_col - 1) + // next "column" is short + free_add -= (nr - 1) - (remainder / 2); + + 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 = myHelper->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; + } + } + + // 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 = myHelper->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 + } // end "linear" simple reduce "31" + else { + } } // 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() ); + } + } +}