X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=d5c4c9ad788c34b02eabf65666a4f478f16e3493;hb=3fa13a8549e50ceb8cb010ff22435aaa0c5a1573;hp=3589a7870609db0a94ccff51a622ac9595b32ff5;hpb=91c92cb54310225231438b4d3bafeb0d1643a7c0;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index 3589a7870..d5c4c9ad7 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2011 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 @@ -26,21 +26,19 @@ #include "StdMeshers_Quadrangle_2D.hxx" -#include "StdMeshers_FaceSide.hxx" - -#include "StdMeshers_QuadrangleParams.hxx" - +#include "SMDS_EdgePosition.hxx" +#include "SMDS_FacePosition.hxx" +#include "SMDS_MeshElement.hxx" +#include "SMDS_MeshNode.hxx" +#include "SMESH_Block.hxx" +#include "SMESH_Comment.hxx" #include "SMESH_Gen.hxx" #include "SMESH_Mesh.hxx" -#include "SMESH_subMesh.hxx" #include "SMESH_MesherHelper.hxx" -#include "SMESH_Block.hxx" -#include "SMESH_Comment.hxx" - -#include "SMDS_MeshElement.hxx" -#include "SMDS_MeshNode.hxx" -#include "SMDS_EdgePosition.hxx" -#include "SMDS_FacePosition.hxx" +#include "SMESH_subMesh.hxx" +#include "StdMeshers_FaceSide.hxx" +#include "StdMeshers_QuadrangleParams.hxx" +#include "StdMeshers_ViscousLayers2D.hxx" #include #include @@ -79,7 +77,8 @@ typedef SMESH_Comment TComm; StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMESH_Gen* gen) - : SMESH_2D_Algo(hypId, studyId, gen) + : SMESH_2D_Algo(hypId, studyId, gen), + myHelper( 0 ) { MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D"); _name = "Quadrangle_2D"; @@ -87,7 +86,7 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, _compatibleHypothesis.push_back("QuadrangleParams"); _compatibleHypothesis.push_back("QuadranglePreference"); _compatibleHypothesis.push_back("TrianglePreference"); - myHelper = 0; + _compatibleHypothesis.push_back("ViscousLayers2D"); } //============================================================================= @@ -193,11 +192,11 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis */ //============================================================================= -bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape)// throw (SALOME_Exception) +bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape) { - // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside - //Unexpect aCatchSalomeException); + const TopoDS_Face& F = TopoDS::Face(aShape); + Handle(Geom_Surface) S = BRep_Tool::Surface(F); SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); aMesh.GetSubMesh(aShape); @@ -205,6 +204,10 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, SMESH_MesherHelper helper (aMesh); myHelper = &helper; + myProxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F ); + if ( !myProxyMesh ) + return false; + _quadraticMesh = myHelper->IsQuadraticSubMesh(aShape); myNeedSmooth = false; @@ -263,9 +266,6 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, int nbhoriz = Min(nbdown, nbup); int nbvertic = Min(nbright, nbleft); - const TopoDS_Face& F = TopoDS::Face(aShape); - Handle(Geom_Surface) S = BRep_Tool::Surface(F); - // internal mesh nodes int i, j, geomFaceID = meshDS->ShapeToIndex(F); for (i = 1; i < nbhoriz - 1; i++) { @@ -307,10 +307,10 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, for (i = ilow; i < iup; i++) { for (j = jlow; j < jup; j++) { const SMDS_MeshNode *a, *b, *c, *d; - a = quad->uv_grid[j * nbhoriz + i].node; - b = quad->uv_grid[j * nbhoriz + i + 1].node; + a = quad->uv_grid[j * nbhoriz + i ].node; + b = quad->uv_grid[j * nbhoriz + i + 1].node; c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node; - d = quad->uv_grid[(j + 1) * nbhoriz + i].node; + d = quad->uv_grid[(j + 1) * nbhoriz + i ].node; SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d); if (face) { meshDS->SetMeshElementOnShape(face, geomFaceID); @@ -831,9 +831,16 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes } if (!E1.IsNull() && !E2.IsNull() && !E3.IsNull()) { - quad->side.push_back(new StdMeshers_FaceSide(F, E1, &aMesh, true, ignoreMediumNodes)); - quad->side.push_back(new StdMeshers_FaceSide(F, E2, &aMesh, true, ignoreMediumNodes)); - quad->side.push_back(new StdMeshers_FaceSide(F, E3, &aMesh, false,ignoreMediumNodes)); + if ( myProxyMesh->GetProxySubMesh( E1 ) || + myProxyMesh->GetProxySubMesh( E2 ) || + myProxyMesh->GetProxySubMesh( E3 ) ) + + quad->side.push_back(new StdMeshers_FaceSide(F, E1, &aMesh, true, + ignoreMediumNodes, myProxyMesh)); + quad->side.push_back(new StdMeshers_FaceSide(F, E2, &aMesh, true, + ignoreMediumNodes, myProxyMesh)); + quad->side.push_back(new StdMeshers_FaceSide(F, E3, &aMesh, false, + ignoreMediumNodes, myProxyMesh)); const vector& UVPSleft = quad->side[0]->GetUVPtStruct(true,0); /* vector& UVPStop = */quad->side[1]->GetUVPtStruct(false,1); /* vector& UVPSright = */quad->side[2]->GetUVPtStruct(true,1); @@ -856,8 +863,8 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes else if (nbEdgesInWire.front() == 4) // exactly 4 edges { for (; edgeIt != edges.end(); ++edgeIt, nbSides++) - quad->side.push_back(new StdMeshers_FaceSide(F, *edgeIt, &aMesh, - nbSidesside.push_back(new StdMeshers_FaceSide(F, *edgeIt, &aMesh, nbSides < TOP_SIDE, + ignoreMediumNodes, myProxyMesh)); } else if (nbEdgesInWire.front() > 4) // more than 4 edges - try to unite some { @@ -883,8 +890,8 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes if ( sideEdges.size() == 1 && BRep_Tool::Degenerated( sideEdges.front() )) degenSides.push_back( nbSides ); - quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, - nbSidesside.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, nbSides < TOP_SIDE, + ignoreMediumNodes, myProxyMesh)); ++nbSides; } if ( !degenSides.empty() && nbSides - degenSides.size() == 4 ) @@ -895,9 +902,9 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes for ( int i = degenSides.size()-1; i > -1; --i ) { - StdMeshers_FaceSide * & degenSide = quad->side[ degenSides[ i ]]; + StdMeshers_FaceSide* degenSide = quad->side[ degenSides[ i ]]; delete degenSide; - quad->side.erase( vector::iterator( & degenSide )); + quad->side.erase( quad->side.begin() + degenSides[ i ] ); } for ( unsigned i = TOP_SIDE; i < quad->side.size(); ++i ) quad->side[i]->Reverse(); @@ -935,7 +942,8 @@ FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMes } } quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, - nbSidesisEdgeOut[0]) { - int j = 0; - for (int i = 0; i < nbhoriz; i++) { // down - int ij = j * nbhoriz + i; - uv_grid[ij].node = uv_e0[i].node; - } + // copy data of face boundary + /*if (! quad->isEdgeOut[0])*/ { + const int j = 0; + for (int i = 0; i < nbhoriz; i++) // down + uv_grid[ j * nbhoriz + i ] = uv_e0[i]; } - if (! quad->isEdgeOut[1]) { - int i = nbhoriz - 1; - for (int j = 0; j < nbvertic; j++) { // right - int ij = j * nbhoriz + i; - uv_grid[ij].node = uv_e1[j].node; - } + /*if (! quad->isEdgeOut[1])*/ { + const int i = nbhoriz - 1; + for (int j = 0; j < nbvertic; j++) // right + uv_grid[ j * nbhoriz + i ] = uv_e1[j]; } - if (! quad->isEdgeOut[2]) { - int j = nbvertic - 1; - for (int i = 0; i < nbhoriz; i++) { // up - int ij = j * nbhoriz + i; - uv_grid[ij].node = uv_e2[i].node; - } + /*if (! quad->isEdgeOut[2])*/ { + const int j = nbvertic - 1; + for (int i = 0; i < nbhoriz; i++) // up + uv_grid[ j * nbhoriz + i ] = uv_e2[i]; } - if (! quad->isEdgeOut[3]) { + /*if (! quad->isEdgeOut[3])*/ { int i = 0; - for (int j = 0; j < nbvertic; j++) { // left - int ij = j * nbhoriz + i; - uv_grid[ij].node = uv_e3[j].node; - } + for (int j = 0; j < nbvertic; j++) // left + uv_grid[ j * nbhoriz + i ] = uv_e3[j]; } - // normalized 2d values on grid + // normalized 2d parameters on grid for (int i = 0; i < nbhoriz; i++) { for (int j = 0; j < nbvertic; j++) { int ij = j * nbhoriz + i; @@ -1316,26 +1316,23 @@ bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, } // 4 --- projection on 2d domain (u,v) - gp_UV a0(uv_e0.front().u, uv_e0.front().v); - gp_UV a1(uv_e0.back().u, uv_e0.back().v); - gp_UV a2(uv_e2.back().u, uv_e2.back().v); - gp_UV a3(uv_e2.front().u, uv_e2.front().v); + gp_UV a0 (uv_e0.front().u, uv_e0.front().v); + gp_UV a1 (uv_e0.back().u, uv_e0.back().v ); + gp_UV a2 (uv_e2.back().u, uv_e2.back().v ); + gp_UV a3 (uv_e2.front().u, uv_e2.front().v); + + for (int i = 0; i < nbhoriz; i++) + { + gp_UV p0( uv_e0[i].u, uv_e0[i].v ); + gp_UV p2( uv_e2[i].u, uv_e2[i].v ); + for (int j = 0; j < nbvertic; j++) + { + gp_UV p1( uv_e1[j].u, uv_e1[j].v ); + gp_UV p3( uv_e3[j].u, uv_e3[j].v ); - for (int i = 0; i < nbhoriz; i++) { - for (int j = 0; j < nbvertic; j++) { int ij = j * nbhoriz + i; double x = uv_grid[ij].x; double y = uv_grid[ij].y; - double param_0 = uv_e0[0].normParam + x * (uv_e0.back().normParam - uv_e0[0].normParam); // sud - double param_2 = uv_e2[0].normParam + x * (uv_e2.back().normParam - uv_e2[0].normParam); // nord - double param_1 = uv_e1[0].normParam + y * (uv_e1.back().normParam - uv_e1[0].normParam); // est - double param_3 = uv_e3[0].normParam + y * (uv_e3.back().normParam - uv_e3[0].normParam); // ouest - - //MESSAGE("params "<side[0]->Value2d(param_0).XY(); - gp_UV p1 = quad->side[1]->Value2d(param_1).XY(); - gp_UV p2 = quad->side[2]->Value2d(param_2).XY(); - gp_UV p3 = quad->side[3]->Value2d(param_3).XY(); gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3); @@ -1353,7 +1350,8 @@ bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, static void ShiftQuad(FaceQuadStruct* quad, const int num, bool) { - StdMeshers_FaceSide* side[4] = { quad->side[0], quad->side[1], quad->side[2], quad->side[3] }; + StdMeshers_FaceSide* side[4] = { quad->side[0], quad->side[1], + quad->side[2], quad->side[3] }; for (int i = BOTTOM_SIDE; i < NB_SIDES; ++i) { int id = (i + num) % NB_SIDES; bool wasForward = (i < TOP_SIDE); @@ -1374,23 +1372,18 @@ static gp_UV CalcUV(double x0, double x1, double y0, double y1, const gp_UV& a0, const gp_UV& a1, const gp_UV& a2, const gp_UV& a3) { - const vector& uv_eb = quad->side[0]->GetUVPtStruct(true,0); - const vector& uv_er = quad->side[1]->GetUVPtStruct(false,1); - const vector& uv_et = quad->side[2]->GetUVPtStruct(true,1); - const vector& uv_el = quad->side[3]->GetUVPtStruct(false,0); + // const vector& uv_eb = quad->side[0]->GetUVPtStruct(true,0); + // const vector& uv_er = quad->side[1]->GetUVPtStruct(false,1); + // const vector& uv_et = quad->side[2]->GetUVPtStruct(true,1); + // const vector& uv_el = quad->side[3]->GetUVPtStruct(false,0); double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); double y = y0 + x * (y1 - y0); - double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam); - double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam); - double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam); - double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam); - - gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY(); - gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY(); - gp_UV p2 = quad->side[TOP_SIDE ]->Value2d(param_t).XY(); - gp_UV p3 = quad->side[LEFT_SIDE ]->Value2d(param_l).XY(); + gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(x).XY(); + gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(y).XY(); + gp_UV p2 = quad->side[TOP_SIDE ]->Value2d(x).XY(); + gp_UV p3 = quad->side[LEFT_SIDE ]->Value2d(y).XY(); gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3); @@ -2152,6 +2145,183 @@ void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS, } } +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) @@ -2171,7 +2341,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, int nt = quad->side[2]->NbPoints(); int nl = quad->side[3]->NbPoints(); - // Simple Reduce 8->6->4->2 (3 steps) Multiple Reduce 8->2 (1 step) + // Simple Reduce 10->8->6->4 (3 steps) Multiple Reduce 10->4 (1 step) // // .-----.-----.-----.-----. .-----.-----.-----.-----. // | / \ | / \ | | / \ | / \ | @@ -2219,7 +2389,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, int ncol_top = nt1 - 1; int ncol_bot = nb1 - 1; // number of rows needed to reduce ncol_bot to ncol_top using simple 3->1 "tree" (see below) - int nrows_tree31 = int( log( ncol_bot / ncol_top ) / log( 3 )); // = log x base 3 + int nrows_tree31 = int( log( (double)(ncol_bot / ncol_top) ) / log((double) 3 )); // = log x base 3 if ( nrows < nrows_tree31 ) MultipleReduce = true; } @@ -2256,7 +2426,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, nl = quad->side[3]->NbPoints(); dh = abs(nb-nt); dv = abs(nr-nl); - int nbh = Max(nb,nt); + int nbh = Max(nb,nt); int nbv = Max(nr,nl); int addh = 0; int addv = 0; @@ -2330,7 +2500,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, gp_XY a3 (uv_et.front().u, uv_et.front().v); int nnn = Min(nr,nl); - // auxilary sequence of XY for creation nodes + // auxilary sequence of XY for creation of nodes // in the bottom part of central domain // it's length must be == nbv-nnn-1 TColgp_SequenceOfXY UVL; @@ -2380,7 +2550,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, for (j=1; jAddFace(NodesL.Value(i,j), NodesL.Value(i+1,j), - NodesL.Value(i+1,j+1), NodesL.Value(i,j+1)); + NodesL.Value(i+1,j+1), NodesL.Value(i,j+1)); if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } } @@ -2435,7 +2605,7 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, for (j=1; jAddFace(NodesR.Value(i,j), NodesR.Value(i+1,j), - NodesR.Value(i+1,j+1), NodesR.Value(i,j+1)); + NodesR.Value(i+1,j+1), NodesR.Value(i,j+1)); if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } } @@ -2532,29 +2702,31 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, nr = quad->side[1]->NbPoints(); nt = quad->side[2]->NbPoints(); nl = quad->side[3]->NbPoints(); - + // number of rows and columns int nrows = nr - 1; // and also == nl - 1 int ncol_top = nt - 1; int ncol_bot = nb - 1; int npair_top = ncol_top / 2; // maximum number of bottom elements for "linear" simple reduce 4->2 - int max_lin = ncol_top + npair_top * 2 * nrows; - // maximum number of bottom elements for "linear" simple reduce 4->2 + int max_lin42 = ncol_top + npair_top * 2 * nrows; + // maximum number of bottom elements for "linear" simple reduce 3->1 int max_lin31 = ncol_top + ncol_top * 2 * nrows; // maximum number of bottom elements for "tree" simple reduce 4->2 int max_tree42 = 0; // number of rows needed to reduce ncol_bot to ncol_top using simple 4->2 "tree" - int nrows_tree42 = int( log2( ncol_bot / ncol_top )); // needed to avoid overflow at pow(2) - if (ncol_top > npair_top * 2 && nrows_tree42 < nrows) { + 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); - int delta = ncol_bot - int( max_tree42 ); - for (int irow = 1; irow < nrows; irow++) { - int nfour = delta / 4; - delta -= nfour * 2; + if ( ncol_top > npair_top * 2 ) { + int delta = ncol_bot - max_tree42; + for (int irow = 1; irow < nrows; irow++) { + int nfour = delta / 4; + delta -= nfour * 2; + } + if (delta <= (ncol_top - npair_top * 2)) + max_tree42 = ncol_bot; } - if (delta <= (ncol_top - npair_top * 2)) - max_tree42 = ncol_bot; } // maximum number of bottom elements for "tree" simple reduce 3->1 //int max_tree31 = ncol_top * pow(3.0, nrows); @@ -2562,7 +2734,8 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, bool is_lin_42 = false; bool is_tree_31 = false; bool is_tree_42 = false; - if (ncol_bot > max_lin) { + int max_lin = max_lin42; + if (ncol_bot > max_lin42) { if (ncol_bot <= max_lin31) { is_lin_31 = true; max_lin = max_lin31; @@ -2598,53 +2771,184 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl) return error(COMPERR_BAD_INPUT_MESH); - // arrays for normalized params - TColStd_SequenceOfReal npb, npr, npt, npl; - for (j = 0; j < nb; j++) { - npb.Append(uv_eb[j].normParam); - } - for (i = 0; i < nr; i++) { - npr.Append(uv_er[i].normParam); - } - for (j = 0; j < nt; j++) { - npt.Append(uv_et[j].normParam); - } - for (i = 0; i < nl; i++) { - npl.Append(uv_el[i].normParam); - } + myHelper->SetElementsOnShape( true ); - // We will ajust new points to this grid - if (!SetNormalizedGrid(aMesh, aShape, quad)) - return false; + gp_UV uv[ UV_SIZE ]; + uv[ UV_A0 ].SetCoord( uv_eb.front().u, uv_eb.front().v); + uv[ UV_A1 ].SetCoord( uv_eb.back().u, uv_eb.back().v ); + uv[ UV_A2 ].SetCoord( uv_et.back().u, uv_et.back().v ); + uv[ UV_A3 ].SetCoord( uv_et.front().u, uv_et.front().v); - // TODO ??? - gp_XY a0 (uv_eb.front().u, uv_eb.front().v); - gp_XY a1 (uv_eb.back().u, uv_eb.back().v); - gp_XY a2 (uv_et.back().u, uv_et.back().v); - gp_XY a3 (uv_et.front().u, uv_et.front().v); - //========================================================= + vector curr_base = uv_eb, next_base; - TColStd_SequenceOfInteger curr_base, next_base; - TColStd_SequenceOfReal curr_par_u, curr_par_v; - TColStd_SequenceOfReal next_par_u, next_par_v; - StdMeshers_Array2OfNode NodesBRD (1,nb, 1,nr); - for (j = 1; j <= nb; j++) { - NodesBRD.SetValue(j, 1, uv_eb[j - 1].node); // bottom - curr_base.Append(j); - next_base.Append(-1); - curr_par_u.Append(uv_eb[j-1].u); - curr_par_v.Append(uv_eb[j-1].v); - next_par_u.Append(0.); - next_par_v.Append(0.); - } - for (j = 1; j <= nt; j++) { - NodesBRD.SetValue(j, nr, uv_et[j - 1].node); // top - } + UVPtStruct nullUVPtStruct; nullUVPtStruct.node = 0; int curr_base_len = nb; int next_base_len = 0; - if (is_tree_42) { + if ( true ) + { // ------------------------------------------------------------------ + // New algorithm implemented by request of IPAL22856 + // "2D quadrangle mesher of reduced type works wrong" + // http://bugtracker.opencascade.com/show_bug.cgi?id=22856 + + // the algorithm is following: all reduces are centred in horizontal + // direction and are distributed among all rows + + if (ncol_bot > max_tree42) { + is_lin_31 = true; + } + else { + if ((ncol_top/3)*3 == ncol_top ) { + is_lin_31 = true; + } + else { + is_lin_42 = true; + } + } + + const int col_top_size = is_lin_42 ? 2 : 1; + const int col_base_size = is_lin_42 ? 4 : 3; + + // Compute nb of "columns" (like in "linear" simple reducing) in all rows + + vector nb_col_by_row; + + int delta_all = nb - nt; + int delta_one_col = nrows * 2; + int nb_col = delta_all / delta_one_col; + int remainder = delta_all - nb_col * delta_one_col; + if (remainder > 0) { + nb_col++; + } + if ( nb_col * col_top_size >= nt ) // == "tree" reducing situation + { + // top row is full (all elements reduced), add "columns" one by one + // in rows below until all bottom elements are reduced + nb_col = ( nt - 1 ) / col_top_size; + nb_col_by_row.resize( nrows, nb_col ); + int nbrows_not_full = nrows - 1; + int cur_top_size = nt - 1; + remainder = delta_all - nb_col * delta_one_col; + while ( remainder > 0 ) + { + delta_one_col = nbrows_not_full * 2; + int nb_col_add = remainder / delta_one_col; + cur_top_size += 2 * nb_col_by_row[ nbrows_not_full ]; + int nb_col_free = cur_top_size / col_top_size - nb_col_by_row[ nbrows_not_full-1 ]; + if ( nb_col_add > nb_col_free ) + nb_col_add = nb_col_free; + for ( int irow = 0; irow < nbrows_not_full; ++irow ) + nb_col_by_row[ irow ] += nb_col_add; + nbrows_not_full --; + remainder -= nb_col_add * delta_one_col; + } + } + else // == "linear" reducing situation + { + nb_col_by_row.resize( nrows, nb_col ); + if (remainder > 0) + for ( int irow = remainder / 2; irow < nrows; ++irow ) + nb_col_by_row[ irow ]--; + } + + // Make elements + + PReduceFunction reduceFunction = & ( is_lin_42 ? reduce42 : reduce31 ); + + const int reduce_grp_size = is_lin_42 ? 4 : 3; + + for (i = 1; i < nr; i++) // layer by layer + { + nb_col = nb_col_by_row[ i-1 ]; + int nb_next = curr_base_len - nb_col * 2; + if (nb_next < nt) nb_next = nt; + + const double y = uv_el[ i ].normParam; + + if ( i + 1 == nr ) // top + { + next_base = uv_et; + } + else + { + next_base.clear(); + next_base.resize( nb_next, nullUVPtStruct ); + next_base.front() = uv_el[i]; + next_base.back() = uv_er[i]; + + // compute normalized param u + double du = 1. / ( nb_next - 1 ); + next_base[0].normParam = 0.; + for ( j = 1; j < nb_next; ++j ) + next_base[j].normParam = next_base[j-1].normParam + du; + } + uv[ UV_L ].SetCoord( next_base.front().u, next_base.front().v ); + uv[ UV_R ].SetCoord( next_base.back().u, next_base.back().v ); + + int free_left = ( curr_base_len - 1 - nb_col * col_base_size ) / 2; + int free_middle = curr_base_len - 1 - nb_col * col_base_size - 2 * free_left; + + // not reduced left elements + for (j = 0; j < free_left; j++) + { + // f (i + 1, j + 1) + const SMDS_MeshNode*& Nf = next_base[++next_base_len].node; + if ( !Nf ) + Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S ); + + myHelper->AddFace(curr_base[ j ].node, + curr_base[ j+1 ].node, + Nf, + next_base[ next_base_len-1 ].node); + } + + for (int icol = 1; icol <= nb_col; icol++) + { + // add "H" + reduceFunction( curr_base, next_base, j, next_base_len, quad, uv, y, myHelper, S ); + + j += reduce_grp_size; + + // elements in the middle of "columns" added for symmetry + if ( free_middle > 0 && ( nb_col % 2 == 0 ) && icol == nb_col / 2 ) + { + for (int imiddle = 1; imiddle <= free_middle; imiddle++) { + // f (i + 1, j + imiddle) + const SMDS_MeshNode*& Nf = next_base[++next_base_len].node; + if ( !Nf ) + Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S ); + + myHelper->AddFace(curr_base[ j-1+imiddle ].node, + curr_base[ j +imiddle ].node, + Nf, + next_base[ next_base_len-1 ].node); + } + j += free_middle; + } + } + + // not reduced right elements + for (; j < curr_base_len-1; j++) { + // f (i + 1, j + 1) + const SMDS_MeshNode*& Nf = next_base[++next_base_len].node; + if ( !Nf ) + Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S ); + + myHelper->AddFace(curr_base[ j ].node, + curr_base[ j+1 ].node, + Nf, + next_base[ next_base_len-1 ].node); + } + + curr_base_len = next_base_len + 1; + next_base_len = 0; + curr_base.swap( next_base ); + } + + } + else if ( is_tree_42 || is_tree_31 ) + { // "tree" simple reduce "42": 2->4->8->16->32->... // // .-------------------------------.-------------------------------. nr @@ -2666,236 +2970,6 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, // .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1 // 1 j nb - for (i = 1; i < nr; i++) { // layer by layer - // left - NodesBRD.SetValue(1, i+1, uv_el[i].node); - next_base.SetValue(++next_base_len, 1); - // right - NodesBRD.SetValue(nb, i+1, uv_er[i].node); - - next_par_u.SetValue(next_base_len, uv_el[i].u); - next_par_v.SetValue(next_base_len, uv_el[i].v); - - // to stop reducing, if number of nodes reaches nt - int delta = curr_base_len - nt; - - //double du = uv_er[i].u - uv_el[i].u; - //double dv = uv_er[i].v - uv_el[i].v; - - // to calculate normalized parameter, we must know number of points in next layer - int nb_four = (curr_base_len - 1) / 4; - int nb_next = nb_four*2 + (curr_base_len - nb_four*4); - if (nb_next < nt) nb_next = nt; - - for (j = 1; j + 4 <= curr_base_len && delta > 0; j += 4, delta -= 2) { - // add one "HH": nodes a,b,c,d,e and faces 1,2,3,4,5,6 - // - // .-----a-----b i + 1 - // |\ 5 | 6 /| - // | \ | / | - // | c--d--e | - // |1 |2 |3 |4 | - // | | | | | - // .--.--.--.--. i - // - // j j+2 j+4 - - double u,v; - - // a (i + 1, j + 2) - const SMDS_MeshNode* Na; - next_base_len++; - next_base.SetValue(next_base_len, curr_base.Value(j + 2)); - if (i + 1 == nr) { // top - Na = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else { - //double norm_par = double(next_base_len - 1)/double(nb_next - 1); - //u = uv_el[i].u + du * norm_par; - //v = uv_el[i].v + dv * norm_par; - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 2)); - //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 2)); - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Na1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Na1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na1); - Na = Na1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - - // b (i + 1, j + 4) - const SMDS_MeshNode* Nb; - next_base_len++; - next_base.SetValue(next_base_len, curr_base.Value(j + 4)); - if (i + 1 == nr) { // top - Nb = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 4 == curr_base_len) { // right - Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - //double norm_par = double(next_base_len - 1)/double(nb_next - 1); - //u = uv_el[i].u + du * norm_par; - //v = uv_el[i].v + dv * norm_par; - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 4)); - //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 4)); - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1); - Nb = Nb1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - - // c - u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 2)) / 2.0; - v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 2)) / 2.0; - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nc, geomFaceID, u, v); - - // d - u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 1)) / 2.0; - v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 1)) / 2.0; - P = S->Value(u,v); - SMDS_MeshNode* Nd = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nd, geomFaceID, u, v); - - // e - u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len)) / 2.0; - v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len)) / 2.0; - P = S->Value(u,v); - SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Ne, geomFaceID, u, v); - - // Faces - SMDS_MeshFace* F1 = 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 = 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 = 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 = 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 = 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 = myHelper->AddFace(Nd, Ne, Nb, Na); - if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID); - } - - // not reduced side elements (if any) - for (; j < curr_base_len; j++) { - // f (i + 1, j + 1) - const SMDS_MeshNode* Nf; - double u,v; - next_base.SetValue(++next_base_len, curr_base.Value(j + 1)); - if (i + 1 == nr) { // top - Nf = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 1 == curr_base_len) { // right - Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - //double norm_par = double(next_base_len - 1)/double(nb_next - 1); - //u = uv_el[i].u + du * norm_par; - //v = uv_el[i].v + dv * norm_par; - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 1)); - //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 1)); - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1); - Nf = Nf1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - SMDS_MeshFace* F1 = 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 "42" - else if (is_tree_31) { // "tree" simple reduce "31": 1->3->9->27->... // // .-----------------------------------------------------. nr @@ -2913,178 +2987,84 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, // .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1 // 1 j nb - for (i = 1; i < nr; i++) { // layer by layer - // left - NodesBRD.SetValue(1, i+1, uv_el[i].node); - next_base.SetValue(++next_base_len, 1); - // right - NodesBRD.SetValue(nb, i+1, uv_er[i].node); + PReduceFunction reduceFunction = & ( is_tree_42 ? reduce42 : reduce31 ); - next_par_u.SetValue(next_base_len, uv_el[i].u); - next_par_v.SetValue(next_base_len, uv_el[i].v); + const int reduce_grp_size = is_tree_42 ? 4 : 3; + for (i = 1; i < nr; i++) // layer by layer + { // to stop reducing, if number of nodes reaches nt int delta = curr_base_len - nt; // to calculate normalized parameter, we must know number of points in next layer - int nb_three = (curr_base_len - 1) / 3; - int nb_next = nb_three + (curr_base_len - nb_three*3); + int nb_reduce_groups = (curr_base_len - 1) / reduce_grp_size; + int nb_next = nb_reduce_groups * (reduce_grp_size-2) + (curr_base_len - nb_reduce_groups*reduce_grp_size); if (nb_next < nt) nb_next = nt; - for (j = 1; j + 3 <= curr_base_len && delta > 0; j += 3, delta -= 2) { - // add one "H": nodes b,c,e and faces 1,2,4,5 - // - // .---------b i + 1 - // |\ 5 /| - // | \ / | - // | c---e | - // |1 |2 |4 | - // | | | | - // .--.---.--. i - // - // j j+1 j+2 j+3 - - double u,v; - - // b (i + 1, j + 3) - const SMDS_MeshNode* Nb; - next_base_len++; - next_base.SetValue(next_base_len, curr_base.Value(j + 3)); - if (i + 1 == nr) { // top - Nb = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 3 == curr_base_len) { // right - Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1); - Nb = Nb1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - - // c and e - double u1 = (curr_par_u.Value(j) + next_par_u.Value(next_base_len - 1)) / 2.0; - double u2 = (curr_par_u.Value(j + 3) + next_par_u.Value(next_base_len)) / 2.0; - double u3 = (u2 - u1) / 3.0; - - double v1 = (curr_par_v.Value(j) + next_par_v.Value(next_base_len - 1)) / 2.0; - double v2 = (curr_par_v.Value(j + 3) + next_par_v.Value(next_base_len)) / 2.0; - double v3 = (v2 - v1) / 3.0; - - // c - u = u1 + u3; - v = v1 + v3; - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nc, geomFaceID, u, v); - - // e - u = u1 + u3 + u3; - v = v1 + v3 + v3; - P = S->Value(u,v); - SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Ne, geomFaceID, u, v); - - // Faces - SMDS_MeshFace* F1 = 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); + const double y = uv_el[ i ].normParam; + + if ( i + 1 == nr ) // top + { + next_base = uv_et; + } + else + { + next_base.clear(); + next_base.resize( nb_next, nullUVPtStruct ); + next_base.front() = uv_el[i]; + next_base.back() = uv_er[i]; + + // compute normalized param u + double du = 1. / ( nb_next - 1 ); + next_base[0].normParam = 0.; + for ( j = 1; j < nb_next; ++j ) + next_base[j].normParam = next_base[j-1].normParam + du; + } + uv[ UV_L ].SetCoord( next_base.front().u, next_base.front().v ); + uv[ UV_R ].SetCoord( next_base.back().u, next_base.back().v ); + + for (j = 0; j+reduce_grp_size < curr_base_len && delta > 0; j+=reduce_grp_size, delta-=2) + { + reduceFunction( curr_base, next_base, j, next_base_len, quad, uv, y, myHelper, S ); } // not reduced side elements (if any) - for (; j < curr_base_len; j++) { + for (; j < curr_base_len-1; j++) + { // f (i + 1, j + 1) - const SMDS_MeshNode* Nf; - double u,v; - next_base.SetValue(++next_base_len, curr_base.Value(j + 1)); - if (i + 1 == nr) { // top - Nf = uv_et[next_base_len - 1].node; - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf); - u = uv_et[next_base_len - 1].u; - v = uv_et[next_base_len - 1].v; - } - else if (j + 1 == curr_base_len) { // right - Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1); - u = uv_er[i].u; - v = uv_er[i].v; - } - else { - { - double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1; - int nearest_node_j = (int)rel; - rel -= nearest_node_j; - int ij = (i + 1 - 1) * nt + (nearest_node_j - 1); - double u1 = quad->uv_grid[ij].u; - double v1 = quad->uv_grid[ij].v; - double u2 = quad->uv_grid[ij + 1].u; - double v2 = quad->uv_grid[ij + 1].v; - double duj = (u2 - u1) * rel; - double dvj = (v2 - v1) * rel; - u = u1 + duj; - v = v1 + dvj; - } - gp_Pnt P = S->Value(u,v); - SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v); - NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1); - Nf = Nf1; - } - next_par_u.SetValue(next_base_len, u); - next_par_v.SetValue(next_base_len, v); - SMDS_MeshFace* F1 = 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; + const SMDS_MeshNode*& Nf = next_base[++next_base_len].node; + if ( !Nf ) + Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S ); + + myHelper->AddFace(curr_base[ j ].node, + curr_base[ j+1 ].node, + Nf, + next_base[ next_base_len-1 ].node); + } + curr_base_len = next_base_len + 1; next_base_len = 0; + curr_base.swap( next_base ); } - } // end "tree" simple reduce "31" - else if (is_lin_42) { + } // end "tree" simple reduce + + else if ( is_lin_42 || is_lin_31 ) { + // "linear" simple reduce "31": 2->6->10->14 + // + // .-----------------------------.-----------------------------. nr + // | \ / | \ / | + // | .---------. | .---------. | + // | | | | | | | + // .---------.---------.---------.---------.---------.---------. + // | / \ / \ | / \ / \ | + // | / .-----. \ | / .-----. \ | i + // | / | | \ | / | | \ | + // .-----.-----.-----.-----.-----.-----.-----.-----.-----.-----. + // | / / \ / \ \ | / / \ / \ \ | + // | / / .-. \ \ | / / .-. \ \ | + // | / / / \ \ \ | / / / \ \ \ | + // .--.----.---.-----.---.-----.-.--.----.---.-----.---.-----.-. 1 + // 1 j nb + // "linear" simple reduce "42": 4->8->12->16 // // .---------------.---------------.---------------.---------------. nr @@ -3120,7 +3100,8 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, if (remainder > 0) { nb_col++; } - int free_left = ((nt - 1) - nb_col * 2) / 2; + const int col_top_size = is_lin_42 ? 2 : 1; + int free_left = ((nt - 1) - nb_col * col_top_size) / 2; free_left += nr - 2; int free_middle = (nr - 2) * 2; if (remainder > 0 && nb_col == 1) { @@ -3137,17 +3118,12 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, //int free_left = 2; //int free_middle = 4; - for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) { // layer by layer - // left - NodesBRD.SetValue(1, i+1, uv_el[i].node); - next_base.SetValue(++next_base_len, 1); - // right - NodesBRD.SetValue(nb, i+1, uv_er[i].node); + PReduceFunction reduceFunction = & ( is_lin_42 ? reduce42 : reduce31 ); - // left - next_par_u.SetValue(next_base_len, uv_el[i].u); - next_par_v.SetValue(next_base_len, uv_el[i].v); + const int reduce_grp_size = is_lin_42 ? 4 : 3; + for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) // layer by layer + { // to calculate normalized parameter, we must know number of points in next layer int nb_next = curr_base_len - nb_col * 2; if (remainder > 0 && i > remainder / 2) @@ -3155,399 +3131,40 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, 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; + const double y = uv_el[ i ].normParam; - // 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 = 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 = 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 = 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 = 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 = 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 = myHelper->AddFace(Nd, Ne, Nb, Na); - if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID); - - j += 4; - - // 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; - } + if ( i + 1 == nr ) // top + { + next_base = uv_et; } - - // 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 "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; + 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 = 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); + 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++) { @@ -3556,104 +3173,10 @@ bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh, // 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; + // 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) { @@ -3668,108 +3191,41 @@ 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 = 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); + 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 = 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; + 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 "31" + + } // end "linear" simple reduce + else { + return false; } } // end Simple Reduce implementation