X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=913f5a129400fcad984ebc8d673f5e71bcde7728;hp=14815906d2b93c90ff58445fee5197d0170e024c;hb=0cf919dde01c300a8f74d697198fba6538560c22;hpb=c3bf92bd87b770fd81631a3853f7f5bb1ac6a4e8 diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index 14815906d..913f5a129 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -31,6 +31,7 @@ using namespace std; #include "StdMeshers_Quadrangle_2D.hxx" #include "SMESH_Gen.hxx" #include "SMESH_Mesh.hxx" +#include "SMESH_subMesh.hxx" #include "SMDS_MeshElement.hxx" #include "SMDS_MeshNode.hxx" @@ -40,11 +41,14 @@ using namespace std; #include #include #include + #include #include #include -#include -#include +#include +#include + +#include #include #include @@ -58,13 +62,13 @@ using namespace std; */ //============================================================================= -StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D(int hypId, - int studyId, SMESH_Gen * gen):SMESH_2D_Algo(hypId, studyId, gen) +StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMESH_Gen* gen) + : SMESH_2D_Algo(hypId, studyId, gen) { - MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D"); - _name = "Quadrangle_2D"; - // _shapeType = TopAbs_FACE; - _shapeType = (1 << TopAbs_FACE); + MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D"); + _name = "Quadrangle_2D"; + // _shapeType = TopAbs_FACE; + _shapeType = (1 << TopAbs_FACE); } //============================================================================= @@ -75,7 +79,7 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D(int hypId, StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D() { - MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D"); + MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D"); } //============================================================================= @@ -89,14 +93,14 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis const TopoDS_Shape& aShape, SMESH_Hypothesis::Hypothesis_Status& aStatus) { - //MESSAGE("StdMeshers_Quadrangle_2D::CheckHypothesis"); + //MESSAGE("StdMeshers_Quadrangle_2D::CheckHypothesis"); - bool isOk = true; - aStatus = SMESH_Hypothesis::HYP_OK; + bool isOk = true; + aStatus = SMESH_Hypothesis::HYP_OK; - // nothing to check + // nothing to check - return isOk; + return isOk; } //============================================================================= @@ -105,71 +109,396 @@ 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) throw (SALOME_Exception) { Unexpect aCatch(SalomeException); - //MESSAGE("StdMeshers_Quadrangle_2D::Compute"); - SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); - SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape); - - FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape); - if (!quad) - return false; - - // --- compute 3D values on points, store points & quadrangles - - int nbdown = quad->nbPts[0]; - int nbright = quad->nbPts[1]; - int nbVertices = nbdown * nbright; - int nbQuad = (nbdown - 1) * (nbright - 1); - //SCRUTE(nbVertices); - //SCRUTE(nbQuad); - - // const TopoDS_Face& FF = TopoDS::Face(aShape); - // bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD); - // TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD)); - const TopoDS_Face & F = TopoDS::Face(aShape); - bool faceIsForward = (F.Orientation() == TopAbs_FORWARD); - Handle(Geom_Surface) S = BRep_Tool::Surface(F); - - for (int i = 1; i < nbdown - 1; i++) - for (int j = 1; j < nbright - 1; j++) // internal points - { - int ij = j * nbdown + i; - double u = quad->uv_grid[ij].u; - double v = quad->uv_grid[ij].v; - gp_Pnt P = S->Value(u, v); - SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(node, F); - quad->uv_grid[ij].node = node; -// Handle (SMDS_FacePosition) fpos -// = new SMDS_FacePosition(theSubMesh->GetId(),i,j); // easier than u,v -// node->SetPosition(fpos); - SMDS_FacePosition* fpos = - dynamic_cast(node->GetPosition().get()); - fpos->SetUParameter(i); - fpos->SetVParameter(j); - } - - // bool isQuadForward = ( faceIsForward == quad->isEdgeForward[0]); - for (int i = 0; i < nbdown - 1; i++) - for (int j = 0; j < nbright - 1; j++) // faces - { - const SMDS_MeshNode *a, *b, *c, *d; - a = quad->uv_grid[j * nbdown + i].node; - b = quad->uv_grid[j * nbdown + i + 1].node; - c = quad->uv_grid[(j + 1) * nbdown + i + 1].node; - d = quad->uv_grid[(j + 1) * nbdown + i].node; - // if (isQuadForward) faceId = meshDS->AddFace(a,b,c,d); - // else faceId = meshDS->AddFace(a,d,c,b); - SMDS_MeshFace * face = meshDS->AddFace(a, b, c, d); - meshDS->SetMeshElementOnShape(face, F); - } - - QuadDelete(quad); - bool isOk = true; - return isOk; + //MESSAGE("StdMeshers_Quadrangle_2D::Compute"); + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + aMesh.GetSubMesh(aShape); + + FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape); + if (!quad) + return false; + + // --- compute 3D values on points, store points & quadrangles + + int nbdown = quad->nbPts[0]; + int nbup = quad->nbPts[2]; +// bool isDownOut = (nbdown > nbup); +// bool isUpOut = (nbdown < nbup); + + int nbright = quad->nbPts[1]; + int nbleft = quad->nbPts[3]; +// bool isRightOut = (nbright > nbleft); +// bool isLeftOut = (nbright < nbleft); + + 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; + for (i = 1; i < nbhoriz - 1; i++) { + for (j = 1; j < nbvertic - 1; j++) { + int ij = j * nbhoriz + i; + double u = quad->uv_grid[ij].u; + double v = quad->uv_grid[ij].v; + gp_Pnt P = S->Value(u, v); + SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(node, F); + quad->uv_grid[ij].node = node; + SMDS_FacePosition* fpos = + dynamic_cast(node->GetPosition().get()); + fpos->SetUParameter(u); + fpos->SetVParameter(v); + } + } + + // mesh faces + + // [2] + // --.--.--.--.--.-- nbvertic + // | | ^ + // | | ^ + // [3] | | ^ j [1] + // | | ^ + // | | ^ + // ---.----.----.--- 0 + // 0 > > > > > > > > nbhoriz + // i + // [0] + + i = 0; + int ilow = 0; + int iup = nbhoriz - 1; + if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; } + + int jlow = 0; + int jup = nbvertic - 1; + if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; } + + // regular quadrangles + // bool isQuadForward = ( faceIsForward == quad->isEdgeForward[0]); + 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; + c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node; + d = quad->uv_grid[(j + 1) * nbhoriz + i].node; + // if (isQuadForward) faceId = meshDS->AddFace(a,b,c,d); + // else faceId = meshDS->AddFace(a,d,c,b); + SMDS_MeshFace * face = meshDS->AddFace(a, b, c, d); + meshDS->SetMeshElementOnShape(face, F); + } + } + + UVPtStruct *uv_e0 = quad->uv_edges[0]; + UVPtStruct *uv_e1 = quad->uv_edges[1]; + UVPtStruct *uv_e2 = quad->uv_edges[2]; + UVPtStruct *uv_e3 = quad->uv_edges[3]; + + double eps = Precision::Confusion(); + + // Boundary quadrangles + + if (quad->isEdgeOut[0]) { + // Down edge is out + // + // |___|___|___|___|___|___| + // | | | | | | | + // |___|___|___|___|___|___| + // | | | | | | | + // |___|___|___|___|___|___| __ first row of the regular grid + // . . . . . . . . . __ down edge nodes + // + // >->->->->->->->->->->->-> -- direction of processing + + int g = 0; // number of last processed node in the regular grid + + // number of last node of the down edge to be processed + int stop = nbdown - 1; + // if right edge is out, we will stop at a node, previous to the last one + if (quad->isEdgeOut[1]) stop--; + + // for each node of the down edge find nearest node + // in the first row of the regular grid and link them + for (i = 0; i < stop; i++) { + const SMDS_MeshNode *a, *b, *c, *d; + a = uv_e0[i].node; + b = uv_e0[i + 1].node; + gp_Pnt pb (b->X(), b->Y(), b->Z()); + + // find node c in the regular grid, which will be linked with node b + int near = g; + if (i == stop - 1) { + // right bound reached, link with the rightmost node + near = iup; + c = quad->uv_grid[nbhoriz + iup].node; + } else { + // find in the grid node c, nearest to the b + double mind = RealLast(); + for (int k = g; k <= iup; k++) { + + const SMDS_MeshNode *nk; + if (k < ilow) // this can be, if left edge is out + nk = uv_e3[1].node; // get node from the left edge + else + nk = quad->uv_grid[nbhoriz + k].node; // get one of middle nodes + + gp_Pnt pnk (nk->X(), nk->Y(), nk->Z()); + double dist = pb.Distance(pnk); + if (dist < mind - eps) { + c = nk; + near = k; + mind = dist; + } else { + break; + } + } + } + + if (near == g) { // make triangle + SMDS_MeshFace* face = meshDS->AddFace(a, b, c); + meshDS->SetMeshElementOnShape(face, F); + } else { // make quadrangle + if (near - 1 < ilow) + d = uv_e3[1].node; + else + d = quad->uv_grid[nbhoriz + near - 1].node; + SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + meshDS->SetMeshElementOnShape(face, F); + + // if node d is not at position g - make additional triangles + if (near - 1 > g) { + for (int k = near - 1; k > g; k--) { + c = quad->uv_grid[nbhoriz + k].node; + if (k - 1 < ilow) + d = uv_e3[1].node; + else + d = quad->uv_grid[nbhoriz + k - 1].node; + SMDS_MeshFace* face = meshDS->AddFace(a, c, d); + meshDS->SetMeshElementOnShape(face, F); + } + } + g = near; + } + } + } else { + if (quad->isEdgeOut[2]) { + // Up edge is out + // + // <-<-<-<-<-<-<-<-<-<-<-<-< -- direction of processing + // + // . . . . . . . . . __ up edge nodes + // ___ ___ ___ ___ ___ ___ __ first row of the regular grid + // | | | | | | | + // |___|___|___|___|___|___| + // | | | | | | | + // |___|___|___|___|___|___| + // | | | | | | | + + int g = nbhoriz - 1; // last processed node in the regular grid + + int stop = 0; + // if left edge is out, we will stop at a second node + if (quad->isEdgeOut[3]) stop++; + + // for each node of the up edge find nearest node + // in the first row of the regular grid and link them + for (i = nbup - 1; i > stop; i--) { + const SMDS_MeshNode *a, *b, *c, *d; + a = uv_e2[i].node; + b = uv_e2[i - 1].node; + gp_Pnt pb (b->X(), b->Y(), b->Z()); + + // find node c in the grid, which will be linked with node b + int near = g; + if (i == stop + 1) { // left bound reached, link with the leftmost node + c = quad->uv_grid[nbhoriz*(nbvertic - 2) + ilow].node; + near = ilow; + } else { + // find node c in the grid, nearest to the b + double mind = RealLast(); + for (int k = g; k >= ilow; k--) { + const SMDS_MeshNode *nk; + if (k > iup) + nk = uv_e1[nbright - 2].node; + else + nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node; + gp_Pnt pnk (nk->X(), nk->Y(), nk->Z()); + double dist = pb.Distance(pnk); + if (dist < mind - eps) { + c = nk; + near = k; + mind = dist; + } else { + break; + } + } + } + + if (near == g) { // make triangle + SMDS_MeshFace* face = meshDS->AddFace(a, b, c); + meshDS->SetMeshElementOnShape(face, F); + } else { // make quadrangle + if (near + 1 > iup) + d = uv_e1[nbright - 2].node; + else + d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node; + SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + meshDS->SetMeshElementOnShape(face, F); + + if (near + 1 < g) { // if d not is at g - make additional triangles + for (int k = near + 1; k < g; k++) { + c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node; + if (k + 1 > iup) + d = uv_e1[nbright - 2].node; + else + d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node; + SMDS_MeshFace* face = meshDS->AddFace(a, c, d); + meshDS->SetMeshElementOnShape(face, F); + } + } + g = near; + } + } + } + } + + // right or left boundary quadrangles + if (quad->isEdgeOut[1]) { +// MESSAGE("right edge is out"); + int g = 0; // last processed node in the grid + int stop = nbright - 1; + if (quad->isEdgeOut[2]) stop--; + for (i = 0; i < stop; i++) { + const SMDS_MeshNode *a, *b, *c, *d; + a = uv_e1[i].node; + b = uv_e1[i + 1].node; + gp_Pnt pb (b->X(), b->Y(), b->Z()); + + // find node c in the grid, nearest to the b + int near = g; + if (i == stop - 1) { // up bondary reached + c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node; + near = jup; + } else { + double mind = RealLast(); + for (int k = g; k <= jup; k++) { + const SMDS_MeshNode *nk; + if (k < jlow) + nk = uv_e0[nbdown - 2].node; + else + nk = quad->uv_grid[nbhoriz*(k + 1) - 2].node; + gp_Pnt pnk (nk->X(), nk->Y(), nk->Z()); + double dist = pb.Distance(pnk); + if (dist < mind - eps) { + c = nk; + near = k; + mind = dist; + } else { + break; + } + } + } + + if (near == g) { // make triangle + SMDS_MeshFace* face = meshDS->AddFace(a, b, c); + meshDS->SetMeshElementOnShape(face, F); + } else { // make quadrangle + if (near - 1 < jlow) + d = uv_e0[nbdown - 2].node; + else + d = quad->uv_grid[nbhoriz*near - 2].node; + SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + meshDS->SetMeshElementOnShape(face, F); + + if (near - 1 > g) { // if d not is at g - make additional triangles + for (int k = near - 1; k > g; k--) { + c = quad->uv_grid[nbhoriz*(k + 1) - 2].node; + if (k - 1 < jlow) + d = uv_e0[nbdown - 2].node; + else + d = quad->uv_grid[nbhoriz*k - 2].node; + SMDS_MeshFace* face = meshDS->AddFace(a, c, d); + meshDS->SetMeshElementOnShape(face, F); + } + } + g = near; + } + } + } else { + if (quad->isEdgeOut[3]) { +// MESSAGE("left edge is out"); + int g = nbvertic - 1; // last processed node in the grid + int stop = 0; + if (quad->isEdgeOut[0]) stop++; + for (i = nbleft - 1; i > stop; i--) { + const SMDS_MeshNode *a, *b, *c, *d; + a = uv_e3[i].node; + b = uv_e3[i - 1].node; + gp_Pnt pb (b->X(), b->Y(), b->Z()); + + // find node c in the grid, nearest to the b + int near = g; + if (i == stop + 1) { // down bondary reached + c = quad->uv_grid[nbhoriz*jlow + 1].node; + near = jlow; + } else { + double mind = RealLast(); + for (int k = g; k >= jlow; k--) { + const SMDS_MeshNode *nk; + if (k > jup) + nk = uv_e2[1].node; + else + nk = quad->uv_grid[nbhoriz*k + 1].node; + gp_Pnt pnk (nk->X(), nk->Y(), nk->Z()); + double dist = pb.Distance(pnk); + if (dist < mind - eps) { + c = nk; + near = k; + mind = dist; + } else { + break; + } + } + } + + if (near == g) { // make triangle + SMDS_MeshFace* face = meshDS->AddFace(a, b, c); + meshDS->SetMeshElementOnShape(face, F); + } else { // make quadrangle + if (near + 1 > jup) + d = uv_e2[1].node; + else + d = quad->uv_grid[nbhoriz*(near + 1) + 1].node; + SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + meshDS->SetMeshElementOnShape(face, F); + + if (near + 1 < g) { // if d not is at g - make additional triangles + for (int k = near + 1; k < g; k++) { + c = quad->uv_grid[nbhoriz*k + 1].node; + if (k + 1 > jup) + d = uv_e2[1].node; + else + d = quad->uv_grid[nbhoriz*(k + 1) + 1].node; + SMDS_MeshFace* face = meshDS->AddFace(a, c, d); + meshDS->SetMeshElementOnShape(face, F); + } + } + g = near; + } + } + } + } + + QuadDelete(quad); + bool isOk = true; + return isOk; } //============================================================================= @@ -178,82 +507,55 @@ bool StdMeshers_Quadrangle_2D::Compute(SMESH_Mesh & aMesh, */ //============================================================================= -FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute(SMESH_Mesh & aMesh, - const TopoDS_Shape & aShape)throw(SALOME_Exception) +FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute + (SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) throw(SALOME_Exception) { Unexpect aCatch(SalomeException); - //MESSAGE("StdMeshers_Quadrangle_2D::ComputeWithoutStore"); - - SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape); - - // const TopoDS_Face& FF = TopoDS::Face(aShape); - // bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD); - // TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD)); - const TopoDS_Face & F = TopoDS::Face(aShape); - bool faceIsForward = (F.Orientation() == TopAbs_FORWARD); - - // verify 1 wire only, with 4 edges, same number of points on opposite edges - - if (NumberOfWires(F) != 1) - { - MESSAGE("only 1 wire by face (quadrangles)"); - return 0; - //throw SALOME_Exception(LOCALIZED("only 1 wire by face (quadrangles)")); - } - // const TopoDS_Wire WW = BRepTools::OuterWire(F); - // TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD)); - const TopoDS_Wire & W = BRepTools::OuterWire(F); - BRepTools_WireExplorer wexp(W, F); - - FaceQuadStruct *quad = new FaceQuadStruct; - for (int i = 0; i < 4; i++) - quad->uv_edges[i] = 0; - quad->uv_grid = 0; - - int nbEdges = 0; - for (wexp.Init(W, F); wexp.More(); wexp.Next()) - { - // const TopoDS_Edge& EE = wexp.Current(); - // TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD)); - const TopoDS_Edge & E = wexp.Current(); - int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes(); - if (nbEdges < 4) - { - quad->edge[nbEdges] = E; - quad->nbPts[nbEdges] = nb + 2; // internal points + 2 extrema - } - nbEdges++; - } - - if (nbEdges != 4) - { - MESSAGE("face must have 4 edges /quadrangles"); - QuadDelete(quad); - return 0; - //throw SALOME_Exception(LOCALIZED("face must have 4 edges /quadrangles")); - } - - if (quad->nbPts[0] != quad->nbPts[2]) - { - MESSAGE("different point number-opposed edge"); - QuadDelete(quad); - return 0; - //throw SALOME_Exception(LOCALIZED("different point number-opposed edge")); - } - - if (quad->nbPts[1] != quad->nbPts[3]) - { - MESSAGE("different point number-opposed edge"); - QuadDelete(quad); - return 0; - //throw SALOME_Exception(LOCALIZED("different point number-opposed edge")); - } - - // set normalized grid on unit square in parametric domain - - SetNormalizedGrid(aMesh, F, quad); - - return quad; + + //SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape); + + const TopoDS_Face & F = TopoDS::Face(aShape); + + // verify 1 wire only, with 4 edges + + if (NumberOfWires(F) != 1) + { + INFOS("only 1 wire by face (quadrangles)"); + return 0; + } + const TopoDS_Wire& W = BRepTools::OuterWire(F); + BRepTools_WireExplorer wexp (W, F); + + FaceQuadStruct *quad = new FaceQuadStruct; + for (int i = 0; i < 4; i++) + quad->uv_edges[i] = 0; + quad->uv_grid = 0; + + int nbEdges = 0; + for (wexp.Init(W, F); wexp.More(); wexp.Next()) + { + const TopoDS_Edge& E = wexp.Current(); + int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes(); + if (nbEdges < 4) + { + quad->edge[nbEdges] = E; + quad->nbPts[nbEdges] = nb + 2; // internal points + 2 extrema + } + nbEdges++; + } + + if (nbEdges != 4) + { + INFOS("face must have 4 edges /quadrangles"); + QuadDelete(quad); + return 0; + } + + // set normalized grid on unit square in parametric domain + + SetNormalizedGrid(aMesh, F, quad); + + return quad; } //============================================================================= @@ -262,21 +564,21 @@ FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute(SMESH_Mesh & aMesh, */ //============================================================================= -void StdMeshers_Quadrangle_2D::QuadDelete(FaceQuadStruct * quad) +void StdMeshers_Quadrangle_2D::QuadDelete (FaceQuadStruct * quad) { - //MESSAGE("StdMeshers_Quadrangle_2D::QuadDelete"); - if (quad) - { - for (int i = 0; i < 4; i++) - { - if (quad->uv_edges[i]) - delete[]quad->uv_edges[i]; - quad->edge[i].Nullify(); - } - if (quad->uv_grid) - delete[]quad->uv_grid; - delete quad; - } + //MESSAGE("StdMeshers_Quadrangle_2D::QuadDelete"); + if (quad) + { + for (int i = 0; i < 4; i++) + { + if (quad->uv_edges[i]) + delete [] quad->uv_edges[i]; + quad->edge[i].Nullify(); + } + if (quad->uv_grid) + delete [] quad->uv_grid; + delete quad; + } } //============================================================================= @@ -285,320 +587,440 @@ void StdMeshers_Quadrangle_2D::QuadDelete(FaceQuadStruct * quad) */ //============================================================================= -void StdMeshers_Quadrangle_2D::SetNormalizedGrid(SMESH_Mesh & aMesh, - const TopoDS_Shape & aShape, FaceQuadStruct * quad) throw(SALOME_Exception) +void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, + const TopoDS_Shape& aShape, + FaceQuadStruct* quad) throw (SALOME_Exception) { Unexpect aCatch(SalomeException); - // Algorithme décrit dans "Génération automatique de maillages" - // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85 - // traitement dans le domaine paramétrique 2d u,v - // transport - projection sur le carré unité - - const TopoDS_Face & F = TopoDS::Face(aShape); - - // 1 --- find orientation of the 4 edges, by test on extrema - - // max min 0 x1 1 - // |<----north-2-------^ a3 -------------> a2 - // | | ^1 1^ - // west-3 east-1 =right | | - // | | ==> | | - // y0 | | y1 | | - // | | |0 0| - // v----south-0--------> a0 -------------> a1 - // min max 0 x0 1 - // =down - // - - Handle(Geom2d_Curve) c2d[4]; - gp_Pnt2d pf[4]; - gp_Pnt2d pl[4]; - for (int i = 0; i < 4; i++) - { - c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i], - F, quad->first[i], quad->last[i]); - pf[i] = c2d[i]->Value(quad->first[i]); - pl[i] = c2d[i]->Value(quad->last[i]); - quad->isEdgeForward[i] = false; - } - - double eps2d = 1.e-3; // *** utiliser plutot TopExp::CommonVertex, puis - // distances si piece fausse - int i = 0; - if ((pf[1].Distance(pl[0]) < eps2d) || (pl[1].Distance(pl[0]) < eps2d)) - { - quad->isEdgeForward[0] = true; - } - else - { - double tmp = quad->first[0]; - quad->first[0] = quad->last[0]; - quad->last[0] = tmp; - pf[0] = c2d[0]->Value(quad->first[0]); - pl[0] = c2d[0]->Value(quad->last[0]); - } - for (int i = 1; i < 4; i++) - { - quad->isEdgeForward[i] = (pf[i].Distance(pl[i - 1]) < eps2d); - if (!quad->isEdgeForward[i]) - { - double tmp = quad->first[i]; - quad->first[i] = quad->last[i]; - quad->last[i] = tmp; - pf[i] = c2d[i]->Value(quad->first[i]); - pl[i] = c2d[i]->Value(quad->last[i]); - //SCRUTE(pf[i].Distance(pl[i-1])); - ASSERT(pf[i].Distance(pl[i - 1]) < eps2d); - } - } - //SCRUTE(pf[0].Distance(pl[3])); - ASSERT(pf[0].Distance(pl[3]) < eps2d); - -// for (int i=0; i<4; i++) -// { -// SCRUTE(quad->isEdgeForward[i]); -// MESSAGE(" -first "<uv_edges[i] = LoadEdgePoints(aMesh, F, - quad->edge[i], quad->first[i], quad->last[i]); - - // quad->isEdgeForward[i]); - } - for (int i = 2; i < 4; i++) - { - quad->uv_edges[i] = LoadEdgePoints(aMesh, F, - quad->edge[i], quad->last[i], quad->first[i]); - - // !quad->isEdgeForward[i]); - } - - // 3 --- 2D normalized values on unit square [0..1][0..1] - - int nbdown = quad->nbPts[0]; - int nbright = quad->nbPts[1]; - quad->uv_grid = new UVPtStruct[nbright * nbdown]; - - UVPtStruct *uv_grid = quad->uv_grid; - UVPtStruct *uv_e0 = quad->uv_edges[0]; - UVPtStruct *uv_e1 = quad->uv_edges[1]; - UVPtStruct *uv_e2 = quad->uv_edges[2]; - UVPtStruct *uv_e3 = quad->uv_edges[3]; - gp_Pnt2d a0 = pf[0]; - gp_Pnt2d a1 = pf[1]; - gp_Pnt2d a2 = pf[2]; - gp_Pnt2d a3 = pf[3]; - - // nodes Id on edges - - int j = 0; - for (int i = 0; i < nbdown; i++) - { - int ij = j * nbdown + i; - uv_grid[ij].node = uv_e0[i].node; - } - i = nbdown - 1; - for (int j = 0; j < nbright; j++) - { - int ij = j * nbdown + i; - uv_grid[ij].node = uv_e1[j].node; - } - j = nbright - 1; - for (int i = 0; i < nbdown; i++) - { - int ij = j * nbdown + i; - uv_grid[ij].node = uv_e2[i].node; - } - i = 0; - for (int j = 0; j < nbright; j++) - { - int ij = j * nbdown + i; - uv_grid[ij].node = uv_e3[j].node; - } - - // normalized 2d values on grid - - for (int i = 0; i < nbdown; i++) - for (int j = 0; j < nbright; j++) - { - int ij = j * nbdown + i; - // --- droite i cste : x = x0 + y(x1-x0) - double x0 = uv_e0[i].normParam; // bas - sud - double x1 = uv_e2[i].normParam; // haut - nord - // --- droite j cste : y = y0 + x(y1-y0) - double y0 = uv_e3[j].normParam; // gauche-ouest - double y1 = uv_e1[j].normParam; // droite - est - // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0) - double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); - double y = y0 + x * (y1 - y0); - uv_grid[ij].x = x; - uv_grid[ij].y = y; - //MESSAGE("-xy-01 "<Value(param_0); - gp_Pnt2d p1 = c2d[1]->Value(param_1); - gp_Pnt2d p2 = c2d[2]->Value(param_2); - gp_Pnt2d p3 = c2d[3]->Value(param_3); - - double u = - (1 - y) * p0.X() + x * p1.X() + y * p2.X() + (1 - x) * p3.X(); - double v = - (1 - y) * p0.Y() + x * p1.Y() + y * p2.Y() + (1 - x) * p3.Y(); - - u -= (1 - x) * (1 - y) * a0.X() + x * (1 - y) * a1.X() + - x * y * a2.X() + (1 - x) * y * a3.X(); - v -= (1 - x) * (1 - y) * a0.Y() + x * (1 - y) * a1.Y() + - x * y * a2.Y() + (1 - x) * y * a3.Y(); - - uv_grid[ij].u = u; - uv_grid[ij].v = v; - - //MESSAGE("-uv- "< a2 + // | | ^1 1^ + // west-3 east-1 =right | | + // | | ==> | | + // y0 | | y1 | | + // | | |0 0| + // v----south-0--------> a0 -------------> a1 + // min max 0 x0 1 + // =down + // + + Handle(Geom2d_Curve) c2d[4]; + gp_Pnt2d pf[4]; + gp_Pnt2d pl[4]; + for (int i = 0; i < 4; i++) + { + c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i], F, + quad->first[i], quad->last[i]); + pf[i] = c2d[i]->Value(quad->first[i]); + pl[i] = c2d[i]->Value(quad->last[i]); + quad->isEdgeForward[i] = false; + } + + double l0f1 = pl[0].SquareDistance(pf[1]); + double l0l1 = pl[0].SquareDistance(pl[1]); + double f0f1 = pf[0].SquareDistance(pf[1]); + double f0l1 = pf[0].SquareDistance(pl[1]); + if ( Min( l0f1, l0l1 ) < Min ( f0f1, f0l1 )) + { + quad->isEdgeForward[0] = true; + } else { + double tmp = quad->first[0]; + quad->first[0] = quad->last[0]; + quad->last[0] = tmp; + pf[0] = c2d[0]->Value(quad->first[0]); + pl[0] = c2d[0]->Value(quad->last[0]); + } + for (int i = 1; i < 4; i++) + { + l0l1 = pl[i - 1].SquareDistance(pl[i]); + l0f1 = pl[i - 1].SquareDistance(pf[i]); + quad->isEdgeForward[i] = ( l0f1 < l0l1 ); + if (!quad->isEdgeForward[i]) + { + double tmp = quad->first[i]; + quad->first[i] = quad->last[i]; + quad->last[i] = tmp; + pf[i] = c2d[i]->Value(quad->first[i]); + pl[i] = c2d[i]->Value(quad->last[i]); + } + } + + // 2 --- load 2d edge points (u,v) with orientation and value on unit square + + bool loadOk = true; + for (int i = 0; i < 2; i++) + { + quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i], + quad->first[i], quad->last[i]); + if (!quad->uv_edges[i]) loadOk = false; + } + + for (int i = 2; i < 4; i++) + { + quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i], + quad->last[i], quad->first[i]); + if (!quad->uv_edges[i]) loadOk = false; + } + + if (!loadOk) + { + INFOS("StdMeshers_Quadrangle_2D::SetNormalizedGrid - LoadEdgePoints failed"); + QuadDelete( quad ); + quad = 0; + return; + } + // 3 --- 2D normalized values on unit square [0..1][0..1] + + int nbhoriz = Min(quad->nbPts[0], quad->nbPts[2]); + int nbvertic = Min(quad->nbPts[1], quad->nbPts[3]); + + quad->isEdgeOut[0] = (quad->nbPts[0] > quad->nbPts[2]); + quad->isEdgeOut[1] = (quad->nbPts[1] > quad->nbPts[3]); + quad->isEdgeOut[2] = (quad->nbPts[2] > quad->nbPts[0]); + quad->isEdgeOut[3] = (quad->nbPts[3] > quad->nbPts[1]); + + quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz]; + + UVPtStruct *uv_grid = quad->uv_grid; + UVPtStruct *uv_e0 = quad->uv_edges[0]; + UVPtStruct *uv_e1 = quad->uv_edges[1]; + UVPtStruct *uv_e2 = quad->uv_edges[2]; + UVPtStruct *uv_e3 = quad->uv_edges[3]; + + // 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 (! 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[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[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; + } + } + + // falsificate "out" edges + if (quad->isEdgeOut[0]) // down + uv_e0 = MakeEdgePoints + (aMesh, F, quad->edge[0], quad->first[0], quad->last[0], nbhoriz - 1); + else if (quad->isEdgeOut[2]) // up + uv_e2 = MakeEdgePoints + (aMesh, F, quad->edge[2], quad->last[2], quad->first[2], nbhoriz - 1); + + if (quad->isEdgeOut[1]) // right + uv_e1 = MakeEdgePoints + (aMesh, F, quad->edge[1], quad->first[1], quad->last[1], nbvertic - 1); + else if (quad->isEdgeOut[3]) // left + uv_e3 = MakeEdgePoints + (aMesh, F, quad->edge[3], quad->last[3], quad->first[3], nbvertic - 1); + + // normalized 2d values on grid + for (int i = 0; i < nbhoriz; i++) + { + for (int j = 0; j < nbvertic; j++) + { + int ij = j * nbhoriz + i; + // --- droite i cste : x = x0 + y(x1-x0) + double x0 = uv_e0[i].normParam; // bas - sud + double x1 = uv_e2[i].normParam; // haut - nord + // --- droite j cste : y = y0 + x(y1-y0) + double y0 = uv_e3[j].normParam; // gauche-ouest + double y1 = uv_e1[j].normParam; // droite - est + // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0) + double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); + double y = y0 + x * (y1 - y0); + uv_grid[ij].x = x; + uv_grid[ij].y = y; + //MESSAGE("-xy-01 "<Value(param_0); + gp_Pnt2d p1 = c2d[1]->Value(param_1); + gp_Pnt2d p2 = c2d[2]->Value(param_2); + gp_Pnt2d p3 = c2d[3]->Value(param_3); + + double u = (1 - y) * p0.X() + x * p1.X() + y * p2.X() + (1 - x) * p3.X(); + double v = (1 - y) * p0.Y() + x * p1.Y() + y * p2.Y() + (1 - x) * p3.Y(); + + u -= (1 - x) * (1 - y) * a0.X() + x * (1 - y) * a1.X() + + x * y * a2.X() + (1 - x) * y * a3.X(); + v -= (1 - x) * (1 - y) * a0.Y() + x * (1 - y) * a1.Y() + + x * y * a2.Y() + (1 - x) * y * a3.Y(); + + uv_grid[ij].u = u; + uv_grid[ij].v = v; + } + } } //============================================================================= /*! - * + * LoadEdgePoints */ //============================================================================= +UVPtStruct* StdMeshers_Quadrangle_2D::LoadEdgePoints (SMESH_Mesh & aMesh, + const TopoDS_Face& F, + const TopoDS_Edge& E, + double first, double last) +// bool isForward) +{ + //MESSAGE("StdMeshers_Quadrangle_2D::LoadEdgePoints"); + + //SMDS_Mesh* meshDS = aMesh.GetMeshDS(); + + // --- IDNodes of first and last Vertex + + TopoDS_Vertex VFirst, VLast; + TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l + + ASSERT(!VFirst.IsNull()); + SMDS_NodeIteratorPtr lid = aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes(); + if (!lid->more()) + { + MESSAGE ( "NO NODE BUILT ON VERTEX" ); + return 0; + } + const SMDS_MeshNode* idFirst = lid->next(); + + ASSERT(!VLast.IsNull()); + lid = aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes(); + if (!lid->more()) + { + MESSAGE ( "NO NODE BUILT ON VERTEX" ); + return 0; + } + const SMDS_MeshNode* idLast = lid->next(); + + // --- edge internal IDNodes (relies on good order storage, not checked) + + map params; + SMDS_NodeIteratorPtr ite = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes(); + + while(ite->more()) + { + const SMDS_MeshNode* node = ite->next(); + const SMDS_EdgePosition* epos = + static_cast(node->GetPosition().get()); + double param = epos->GetUParameter(); + params[param] = node; + } + + int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes(); + if (nbPoints != params.size()) + { + MESSAGE( "BAD NODE ON EDGE POSITIONS" ); + return 0; + } + UVPtStruct* uvslf = new UVPtStruct[nbPoints + 2]; + + double f, l; + Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l); + + bool isForward = (((l - f) * (last - first)) > 0); + double paramin = 0; + double paramax = 0; + if (isForward) + { + paramin = f; + paramax = l; + gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward + uvslf[0].x = p.X(); + uvslf[0].y = p.Y(); + uvslf[0].param = f; + uvslf[0].node = idFirst; + //MESSAGE("__ f "<::iterator itp = params.begin(); + for (int i = 1; i <= nbPoints; i++) // nbPoints internal + { + double param = (*itp).first; + gp_Pnt2d p = C2d->Value(param); + uvslf[i].x = p.X(); + uvslf[i].y = p.Y(); + uvslf[i].param = param; + uvslf[i].node = (*itp).second; + //MESSAGE("__ "<Value(l); // last point = Vertex Reversed + uvslf[nbPoints + 1].x = p.X(); + uvslf[nbPoints + 1].y = p.Y(); + uvslf[nbPoints + 1].param = l; + uvslf[nbPoints + 1].node = idLast; + //MESSAGE("__ l "<Value(l); // first point = Vertex Reversed + uvslf[0].x = p.X(); + uvslf[0].y = p.Y(); + uvslf[0].param = l; + uvslf[0].node = idLast; + //MESSAGE("__ l "<::reverse_iterator itp = params.rbegin(); + + for (int j = nbPoints; j >= 1; j--) // nbPoints internal + { + double param = (*itp).first; + int i = nbPoints + 1 - j; + gp_Pnt2d p = C2d->Value(param); + uvslf[i].x = p.X(); + uvslf[i].y = p.Y(); + uvslf[i].param = param; + uvslf[i].node = (*itp).second; + //MESSAGE("__ "<Value(f); // last point = Vertex Forward + uvslf[nbPoints + 1].x = p.X(); + uvslf[nbPoints + 1].y = p.Y(); + uvslf[nbPoints + 1].param = f; + uvslf[nbPoints + 1].node = idFirst; + //MESSAGE("__ f "<GetSubMeshDS()->GetNodes(); - const SMDS_MeshNode * idFirst = lid->next(); - - ASSERT(!VLast.IsNull()); - lid=aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes(); - const SMDS_MeshNode * idLast = lid->next(); - - // --- edge internal IDNodes (relies on good order storage, not checked) - - int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes(); - //SCRUTE(nbPoints); - UVPtStruct *uvslf = new UVPtStruct[nbPoints + 2]; - - double f, l; - Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l); - - map params; - SMDS_NodeIteratorPtr ite= aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes(); - - while(ite->more()) - { - const SMDS_MeshNode * node = ite->next(); - const SMDS_EdgePosition* epos = - static_cast(node->GetPosition().get()); - double param = epos->GetUParameter(); - params[param] = node; - } - - bool isForward = (((l - f) * (last - first)) > 0); - double paramin = 0; - double paramax = 0; - if (isForward) - { - paramin = f; - paramax = l; - gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward - uvslf[0].x = p.X(); - uvslf[0].y = p.Y(); - uvslf[0].param = f; - uvslf[0].node = idFirst; - //MESSAGE("__ f "<::iterator itp = params.begin(); - for (int i = 1; i <= nbPoints; i++) // nbPoints internal - { - double param = (*itp).first; - gp_Pnt2d p = C2d->Value(param); - uvslf[i].x = p.X(); - uvslf[i].y = p.Y(); - uvslf[i].param = param; - uvslf[i].node = (*itp).second; - //MESSAGE("__ "<Value(l); // last point = Vertex Reversed - uvslf[nbPoints + 1].x = p.X(); - uvslf[nbPoints + 1].y = p.Y(); - uvslf[nbPoints + 1].param = l; - uvslf[nbPoints + 1].node = idLast; - //MESSAGE("__ l "<Value(l); // first point = Vertex Reversed - uvslf[0].x = p.X(); - uvslf[0].y = p.Y(); - uvslf[0].param = l; - uvslf[0].node = idLast; - //MESSAGE("__ l "<::reverse_iterator itp = params.rbegin(); - for (int j = nbPoints; j >= 1; j--) // nbPoints internal - { - double param = (*itp).first; - int i = nbPoints + 1 - j; - gp_Pnt2d p = C2d->Value(param); - uvslf[i].x = p.X(); - uvslf[i].y = p.Y(); - uvslf[i].param = param; - uvslf[i].node = (*itp).second; - //MESSAGE("__ "<Value(f); // last point = Vertex Forward - uvslf[nbPoints + 1].x = p.X(); - uvslf[nbPoints + 1].y = p.Y(); - uvslf[nbPoints + 1].param = f; - uvslf[nbPoints + 1].node = idFirst; - //MESSAGE("__ f "< params; + + // --- edge internal points + double fi, li; + Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, fi, li); + if (!Curve.IsNull()) { + try { + GeomAdaptor_Curve C3d (Curve); + double length = EdgeLength(E); + double eltSize = length / nb_segm; + GCPnts_UniformAbscissa Discret (C3d, eltSize, fi, li); + if (!Discret.IsDone()) return false; + int NbPoints = Discret.NbPoints(); + for (int i = 1; i <= NbPoints; i++) { + double param = Discret.Parameter(i); + params.push_back(param); + } + } + catch (Standard_Failure) { + return 0; + } + } + else + { + // Edge is a degenerated Edge + BRep_Tool::Range(E, fi, li); + double du = (li - fi) / nb_segm; + for (int i = 1; i <= nb_segm + 1; i++) + { + double param = fi + (i - 1) * du; + params.push_back(param); + } + } + + double f, l; + Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l); + ASSERT(f != l); + + bool isForward = (((l - f) * (last - first)) > 0); + if (isForward) { + list::iterator itU = params.begin(); + for (int i = 0; i <= nb_segm; i++) // nbPoints internal + { + double param = *itU; + gp_Pnt2d p = C2d->Value(param); + uvslf[i].x = p.X(); + uvslf[i].y = p.Y(); + uvslf[i].param = param; + uvslf[i].normParam = (param - f) / (l - f); + itU++; + } + } else { + list::reverse_iterator itU = params.rbegin(); + for (int j = nb_segm; j >= 0; j--) // nbPoints internal + { + double param = *itU; + int i = nb_segm - j; + gp_Pnt2d p = C2d->Value(param); + uvslf[i].x = p.X(); + uvslf[i].y = p.Y(); + uvslf[i].param = param; + uvslf[i].normParam = (param - l) / (f - l); + itU++; + } + } + + return uvslf; } + //============================================================================= /*! *