-// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE
//
-// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
-// This library is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
//
-// This library is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// Lesser General Public License for more details.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
-// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
// File : StdMeshers_Quadrangle_2D.cxx
// Author : Paul RASCLE, EDF
#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 <BRep_Tool.hxx>
#include <Geom_Surface.hxx>
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";
_compatibleHypothesis.push_back("QuadrangleParams");
_compatibleHypothesis.push_back("QuadranglePreference");
_compatibleHypothesis.push_back("TrianglePreference");
- myTool = 0;
+ _compatibleHypothesis.push_back("ViscousLayers2D");
}
//=============================================================================
*/
//=============================================================================
-bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
- const TopoDS_Shape& aShape)// throw (SALOME_Exception)
+bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
+ const TopoDS_Shape& aShape)
{
- // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
- //Unexpect aCatchSalomeException);
+ const TopoDS_Face& F = TopoDS::Face(aShape);
+ Handle(Geom_Surface) S = BRep_Tool::Surface(F);
SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
aMesh.GetSubMesh(aShape);
SMESH_MesherHelper helper (aMesh);
- myTool = &helper;
+ myHelper = &helper;
+
+ myProxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
+ if ( !myProxyMesh )
+ return false;
- _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
+ _quadraticMesh = myHelper->IsQuadraticSubMesh(aShape);
+ myNeedSmooth = false;
FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape);
std::auto_ptr<FaceQuadStruct> quadDeleter (quad); // to delete quad at exit from Compute()
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;
}
}
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;
}
}
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++) {
for (i = ilow; i < iup; i++) {
for (j = jlow; j < jup; j++) {
const SMDS_MeshNode *a, *b, *c, *d;
- a = quad->uv_grid[j * nbhoriz + i].node;
- b = quad->uv_grid[j * nbhoriz + i + 1].node;
+ a = quad->uv_grid[j * nbhoriz + i ].node;
+ b = quad->uv_grid[j * nbhoriz + i + 1].node;
c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
- d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
- SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
+ d = quad->uv_grid[(j + 1) * nbhoriz + i ].node;
+ SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
if (face) {
meshDS->SetMeshElementOnShape(face, geomFaceID);
}
}
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
//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 {
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);
}
}
}
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
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 {
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);
}
}
}
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
//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 {
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);
}
}
}
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
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 {
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);
}
}
}
}
+ if ( myNeedSmooth )
+ Smooth( quad );
+
bool isOk = true;
return isOk;
}
}
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<UVPtStruct>& UVPSleft = quad->side[0]->GetUVPtStruct(true,0);
/* vector<UVPtStruct>& UVPStop = */quad->side[1]->GetUVPtStruct(false,1);
/* vector<UVPtStruct>& UVPSright = */quad->side[2]->GetUVPtStruct(true,1);
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<TOP_SIDE, ignoreMediumNodes));
+ quad->side.push_back(new StdMeshers_FaceSide(F, *edgeIt, &aMesh, nbSides < TOP_SIDE,
+ ignoreMediumNodes, myProxyMesh));
}
- else if (nbEdgesInWire.front() > 4) { // more than 4 edges - try to unite some
+ else if (nbEdgesInWire.front() > 4) // more than 4 edges - try to unite some
+ {
list< TopoDS_Edge > sideEdges;
+ vector< int > degenSides;
while (!edges.empty()) {
sideEdges.clear();
sideEdges.splice(sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
sideEdges.splice(sideEdges.begin(), edges, --edges.end());
}
}
- quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh,
- nbSides<TOP_SIDE, ignoreMediumNodes));
+ if ( sideEdges.size() == 1 && BRep_Tool::Degenerated( sideEdges.front() ))
+ degenSides.push_back( nbSides );
+
+ quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, nbSides < TOP_SIDE,
+ ignoreMediumNodes, myProxyMesh));
++nbSides;
}
+ if ( !degenSides.empty() && nbSides - degenSides.size() == 4 )
+ {
+ myNeedSmooth = true;
+ for ( unsigned i = TOP_SIDE; i < quad->side.size(); ++i )
+ quad->side[i]->Reverse();
+
+ for ( int i = degenSides.size()-1; i > -1; --i )
+ {
+ StdMeshers_FaceSide* degenSide = quad->side[ degenSides[ i ]];
+ delete degenSide;
+ quad->side.erase( quad->side.begin() + degenSides[ i ] );
+ }
+ for ( unsigned i = TOP_SIDE; i < quad->side.size(); ++i )
+ quad->side[i]->Reverse();
+
+ nbSides -= degenSides.size();
+ }
// issue 20222. Try to unite only edges shared by two same faces
if (nbSides < 4) {
// delete found sides
}
}
quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh,
- nbSides<TOP_SIDE, ignoreMediumNodes));
+ nbSides < TOP_SIDE,
+ ignoreMediumNodes, myProxyMesh));
++nbSides;
}
}
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;
//return error("Can't find nodes on sides");
return error(COMPERR_BAD_INPUT_MESH);
- // nodes Id on "in" edges
- if (! quad->isEdgeOut[0]) {
- int j = 0;
- for (int i = 0; i < nbhoriz; i++) { // down
- int ij = j * nbhoriz + i;
- uv_grid[ij].node = uv_e0[i].node;
- }
+ if ( myNeedSmooth )
+ UpdateDegenUV( quad );
+
+ // copy data of face boundary
+ /*if (! quad->isEdgeOut[0])*/ {
+ const int j = 0;
+ for (int i = 0; i < nbhoriz; i++) // down
+ uv_grid[ j * nbhoriz + i ] = uv_e0[i];
}
- if (! quad->isEdgeOut[1]) {
- int i = nbhoriz - 1;
- for (int j = 0; j < nbvertic; j++) { // right
- int ij = j * nbhoriz + i;
- uv_grid[ij].node = uv_e1[j].node;
- }
+ /*if (! quad->isEdgeOut[1])*/ {
+ const int i = nbhoriz - 1;
+ for (int j = 0; j < nbvertic; j++) // right
+ uv_grid[ j * nbhoriz + i ] = uv_e1[j];
}
- if (! quad->isEdgeOut[2]) {
- int j = nbvertic - 1;
- for (int i = 0; i < nbhoriz; i++) { // up
- int ij = j * nbhoriz + i;
- uv_grid[ij].node = uv_e2[i].node;
- }
+ /*if (! quad->isEdgeOut[2])*/ {
+ const int j = nbvertic - 1;
+ for (int i = 0; i < nbhoriz; i++) // up
+ uv_grid[ j * nbhoriz + i ] = uv_e2[i];
}
- if (! quad->isEdgeOut[3]) {
+ /*if (! quad->isEdgeOut[3])*/ {
int i = 0;
- for (int j = 0; j < nbvertic; j++) { // left
- int ij = j * nbhoriz + i;
- uv_grid[ij].node = uv_e3[j].node;
- }
+ for (int j = 0; j < nbvertic; j++) // left
+ uv_grid[ j * nbhoriz + i ] = uv_e3[j];
}
- // normalized 2d values on grid
+ // normalized 2d parameters on grid
for (int i = 0; i < nbhoriz; i++) {
for (int j = 0; j < nbvertic; j++) {
int ij = j * nbhoriz + i;
}
// 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 "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
- gp_UV p0 = quad->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);
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);
const gp_UV& a0, const gp_UV& a1,
const gp_UV& a2, const gp_UV& a3)
{
- const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
- const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
- const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
- const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
+ // const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
+ // const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
+ // const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
+ // const vector<UVPtStruct>& 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);
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:"<<endl;
TColStd_SequenceOfReal npb, npr, npt, npl;
for (j=1; j<nl; j++) {
if (WisF) {
SMDS_MeshFace* F =
- myTool->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);
}
for (j=1; j<nr; j++) {
if (WisF) {
SMDS_MeshFace* F =
- myTool->AddFace(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);
}
for (j=1; j<nbv; j++) {
if (WisF) {
SMDS_MeshFace* F =
- myTool->AddFace(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);
}
for (i=1; i<nb; i++) {
if (WisF) {
SMDS_MeshFace* F =
- myTool->AddFace(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);
}
for (i=1; i<nb; i++) {
if (WisF) {
SMDS_MeshFace* F =
- myTool->AddFace(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);
}
for (i=1; i<nt; i++) {
if (WisF) {
SMDS_MeshFace* F =
- myTool->AddFace(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);
}
gp_Pnt d(theNode4->X(),theNode4->Y(),theNode4->Z());
SMDS_MeshFace* face;
if (a.Distance(c) > b.Distance(d)){
- face = myTool->AddFace(theNode2, theNode4 , theNode1);
+ face = myHelper->AddFace(theNode2, theNode4 , theNode1);
if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
- face = myTool->AddFace(theNode2, theNode3, theNode4);
+ face = myHelper->AddFace(theNode2, theNode3, theNode4);
if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
}
else{
- face = myTool->AddFace(theNode1, theNode2 ,theNode3);
+ face = myHelper->AddFace(theNode1, theNode2 ,theNode3);
if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
- face = myTool->AddFace(theNode1, theNode3, theNode4);
+ face = myHelper->AddFace(theNode1, theNode3, theNode4);
if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
}
}
+namespace
+{
+ enum uvPos { UV_A0, UV_A1, UV_A2, UV_A3, UV_B, UV_R, UV_T, UV_L, UV_SIZE };
+
+ inline SMDS_MeshNode* makeNode( UVPtStruct & uvPt,
+ const double y,
+ FaceQuadStruct* quad,
+ const gp_UV* UVs,
+ SMESH_MesherHelper* helper,
+ Handle(Geom_Surface) S)
+ {
+ const vector<UVPtStruct>& uv_eb = quad->side[BOTTOM_SIDE]->GetUVPtStruct();
+ const vector<UVPtStruct>& 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<UVPtStruct>& curr_base,
+ vector<UVPtStruct>& 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<UVPtStruct>& curr_base,
+ vector<UVPtStruct>& 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<UVPtStruct>& curr_base,
+ vector<UVPtStruct>& 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)
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)
//
// .-----.-----.-----.-----. .-----.-----.-----.-----.
// | / \ | / \ | | / \ | / \ |
int nb1 = nb;
int nr1 = nr;
int nt1 = nt;
- int nl1 = nl;
if (nr == nl) {
if (nb < nt) {
}
}
else if (nb == nt) {
- nl1 = nb; // and == nt
nr1 = nb; // and == nt
if (nl < nr) {
nt1 = nl;
}
// number of rows and columns
- int nrows = nr1 - 1; // and also == nl1 - 1
+ int nrows = nr1 - 1;
int ncol_top = nt1 - 1;
int ncol_bot = nb1 - 1;
- // maximum number of bottom elements for "tree" simple reduce 3->1
- int max_tree31 = ncol_top * pow(3.0, nrows);
- if (ncol_bot > max_tree31)
+ // number of rows needed to reduce ncol_bot to ncol_top using simple 3->1 "tree" (see below)
+ int nrows_tree31 = int( log( (double)(ncol_bot / ncol_top) ) / log((double) 3 )); // = log x base 3
+ if ( nrows < nrows_tree31 )
MultipleReduce = true;
}
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;
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++) {
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;
for (i=1; i<=dl; i++) {
for (j=1; j<nl; j++) {
SMDS_MeshFace* F =
- myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
- NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
+ myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
+ NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
}
}
for (i=1; i<=dr; i++) {
for (j=1; j<nr; j++) {
SMDS_MeshFace* F =
- myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
- NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
+ myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
+ NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
}
}
for (i=1; i<nb; i++) {
for (j=1; j<nbv; j++) {
SMDS_MeshFace* F =
- myTool->AddFace(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);
}
nr = quad->side[1]->NbPoints();
nt = quad->side[2]->NbPoints();
nl = quad->side[3]->NbPoints();
-
+
// number of rows and columns
int nrows = nr - 1; // and also == nl - 1
int ncol_top = nt - 1;
int ncol_bot = nb - 1;
int npair_top = ncol_top / 2;
// maximum number of bottom elements for "linear" simple reduce 4->2
- int max_lin = ncol_top + npair_top * 2 * nrows;
- // maximum number of bottom elements for "linear" simple reduce 4->2
+ int max_lin42 = ncol_top + npair_top * 2 * nrows;
+ // maximum number of bottom elements for "linear" simple reduce 3->1
int max_lin31 = ncol_top + ncol_top * 2 * nrows;
// maximum number of bottom elements for "tree" simple reduce 4->2
- int max_tree42 = npair_top * pow(2.0, nrows + 1);
- if (ncol_top > npair_top * 2) {
- int delta = ncol_bot - max_tree42;
- for (int irow = 1; irow < nrows; irow++) {
- int nfour = delta / 4;
- delta -= nfour * 2;
+ int max_tree42 = 0;
+ // number of rows needed to reduce ncol_bot to ncol_top using simple 4->2 "tree"
+ int nrows_tree42 = int( log( (double)(ncol_bot / ncol_top) )/log((double)2) ); // needed to avoid overflow at pow(2) while computing max_tree42
+ if (nrows_tree42 < nrows) {
+ max_tree42 = npair_top * pow(2.0, nrows + 1);
+ if ( ncol_top > npair_top * 2 ) {
+ int delta = ncol_bot - max_tree42;
+ for (int irow = 1; irow < nrows; irow++) {
+ int nfour = delta / 4;
+ delta -= nfour * 2;
+ }
+ if (delta <= (ncol_top - npair_top * 2))
+ max_tree42 = ncol_bot;
}
- if (delta <= (ncol_top - npair_top * 2))
- max_tree42 = ncol_bot;
}
// maximum number of bottom elements for "tree" simple reduce 3->1
//int max_tree31 = ncol_top * pow(3.0, nrows);
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;
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<UVPtStruct> 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<int> 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
// .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1
// 1 j nb
- for (i = 1; i < nr; i++) { // layer by layer
- // left
- NodesBRD.SetValue(1, i+1, uv_el[i].node);
- next_base.SetValue(++next_base_len, 1);
- // right
- NodesBRD.SetValue(nb, i+1, uv_er[i].node);
-
- next_par_u.SetValue(next_base_len, uv_el[i].u);
- next_par_v.SetValue(next_base_len, uv_el[i].v);
-
- // to stop reducing, if number of nodes reaches nt
- int delta = curr_base_len - nt;
-
- //double du = uv_er[i].u - uv_el[i].u;
- //double dv = uv_er[i].v - uv_el[i].v;
-
- // to calculate normalized parameter, we must know number of points in next layer
- int nb_four = (curr_base_len - 1) / 4;
- int nb_next = nb_four*2 + (curr_base_len - nb_four*4);
- if (nb_next < nt) nb_next = nt;
-
- for (j = 1; j + 4 <= curr_base_len && delta > 0; j += 4, delta -= 2) {
- // add one "HH": nodes a,b,c,d,e and faces 1,2,3,4,5,6
- //
- // .-----a-----b i + 1
- // |\ 5 | 6 /|
- // | \ | / |
- // | c--d--e |
- // |1 |2 |3 |4 |
- // | | | | |
- // .--.--.--.--. i
- //
- // j j+2 j+4
-
- double u,v;
-
- // a (i + 1, j + 2)
- const SMDS_MeshNode* Na;
- next_base_len++;
- next_base.SetValue(next_base_len, curr_base.Value(j + 2));
- if (i + 1 == nr) { // top
- Na = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else {
- //double norm_par = double(next_base_len - 1)/double(nb_next - 1);
- //u = uv_el[i].u + du * norm_par;
- //v = uv_el[i].v + dv * norm_par;
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 2));
- //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 2));
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Na1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Na1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na1);
- Na = Na1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
-
- // b (i + 1, j + 4)
- const SMDS_MeshNode* Nb;
- next_base_len++;
- next_base.SetValue(next_base_len, curr_base.Value(j + 4));
- if (i + 1 == nr) { // top
- Nb = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else if (j + 4 == curr_base_len) { // right
- Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
- u = uv_er[i].u;
- v = uv_er[i].v;
- }
- else {
- //double norm_par = double(next_base_len - 1)/double(nb_next - 1);
- //u = uv_el[i].u + du * norm_par;
- //v = uv_el[i].v + dv * norm_par;
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 4));
- //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 4));
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
- Nb = Nb1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
-
- // c
- u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 2)) / 2.0;
- v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 2)) / 2.0;
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
-
- // d
- u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 1)) / 2.0;
- v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 1)) / 2.0;
- P = S->Value(u,v);
- SMDS_MeshNode* Nd = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nd, geomFaceID, u, v);
-
- // e
- u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len)) / 2.0;
- v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len)) / 2.0;
- P = S->Value(u,v);
- SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
-
- // Faces
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
- NodesBRD.Value(curr_base.Value(j + 1), i),
- Nc,
- NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
-
- SMDS_MeshFace* F2 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
- NodesBRD.Value(curr_base.Value(j + 2), i),
- Nd, Nc);
- if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
-
- SMDS_MeshFace* F3 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
- NodesBRD.Value(curr_base.Value(j + 3), i),
- Ne, Nd);
- if (F3) meshDS->SetMeshElementOnShape(F3, geomFaceID);
-
- SMDS_MeshFace* F4 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 3), i),
- NodesBRD.Value(curr_base.Value(j + 4), i),
- Nb, Ne);
- if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
-
- SMDS_MeshFace* F5 = myTool->AddFace(Nc, Nd, Na,
- NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
- if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
-
- SMDS_MeshFace* F6 = myTool->AddFace(Nd, Ne, Nb, Na);
- if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID);
- }
-
- // not reduced side elements (if any)
- for (; j < curr_base_len; j++) {
- // f (i + 1, j + 1)
- const SMDS_MeshNode* Nf;
- double u,v;
- next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
- if (i + 1 == nr) { // top
- Nf = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else if (j + 1 == curr_base_len) { // right
- Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
- u = uv_er[i].u;
- v = uv_er[i].v;
- }
- else {
- //double norm_par = double(next_base_len - 1)/double(nb_next - 1);
- //u = uv_el[i].u + du * norm_par;
- //v = uv_el[i].v + dv * norm_par;
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 1));
- //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 1));
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
- Nf = Nf1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i),
- NodesBRD.Value(curr_base.Value(j + 1), i),
- NodesBRD.Value(next_base.Value(next_base_len), i + 1),
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
- }
-
- curr_base_len = next_base_len;
- curr_base = next_base;
- curr_par_u = next_par_u;
- curr_par_v = next_par_v;
- next_base_len = 0;
- }
- } // end "tree" simple reduce "42"
- else if (is_tree_31) {
// "tree" simple reduce "31": 1->3->9->27->...
//
// .-----------------------------------------------------. nr
// .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1
// 1 j nb
- for (i = 1; i < nr; i++) { // layer by layer
- // left
- NodesBRD.SetValue(1, i+1, uv_el[i].node);
- next_base.SetValue(++next_base_len, 1);
- // right
- NodesBRD.SetValue(nb, i+1, uv_er[i].node);
+ PReduceFunction reduceFunction = & ( is_tree_42 ? reduce42 : reduce31 );
- next_par_u.SetValue(next_base_len, uv_el[i].u);
- next_par_v.SetValue(next_base_len, uv_el[i].v);
+ const int reduce_grp_size = is_tree_42 ? 4 : 3;
+ for (i = 1; i < nr; i++) // layer by layer
+ {
// to stop reducing, if number of nodes reaches nt
int delta = curr_base_len - nt;
// to calculate normalized parameter, we must know number of points in next layer
- int nb_three = (curr_base_len - 1) / 3;
- int nb_next = nb_three + (curr_base_len - nb_three*3);
+ int nb_reduce_groups = (curr_base_len - 1) / reduce_grp_size;
+ int nb_next = nb_reduce_groups * (reduce_grp_size-2) + (curr_base_len - nb_reduce_groups*reduce_grp_size);
if (nb_next < nt) nb_next = nt;
- for (j = 1; j + 3 <= curr_base_len && delta > 0; j += 3, delta -= 2) {
- // add one "H": nodes b,c,e and faces 1,2,4,5
- //
- // .---------b i + 1
- // |\ 5 /|
- // | \ / |
- // | c---e |
- // |1 |2 |4 |
- // | | | |
- // .--.---.--. i
- //
- // j j+1 j+2 j+3
-
- double u,v;
-
- // b (i + 1, j + 3)
- const SMDS_MeshNode* Nb;
- next_base_len++;
- next_base.SetValue(next_base_len, curr_base.Value(j + 3));
- if (i + 1 == nr) { // top
- Nb = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else if (j + 3 == curr_base_len) { // right
- Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
- u = uv_er[i].u;
- v = uv_er[i].v;
- }
- else {
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
- Nb = Nb1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
-
- // c and e
- double u1 = (curr_par_u.Value(j) + next_par_u.Value(next_base_len - 1)) / 2.0;
- double u2 = (curr_par_u.Value(j + 3) + next_par_u.Value(next_base_len)) / 2.0;
- double u3 = (u2 - u1) / 3.0;
-
- double v1 = (curr_par_v.Value(j) + next_par_v.Value(next_base_len - 1)) / 2.0;
- double v2 = (curr_par_v.Value(j + 3) + next_par_v.Value(next_base_len)) / 2.0;
- double v3 = (v2 - v1) / 3.0;
-
- // c
- u = u1 + u3;
- v = v1 + v3;
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
-
- // e
- u = u1 + u3 + u3;
- v = v1 + v3 + v3;
- P = S->Value(u,v);
- SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
-
- // Faces
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
- NodesBRD.Value(curr_base.Value(j + 1), i),
- Nc,
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
-
- SMDS_MeshFace* F2 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
- NodesBRD.Value(curr_base.Value(j + 2), i),
- Ne, Nc);
- if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
-
- SMDS_MeshFace* F4 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
- NodesBRD.Value(curr_base.Value(j + 3), i),
- Nb, Ne);
- if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
-
- SMDS_MeshFace* F5 = myTool->AddFace(Nc, Ne, Nb,
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
+ const double y = uv_el[ i ].normParam;
+
+ if ( i + 1 == nr ) // top
+ {
+ next_base = uv_et;
+ }
+ else
+ {
+ next_base.clear();
+ next_base.resize( nb_next, nullUVPtStruct );
+ next_base.front() = uv_el[i];
+ next_base.back() = uv_er[i];
+
+ // compute normalized param u
+ double du = 1. / ( nb_next - 1 );
+ next_base[0].normParam = 0.;
+ for ( j = 1; j < nb_next; ++j )
+ next_base[j].normParam = next_base[j-1].normParam + du;
+ }
+ uv[ UV_L ].SetCoord( next_base.front().u, next_base.front().v );
+ uv[ UV_R ].SetCoord( next_base.back().u, next_base.back().v );
+
+ for (j = 0; j+reduce_grp_size < curr_base_len && delta > 0; j+=reduce_grp_size, delta-=2)
+ {
+ reduceFunction( curr_base, next_base, j, next_base_len, quad, uv, y, myHelper, S );
}
// not reduced side elements (if any)
- for (; j < curr_base_len; j++) {
+ for (; j < curr_base_len-1; j++)
+ {
// f (i + 1, j + 1)
- const SMDS_MeshNode* Nf;
- double u,v;
- next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
- if (i + 1 == nr) { // top
- Nf = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else if (j + 1 == curr_base_len) { // right
- Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
- u = uv_er[i].u;
- v = uv_er[i].v;
- }
- else {
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
- Nf = Nf1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i),
- NodesBRD.Value(curr_base.Value(j + 1), i),
- NodesBRD.Value(next_base.Value(next_base_len), i + 1),
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
- }
-
- curr_base_len = next_base_len;
- curr_base = next_base;
- curr_par_u = next_par_u;
- curr_par_v = next_par_v;
+ const SMDS_MeshNode*& Nf = next_base[++next_base_len].node;
+ if ( !Nf )
+ Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S );
+
+ myHelper->AddFace(curr_base[ j ].node,
+ curr_base[ j+1 ].node,
+ Nf,
+ next_base[ next_base_len-1 ].node);
+ }
+ curr_base_len = next_base_len + 1;
next_base_len = 0;
+ curr_base.swap( next_base );
}
- } // end "tree" simple reduce "31"
- else if (is_lin_42) {
+ } // end "tree" simple reduce
+
+ else if ( is_lin_42 || is_lin_31 ) {
+ // "linear" simple reduce "31": 2->6->10->14
+ //
+ // .-----------------------------.-----------------------------. nr
+ // | \ / | \ / |
+ // | .---------. | .---------. |
+ // | | | | | | |
+ // .---------.---------.---------.---------.---------.---------.
+ // | / \ / \ | / \ / \ |
+ // | / .-----. \ | / .-----. \ | i
+ // | / | | \ | / | | \ |
+ // .-----.-----.-----.-----.-----.-----.-----.-----.-----.-----.
+ // | / / \ / \ \ | / / \ / \ \ |
+ // | / / .-. \ \ | / / .-. \ \ |
+ // | / / / \ \ \ | / / / \ \ \ |
+ // .--.----.---.-----.---.-----.-.--.----.---.-----.---.-----.-. 1
+ // 1 j nb
+
// "linear" simple reduce "42": 4->8->12->16
//
// .---------------.---------------.---------------.---------------. nr
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) {
//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)
nb_next += 2;
if (nb_next < nt) nb_next = nt;
+ const double y = uv_el[ i ].normParam;
+
+ if ( i + 1 == nr ) // top
+ {
+ next_base = uv_et;
+ }
+ else
+ {
+ next_base.clear();
+ next_base.resize( nb_next, nullUVPtStruct );
+ next_base.front() = uv_el[i];
+ next_base.back() = uv_er[i];
+
+ // compute normalized param u
+ double du = 1. / ( nb_next - 1 );
+ next_base[0].normParam = 0.;
+ for ( j = 1; j < nb_next; ++j )
+ next_base[j].normParam = next_base[j-1].normParam + du;
+ }
+ uv[ UV_L ].SetCoord( next_base.front().u, next_base.front().v );
+ uv[ UV_R ].SetCoord( next_base.back().u, next_base.back().v );
+
// not reduced left elements
- for (j = 1; j <= free_left; j++) {
+ for (j = 0; j < free_left; j++)
+ {
// f (i + 1, j + 1)
- const SMDS_MeshNode* Nf;
- double u,v;
- next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
- if (i + 1 == nr) { // top
- Nf = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else {
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
- Nf = Nf1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i),
- NodesBRD.Value(curr_base.Value(j + 1), i),
- NodesBRD.Value(next_base.Value(next_base_len), i + 1),
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
+ const SMDS_MeshNode*& Nf = next_base[++next_base_len].node;
+ if ( !Nf )
+ Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S );
+
+ myHelper->AddFace(curr_base[ j ].node,
+ curr_base[ j+1 ].node,
+ Nf,
+ next_base[ next_base_len-1 ].node);
}
for (int icol = 1; icol <= nb_col; icol++) {
// stop short "column"
break;
- // add one "HH": nodes a,b,c,d,e and faces 1,2,3,4,5,6
- //
- // .-----a-----b i + 1
- // |\ 5 | 6 /|
- // | \ | / |
- // | c--d--e |
- // |1 |2 |3 |4 |
- // | | | | |
- // .--.--.--.--. i
- //
- // j j+2 j+4
-
- double u,v;
-
- // a (i + 1, j + 2)
- const SMDS_MeshNode* Na;
- next_base_len++;
- next_base.SetValue(next_base_len, curr_base.Value(j + 2));
- if (i + 1 == nr) { // top
- Na = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else {
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Na1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Na1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na1);
- Na = Na1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
-
- // b (i + 1, j + 4)
- const SMDS_MeshNode* Nb;
- next_base_len++;
- next_base.SetValue(next_base_len, curr_base.Value(j + 4));
- if (i + 1 == nr) { // top
- Nb = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else if (j + 4 == curr_base_len) { // right
- Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
- u = uv_er[i].u;
- v = uv_er[i].v;
- }
- else {
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
- Nb = Nb1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
-
- // c
- u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 2)) / 2.0;
- v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 2)) / 2.0;
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
-
- // d
- u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 1)) / 2.0;
- v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 1)) / 2.0;
- P = S->Value(u,v);
- SMDS_MeshNode* Nd = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nd, geomFaceID, u, v);
-
- // e
- u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len)) / 2.0;
- v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len)) / 2.0;
- P = S->Value(u,v);
- SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
-
- // Faces
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
- NodesBRD.Value(curr_base.Value(j + 1), i),
- Nc,
- NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
-
- SMDS_MeshFace* F2 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
- NodesBRD.Value(curr_base.Value(j + 2), i),
- Nd, Nc);
- if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
-
- SMDS_MeshFace* F3 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
- NodesBRD.Value(curr_base.Value(j + 3), i),
- Ne, Nd);
- if (F3) meshDS->SetMeshElementOnShape(F3, geomFaceID);
-
- SMDS_MeshFace* F4 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 3), i),
- NodesBRD.Value(curr_base.Value(j + 4), i),
- Nb, Ne);
- if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
-
- SMDS_MeshFace* F5 = myTool->AddFace(Nc, Nd, Na,
- NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
- if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
-
- SMDS_MeshFace* F6 = myTool->AddFace(Nd, Ne, Nb, Na);
- if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID);
-
- j += 4;
+ // add "H"
+ reduceFunction( curr_base, next_base, j, next_base_len, quad, uv, y, myHelper, S );
+
+ j += reduce_grp_size;
// not reduced middle elements
if (icol < nb_col) {
for (int imiddle = 1; imiddle <= free_add; imiddle++) {
// f (i + 1, j + imiddle)
- const SMDS_MeshNode* Nf;
- double u,v;
- next_base.SetValue(++next_base_len, curr_base.Value(j + imiddle));
- if (i + 1 == nr) { // top
- Nf = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else if (j + imiddle == curr_base_len) { // right
- Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
- u = uv_er[i].u;
- v = uv_er[i].v;
- }
- else {
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
- Nf = Nf1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j - 1 + imiddle), i),
- NodesBRD.Value(curr_base.Value(j + imiddle), i),
- NodesBRD.Value(next_base.Value(next_base_len), i + 1),
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
+ const SMDS_MeshNode*& Nf = next_base[++next_base_len].node;
+ if ( !Nf )
+ Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S );
+
+ myHelper->AddFace(curr_base[ j-1+imiddle ].node,
+ curr_base[ j +imiddle ].node,
+ Nf,
+ next_base[ next_base_len-1 ].node);
}
j += free_add;
}
}
// not reduced right elements
- for (; j < curr_base_len; j++) {
+ for (; j < curr_base_len-1; j++) {
// f (i + 1, j + 1)
- const SMDS_MeshNode* Nf;
- double u,v;
- next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
- if (i + 1 == nr) { // top
- Nf = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else if (j + 1 == curr_base_len) { // right
- Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
- u = uv_er[i].u;
- v = uv_er[i].v;
- }
- else {
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
- Nf = Nf1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i),
- NodesBRD.Value(curr_base.Value(j + 1), i),
- NodesBRD.Value(next_base.Value(next_base_len), i + 1),
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
- }
-
- curr_base_len = next_base_len;
- curr_base = next_base;
- curr_par_u = next_par_u;
- curr_par_v = next_par_v;
+ const SMDS_MeshNode*& Nf = next_base[++next_base_len].node;
+ if ( !Nf )
+ Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S );
+
+ myHelper->AddFace(curr_base[ j ].node,
+ curr_base[ j+1 ].node,
+ Nf,
+ next_base[ next_base_len-1 ].node);
+ }
+
+ curr_base_len = next_base_len + 1;
next_base_len = 0;
+ curr_base.swap( next_base );
}
- } // end "linear" simple reduce "42"
- else if (is_lin_31) {
- // "linear" simple reduce "31": 2->6->10->14
- //
- // .-----------------------------.-----------------------------. nr
- // | \ / | \ / |
- // | .---------. | .---------. |
- // | | | | | | |
- // .---------.---------.---------.---------.---------.---------.
- // | / \ / \ | / \ / \ |
- // | / .-----. \ | / .-----. \ | i
- // | / | | \ | / | | \ |
- // .-----.-----.-----.-----.-----.-----.-----.-----.-----.-----.
- // | / / \ / \ \ | / / \ / \ \ |
- // | / / .-. \ \ | / / .-. \ \ |
- // | / / / \ \ \ | / / / \ \ \ |
- // .--.----.---.-----.---.-----.-.--.----.---.-----.---.-----.-. 1
- // 1 j nb
- int delta_all = nb - nt;
- int delta_one_col = (nr - 1) * 2;
- int nb_col = delta_all / delta_one_col;
- int remainder = delta_all - nb_col * delta_one_col;
- if (remainder > 0) {
- nb_col++;
- }
- int free_left = ((nt - 1) - nb_col) / 2;
- free_left += nr - 2;
- int free_middle = (nr - 2) * 2;
- if (remainder > 0 && nb_col == 1) {
- int nb_rows_short_col = remainder / 2;
- int nb_rows_thrown = (nr - 1) - nb_rows_short_col;
- free_left -= nb_rows_thrown;
- }
+ } // end "linear" simple reduce
- for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) { // layer by layer
- // left
- NodesBRD.SetValue(1, i+1, uv_el[i].node);
- next_base.SetValue(++next_base_len, 1);
- // right
- NodesBRD.SetValue(nb, i+1, uv_er[i].node);
+ else {
+ return false;
+ }
+ } // end Simple Reduce implementation
- // left
- next_par_u.SetValue(next_base_len, uv_el[i].u);
- next_par_v.SetValue(next_base_len, uv_el[i].v);
+ bool isOk = true;
+ return isOk;
+}
- // to calculate normalized parameter, we must know number of points in next layer
- int nb_next = curr_base_len - nb_col * 2;
- if (remainder > 0 && i > remainder / 2)
- // take into account short "column"
- nb_next += 2;
- if (nb_next < nt) nb_next = nt;
+//================================================================================
+namespace // data for smoothing
+{
+ struct TSmoothNode;
+ // --------------------------------------------------------------------------------
+ /*!
+ * \brief Structure used to check validity of node position after smoothing.
+ * It holds two nodes connected to a smoothed node and belonging to
+ * one mesh face
+ */
+ struct TTriangle
+ {
+ TSmoothNode* _n1;
+ TSmoothNode* _n2;
+ TTriangle( TSmoothNode* n1=0, TSmoothNode* n2=0 ): _n1(n1), _n2(n2) {}
+
+ inline bool IsForward( gp_UV uv ) const;
+ };
+ // --------------------------------------------------------------------------------
+ /*!
+ * \brief Data of a smoothed node
+ */
+ struct TSmoothNode
+ {
+ gp_XY _uv;
+ vector< TTriangle > _triangles; // if empty, then node is not movable
+ };
+ // --------------------------------------------------------------------------------
+ inline bool TTriangle::IsForward( gp_UV uv ) const
+ {
+ gp_Vec2d v1( uv, _n1->_uv ), v2( uv, _n2->_uv );
+ double d = v1 ^ v2;
+ return d > 1e-100;
+ }
+}
- // not reduced left elements
- for (j = 1; j <= free_left; j++) {
- // f (i + 1, j + 1)
- const SMDS_MeshNode* Nf;
- double u,v;
- next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
- if (i + 1 == nr) { // top
- Nf = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else {
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
- Nf = Nf1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i),
- NodesBRD.Value(curr_base.Value(j + 1), i),
- NodesBRD.Value(next_base.Value(next_base_len), i + 1),
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
- }
+//================================================================================
+/*!
+ * \brief Set UV of nodes on degenerated VERTEXes in the middle of degenerated EDGE
+ *
+ * WARNING: this method must be called AFTER retrieving UVPtStruct's from quad
+ */
+//================================================================================
- for (int icol = 1; icol <= nb_col; icol++) {
+void StdMeshers_Quadrangle_2D::UpdateDegenUV(FaceQuadStruct* quad)
+{
+ for ( unsigned i = 0; i < quad->side.size(); ++i )
+ {
+ StdMeshers_FaceSide* side = quad->side[i];
+ const vector<UVPtStruct>& 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<UVPtStruct>& uvVec2 = side2->GetUVPtStruct();
+ int degenInd2 = -1;
+ if ( uvVec[ degenInd ].node == uvVec2[0].node )
+ degenInd2 = 0;
+ else if ( uvVec[ degenInd ].node == uvVec2.back().node )
+ degenInd2 = uvVec2.size() - 1;
+ else
+ throw SALOME_Exception( LOCALIZED( "Logical error" ));
- if (remainder > 0 && icol == nb_col && i > remainder / 2)
- // stop short "column"
- break;
+ // move UV in the middle
+ uvPtStruct& uv1 = const_cast<uvPtStruct&>( uvVec [ degenInd ]);
+ uvPtStruct& uv2 = const_cast<uvPtStruct&>( uvVec2[ degenInd2 ]);
+ uv1.u = uv2.u = 0.5 * ( uv1.u + uv2.u );
+ uv1.v = uv2.v = 0.5 * ( uv1.v + uv2.v );
+ }
+}
- // add one "H": nodes b,c,e and faces 1,2,4,5
- //
- // .---------b i + 1
- // |\ 5 /|
- // | \ / |
- // | c---e |
- // |1 |2 |4 |
- // | | | |
- // .--.---.--. i
- //
- // j j+1 j+2 j+3
-
- double u,v;
-
- // b (i + 1, j + 3)
- const SMDS_MeshNode* Nb;
- next_base_len++;
- next_base.SetValue(next_base_len, curr_base.Value(j + 3));
- if (i + 1 == nr) { // top
- Nb = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else if (j + 3 == curr_base_len) { // right
- Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
- u = uv_er[i].u;
- v = uv_er[i].v;
- }
- else {
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
- Nb = Nb1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
-
- // c and d
- double u1 = (curr_par_u.Value(j) + next_par_u.Value(next_base_len - 1)) / 2.0;
- double u2 = (curr_par_u.Value(j + 3) + next_par_u.Value(next_base_len)) / 2.0;
- double u3 = (u2 - u1) / 3.0;
-
- double v1 = (curr_par_v.Value(j) + next_par_v.Value(next_base_len - 1)) / 2.0;
- double v2 = (curr_par_v.Value(j + 3) + next_par_v.Value(next_base_len)) / 2.0;
- double v3 = (v2 - v1) / 3.0;
-
- // c
- u = u1 + u3;
- v = v1 + v3;
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
-
- // e
- u = u1 + u3 + u3;
- v = v1 + v3 + v3;
- P = S->Value(u,v);
- SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
-
- // Faces
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
- NodesBRD.Value(curr_base.Value(j + 1), i),
- Nc,
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
-
- SMDS_MeshFace* F2 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
- NodesBRD.Value(curr_base.Value(j + 2), i),
- Ne, Nc);
- if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
-
- SMDS_MeshFace* F4 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
- NodesBRD.Value(curr_base.Value(j + 3), i),
- Nb, Ne);
- if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
-
- SMDS_MeshFace* F5 = myTool->AddFace(Nc, Ne, Nb,
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
-
- j += 3;
+//================================================================================
+/*!
+ * \brief Perform smoothing of 2D elements on a FACE with ignored degenerated EDGE
+ */
+//================================================================================
- // not reduced middle elements
- if (icol < nb_col) {
- if (remainder > 0 && icol == nb_col - 1 && i > remainder / 2)
- // pass middle elements before stopped short "column"
- break;
+void StdMeshers_Quadrangle_2D::Smooth (FaceQuadStruct* quad)
+{
+ if ( !myNeedSmooth ) return;
- int free_add = free_middle;
- if (remainder > 0 && icol == nb_col - 1)
- // next "column" is short
- free_add -= (nr - 1) - (remainder / 2);
+ // Get nodes to smooth
- for (int imiddle = 1; imiddle <= free_add; imiddle++) {
- // f (i + 1, j + imiddle)
- const SMDS_MeshNode* Nf;
- double u,v;
- next_base.SetValue(++next_base_len, curr_base.Value(j + imiddle));
- if (i + 1 == nr) { // top
- Nf = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else if (j + imiddle == curr_base_len) { // right
- Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
- u = uv_er[i].u;
- v = uv_er[i].v;
- }
- else {
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
- Nf = Nf1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j - 1 + imiddle), i),
- NodesBRD.Value(curr_base.Value(j + imiddle), i),
- NodesBRD.Value(next_base.Value(next_base_len), i + 1),
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
- }
- j += free_add;
- }
- }
+ typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap;
+ TNo2SmooNoMap smooNoMap;
- // not reduced right elements
- for (; j < curr_base_len; j++) {
- // f (i + 1, j + 1)
- const SMDS_MeshNode* Nf;
- double u,v;
- next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
- if (i + 1 == nr) { // top
- Nf = uv_et[next_base_len - 1].node;
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
- u = uv_et[next_base_len - 1].u;
- v = uv_et[next_base_len - 1].v;
- }
- else if (j + 1 == curr_base_len) { // right
- Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
- u = uv_er[i].u;
- v = uv_er[i].v;
- }
- else {
- {
- double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
- int nearest_node_j = (int)rel;
- rel -= nearest_node_j;
- int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
- double u1 = quad->uv_grid[ij].u;
- double v1 = quad->uv_grid[ij].v;
- double u2 = quad->uv_grid[ij + 1].u;
- double v2 = quad->uv_grid[ij + 1].v;
- double duj = (u2 - u1) * rel;
- double dvj = (v2 - v1) * rel;
- u = u1 + duj;
- v = v1 + dvj;
- }
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
- NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
- Nf = Nf1;
- }
- next_par_u.SetValue(next_base_len, u);
- next_par_v.SetValue(next_base_len, v);
- SMDS_MeshFace* F1 = myTool->AddFace(NodesBRD.Value(curr_base.Value(j), i),
- NodesBRD.Value(curr_base.Value(j + 1), i),
- NodesBRD.Value(next_base.Value(next_base_len), i + 1),
- NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
- if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
- }
-
- curr_base_len = next_base_len;
- curr_base = next_base;
- curr_par_u = next_par_u;
- curr_par_v = next_par_v;
- next_base_len = 0;
- }
- } // end "linear" simple reduce "31"
- else {
+ const TopoDS_Face& geomFace = TopoDS::Face( myHelper->GetSubShape() );
+ SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
+ SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace );
+ SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes();
+ while ( nIt->more() ) // loop on nodes bound to a FACE
+ {
+ const SMDS_MeshNode* node = nIt->next();
+ TSmoothNode & sNode = smooNoMap[ node ];
+ sNode._uv = myHelper->GetNodeUV( geomFace, node );
+
+ // set sNode._triangles
+ SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face );
+ while ( fIt->more() )
+ {
+ const SMDS_MeshElement* face = fIt->next();
+ const int nbN = face->NbCornerNodes();
+ const int nInd = face->GetNodeIndex( node );
+ const int prevInd = myHelper->WrapIndex( nInd - 1, nbN );
+ const int nextInd = myHelper->WrapIndex( nInd + 1, nbN );
+ const SMDS_MeshNode* prevNode = face->GetNode( prevInd );
+ const SMDS_MeshNode* nextNode = face->GetNode( nextInd );
+ sNode._triangles.push_back( TTriangle( & smooNoMap[ prevNode ],
+ & smooNoMap[ nextNode ]));
}
- } // end Simple Reduce implementation
+ }
+ // set _uv of smooth nodes on FACE boundary
+ for ( unsigned i = 0; i < quad->side.size(); ++i )
+ {
+ const vector<UVPtStruct>& uvVec = quad->side[i]->GetUVPtStruct();
+ for ( unsigned j = 0; j < uvVec.size(); ++j )
+ {
+ TSmoothNode & sNode = smooNoMap[ uvVec[j].node ];
+ sNode._uv.SetCoord( uvVec[j].u, uvVec[j].v );
+ }
+ }
- bool isOk = true;
- return isOk;
+ // define refernce orientation in 2D
+ TNo2SmooNoMap::iterator n2sn = smooNoMap.begin();
+ for ( ; n2sn != smooNoMap.end(); ++n2sn )
+ if ( !n2sn->second._triangles.empty() )
+ break;
+ if ( n2sn == smooNoMap.end() ) return;
+ const TSmoothNode & sampleNode = n2sn->second;
+ const bool refForward = ( sampleNode._triangles[0].IsForward( sampleNode._uv ));
+
+ // Smoothing
+
+ for ( int iLoop = 0; iLoop < 5; ++iLoop )
+ {
+ for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn )
+ {
+ TSmoothNode& sNode = n2sn->second;
+ if ( sNode._triangles.empty() )
+ continue; // not movable node
+
+ // compute a new UV
+ gp_XY newUV (0,0);
+ for ( unsigned i = 0; i < sNode._triangles.size(); ++i )
+ newUV += sNode._triangles[i]._n1->_uv;
+ newUV /= sNode._triangles.size();
+
+ // check validity of the newUV
+ bool isValid = true;
+ for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i )
+ isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward );
+
+ if ( isValid )
+ sNode._uv = newUV;
+ }
+ }
+
+ // Set new XYZ to the smoothed nodes
+
+ Handle(Geom_Surface) surface = BRep_Tool::Surface( geomFace );
+
+ for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn )
+ {
+ TSmoothNode& sNode = n2sn->second;
+ if ( sNode._triangles.empty() )
+ continue; // not movable node
+
+ SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( n2sn->first );
+ gp_Pnt xyz = surface->Value( sNode._uv.X(), sNode._uv.Y() );
+ meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() );
+
+ // store the new UV
+ node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( sNode._uv.X(), sNode._uv.Y() )));
+ }
+
+ // Move medium nodes in quadratic mesh
+ if ( _quadraticMesh )
+ {
+ const TLinkNodeMap& links = myHelper->GetTLinkNodeMap();
+ TLinkNodeMap::const_iterator linkIt = links.begin();
+ for ( ; linkIt != links.end(); ++linkIt )
+ {
+ const SMESH_TLink& link = linkIt->first;
+ SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( linkIt->second );
+
+ if ( node->getshapeId() != myHelper->GetSubShapeID() )
+ continue; // medium node is on EDGE or VERTEX
+
+ gp_XY uv1 = myHelper->GetNodeUV( geomFace, link.node1(), node );
+ gp_XY uv2 = myHelper->GetNodeUV( geomFace, link.node2(), node );
+
+ gp_XY uv = myHelper->GetMiddleUV( surface, uv1, uv2 );
+ node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() )));
+
+ gp_Pnt xyz = surface->Value( uv.X(), uv.Y() );
+ meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() );
+ }
+ }
}
