X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=b035da0ddd737e2f1ab7be8605caa5721fceedd9;hp=913f5a129400fcad984ebc8d673f5e71bcde7728;hb=57b43b4d010e2d0a1529d3c131bbb9d416e63258;hpb=e90b1a5ede07f874fbfe9ba40319d15c5be6a489 diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index 913f5a129..b035da0dd 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -47,14 +47,26 @@ using namespace std; #include #include #include +#include #include #include #include +#include +#include #include "utilities.h" #include "Utils_ExceptHandlers.hxx" +#ifndef StdMeshers_Array2OfNode_HeaderFile +#define StdMeshers_Array2OfNode_HeaderFile +typedef const SMDS_MeshNode* SMDS_MeshNodePtr; +#include +DEFINE_BASECOLLECTION (StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr) +DEFINE_ARRAY2(StdMeshers_Array2OfNode, + StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr) +#endif + //============================================================================= /*! @@ -67,8 +79,9 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMES { MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D"); _name = "Quadrangle_2D"; - // _shapeType = TopAbs_FACE; _shapeType = (1 << TopAbs_FACE); + _compatibleHypothesis.push_back("QuadranglePreference"); + myTool = 0; } //============================================================================= @@ -80,6 +93,8 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMES StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D() { MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D"); + if ( myTool ) + delete myTool; } //============================================================================= @@ -89,16 +104,16 @@ StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D() //============================================================================= bool StdMeshers_Quadrangle_2D::CheckHypothesis - (SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape, + (SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, SMESH_Hypothesis::Hypothesis_Status& aStatus) { - //MESSAGE("StdMeshers_Quadrangle_2D::CheckHypothesis"); - bool isOk = true; aStatus = SMESH_Hypothesis::HYP_OK; - // nothing to check + // there is only one compatible Hypothesis so far + const list &hyps = GetUsedHypothesis(aMesh, aShape); + myQuadranglePreference = hyps.size() > 0; return isOk; } @@ -117,21 +132,48 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); aMesh.GetSubMesh(aShape); - FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape); - if (!quad) + if ( !myTool ) + myTool = new SMESH_MesherHelper(aMesh); + _quadraticMesh = myTool->IsQuadraticSubMesh(aShape); + + //FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape); + FaceQuadStruct* quad = CheckNbEdges(aMesh, aShape); + + if (!quad) { + delete myTool; myTool = 0; return false; + } + + if(myQuadranglePreference) { + int n1 = quad->nbPts[0]; + int n2 = quad->nbPts[1]; + int n3 = quad->nbPts[2]; + int n4 = quad->nbPts[3]; + int nfull = n1+n2+n3+n4; + int ntmp = nfull/2; + ntmp = ntmp*2; + if( nfull==ntmp && ( (n1!=n3) || (n2!=n4) ) ) { + // special path for using only quandrangle faces + bool ok = ComputeQuadPref(aMesh, aShape, quad); + delete myTool; myTool = 0; + return ok; + } + } + + // set normalized grid on unit square in parametric domain + SetNormalizedGrid(aMesh, aShape, quad); + if (!quad) { + delete myTool; myTool = 0; + 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); @@ -140,7 +182,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, Handle(Geom_Surface) S = BRep_Tool::Surface(F); // internal mesh nodes - int i, j; + int i, j, geomFaceID = meshDS->ShapeToIndex( F ); for (i = 1; i < nbhoriz - 1; i++) { for (j = 1; j < nbvertic - 1; j++) { int ij = j * nbhoriz + i; @@ -148,15 +190,11 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, 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); + meshDS->SetNodeOnFace(node, geomFaceID, u, v); quad->uv_grid[ij].node = node; - SMDS_FacePosition* fpos = - dynamic_cast(node->GetPosition().get()); - fpos->SetUParameter(u); - fpos->SetVParameter(v); } } - + // mesh faces // [2] @@ -170,18 +208,17 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, // 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; @@ -189,13 +226,12 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, b = quad->uv_grid[j * nbhoriz + i + 1].node; c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node; d = quad->uv_grid[(j + 1) * nbhoriz + i].node; - // 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); + //SMDS_MeshFace * face = meshDS->AddFace(a, b, c, d); + SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + meshDS->SetMeshElementOnShape(face, geomFaceID); } } - + UVPtStruct *uv_e0 = quad->uv_edges[0]; UVPtStruct *uv_e1 = quad->uv_edges[1]; UVPtStruct *uv_e2 = quad->uv_edges[2]; @@ -204,7 +240,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, double eps = Precision::Confusion(); // Boundary quadrangles - + if (quad->isEdgeOut[0]) { // Down edge is out // @@ -216,14 +252,14 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, // . . . . . . . . . __ 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++) { @@ -231,18 +267,19 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, 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 { + } + 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 @@ -262,15 +299,18 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = meshDS->AddFace(a, b, c); - meshDS->SetMeshElementOnShape(face, F); - } else { // make quadrangle + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c); + SMDS_MeshFace* face = myTool->AddFace(a, b, c); + meshDS->SetMeshElementOnShape(face, geomFaceID); + } + 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); + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + meshDS->SetMeshElementOnShape(face, geomFaceID); // if node d is not at position g - make additional triangles if (near - 1 > g) { @@ -280,8 +320,9 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, d = uv_e3[1].node; else d = quad->uv_grid[nbhoriz + k - 1].node; - SMDS_MeshFace* face = meshDS->AddFace(a, c, d); - meshDS->SetMeshElementOnShape(face, F); + //SMDS_MeshFace* face = meshDS->AddFace(a, c, d); + SMDS_MeshFace* face = myTool->AddFace(a, c, d); + meshDS->SetMeshElementOnShape(face, geomFaceID); } } g = near; @@ -342,15 +383,18 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = meshDS->AddFace(a, b, c); - meshDS->SetMeshElementOnShape(face, F); - } else { // make quadrangle + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c); + SMDS_MeshFace* face = myTool->AddFace(a, b, c); + meshDS->SetMeshElementOnShape(face, geomFaceID); + } + 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); + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + meshDS->SetMeshElementOnShape(face, geomFaceID); if (near + 1 < g) { // if d not is at g - make additional triangles for (int k = near + 1; k < g; k++) { @@ -359,8 +403,9 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, d = uv_e1[nbright - 2].node; else d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node; - SMDS_MeshFace* face = meshDS->AddFace(a, c, d); - meshDS->SetMeshElementOnShape(face, F); + //SMDS_MeshFace* face = meshDS->AddFace(a, c, d); + SMDS_MeshFace* face = myTool->AddFace(a, c, d); + meshDS->SetMeshElementOnShape(face, geomFaceID); } } g = near; @@ -407,15 +452,18 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = meshDS->AddFace(a, b, c); - meshDS->SetMeshElementOnShape(face, F); - } else { // make quadrangle + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c); + SMDS_MeshFace* face = myTool->AddFace(a, b, c); + meshDS->SetMeshElementOnShape(face, geomFaceID); + } + 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); + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + meshDS->SetMeshElementOnShape(face, geomFaceID); if (near - 1 > g) { // if d not is at g - make additional triangles for (int k = near - 1; k > g; k--) { @@ -424,8 +472,9 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, d = uv_e0[nbdown - 2].node; else d = quad->uv_grid[nbhoriz*k - 2].node; - SMDS_MeshFace* face = meshDS->AddFace(a, c, d); - meshDS->SetMeshElementOnShape(face, F); + //SMDS_MeshFace* face = meshDS->AddFace(a, c, d); + SMDS_MeshFace* face = myTool->AddFace(a, c, d); + meshDS->SetMeshElementOnShape(face, geomFaceID); } } g = near; @@ -469,15 +518,18 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = meshDS->AddFace(a, b, c); - meshDS->SetMeshElementOnShape(face, F); - } else { // make quadrangle + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c); + SMDS_MeshFace* face = myTool->AddFace(a, b, c); + meshDS->SetMeshElementOnShape(face, geomFaceID); + } + 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); + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + meshDS->SetMeshElementOnShape(face, geomFaceID); if (near + 1 < g) { // if d not is at g - make additional triangles for (int k = near + 1; k < g; k++) { @@ -486,8 +538,9 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, d = uv_e2[1].node; else d = quad->uv_grid[nbhoriz*(k + 1) + 1].node; - SMDS_MeshFace* face = meshDS->AddFace(a, c, d); - meshDS->SetMeshElementOnShape(face, F); + //SMDS_MeshFace* face = meshDS->AddFace(a, c, d); + SMDS_MeshFace* face = myTool->AddFace(a, c, d); + meshDS->SetMeshElementOnShape(face, geomFaceID); } } g = near; @@ -497,63 +550,88 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } QuadDelete(quad); + delete myTool; myTool = 0; + bool isOk = true; return isOk; } + //============================================================================= /*! * */ //============================================================================= -FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute - (SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) throw(SALOME_Exception) +FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMesh, + const TopoDS_Shape & aShape) + throw(SALOME_Exception) { Unexpect aCatch(SalomeException); - //SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape); - const TopoDS_Face & F = TopoDS::Face(aShape); // verify 1 wire only, with 4 edges - if (NumberOfWires(F) != 1) - { + 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; + 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()) - { + 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) - { + if (nbEdges < 4) { quad->edge[nbEdges] = E; - quad->nbPts[nbEdges] = nb + 2; // internal points + 2 extrema + if(!_quadraticMesh) { + quad->nbPts[nbEdges] = nb + 2; // internal points + 2 extrema + } + else { + int tmp = nb/2; + quad->nbPts[nbEdges] = tmp + 2; // internal not medium points + 2 extrema + } } nbEdges++; } - if (nbEdges != 4) - { + if (nbEdges != 4) { INFOS("face must have 4 edges /quadrangles"); QuadDelete(quad); return 0; } - // set normalized grid on unit square in parametric domain + return quad; +} - SetNormalizedGrid(aMesh, F, quad); +//============================================================================= +/*! + * CheckAnd2Dcompute + */ +//============================================================================= + +FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute + (SMESH_Mesh & aMesh, + const TopoDS_Shape & aShape, + const bool CreateQuadratic) throw(SALOME_Exception) +{ + Unexpect aCatch(SalomeException); + + _quadraticMesh = CreateQuadratic; + + FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape); + + if(!quad) return 0; + + // set normalized grid on unit square in parametric domain + SetNormalizedGrid(aMesh, aShape, quad); return quad; } @@ -802,6 +880,685 @@ void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, } } + +//======================================================================= +//function : ShiftQuad +//purpose : auxilary function for ComputeQuadPref +//======================================================================= +static void ShiftQuad(FaceQuadStruct* quad, const int num, bool WisF) +{ + if(num>3) return; + int i; + for(i=1; i<=num; i++) { + int nbPts3 = quad->nbPts[0]; + quad->nbPts[0] = quad->nbPts[1]; + quad->nbPts[1] = quad->nbPts[2]; + quad->nbPts[2] = quad->nbPts[3]; + quad->nbPts[3] = nbPts3; + TopoDS_Edge edge3 = quad->edge[0]; + quad->edge[0] = quad->edge[1]; + quad->edge[1] = quad->edge[2]; + quad->edge[2] = quad->edge[3]; + quad->edge[3] = edge3; + double first3 = quad->first[0]; + quad->first[0] = quad->first[1]; + quad->first[1] = quad->first[2]; + quad->first[2] = quad->first[3]; + quad->first[3] = first3; + double last3 = quad->last[0]; + quad->last[0] = quad->last[1]; + quad->last[1] = quad->last[2]; + quad->last[2] = quad->last[3]; + quad->last[3] = last3; + bool isEdgeForward3 = quad->isEdgeForward[0]; + quad->isEdgeForward[0] = quad->isEdgeForward[1]; + quad->isEdgeForward[1] = quad->isEdgeForward[2]; + quad->isEdgeForward[2] = quad->isEdgeForward[3]; + quad->isEdgeForward[3] = isEdgeForward3; + bool isEdgeOut3 = quad->isEdgeOut[0]; + quad->isEdgeOut[0] = quad->isEdgeOut[1]; + quad->isEdgeOut[1] = quad->isEdgeOut[2]; + quad->isEdgeOut[2] = quad->isEdgeOut[3]; + quad->isEdgeOut[3] = isEdgeOut3; + UVPtStruct* uv_edges3 = quad->uv_edges[0]; + quad->uv_edges[0] = quad->uv_edges[1]; + quad->uv_edges[1] = quad->uv_edges[2]; + quad->uv_edges[2] = quad->uv_edges[3]; + quad->uv_edges[3] = uv_edges3; + } + if(!WisF) { + // replacement left and right edges + int nbPts3 = quad->nbPts[1]; + quad->nbPts[1] = quad->nbPts[3]; + quad->nbPts[3] = nbPts3; + TopoDS_Edge edge3 = quad->edge[1]; + quad->edge[1] = quad->edge[3]; + quad->edge[3] = edge3; + double first3 = quad->first[1]; + quad->first[1] = quad->first[3]; + quad->first[3] = first3; + double last3 = quad->last[1]; + quad->last[1] = quad->last[2]; + quad->last[3] = last3; + bool isEdgeForward3 = quad->isEdgeForward[1]; + quad->isEdgeForward[1] = quad->isEdgeForward[3]; + quad->isEdgeForward[3] = isEdgeForward3; + bool isEdgeOut3 = quad->isEdgeOut[1]; + quad->isEdgeOut[1] = quad->isEdgeOut[3]; + quad->isEdgeOut[3] = isEdgeOut3; + UVPtStruct* uv_edges3 = quad->uv_edges[1]; + quad->uv_edges[1] = quad->uv_edges[3]; + quad->uv_edges[3] = uv_edges3; + } +} + + +//======================================================================= +//function : CalcUV +//purpose : auxilary function for ComputeQuadPref +//======================================================================= +static gp_XY CalcUV(double x0, double x1, double y0, double y1, + FaceQuadStruct* quad, + const gp_Pnt2d& a0, const gp_Pnt2d& a1, + const gp_Pnt2d& a2, const gp_Pnt2d& a3, + const Handle(Geom2d_Curve)& c2db, + const Handle(Geom2d_Curve)& c2dr, + const Handle(Geom2d_Curve)& c2dt, + const Handle(Geom2d_Curve)& c2dl) +{ + int nb = quad->nbPts[0]; + int nr = quad->nbPts[1]; + int nt = quad->nbPts[2]; + int nl = quad->nbPts[3]; + + UVPtStruct* uv_eb = quad->uv_edges[0]; + UVPtStruct* uv_er = quad->uv_edges[1]; + UVPtStruct* uv_et = quad->uv_edges[2]; + UVPtStruct* uv_el = quad->uv_edges[3]; + + double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); + double y = y0 + x * (y1 - y0); + + double param_b = uv_eb[0].param + x * (uv_eb[nb-1].param - uv_eb[0].param); + double param_t = uv_et[0].param + x * (uv_et[nt-1].param - uv_et[0].param); + double param_r = uv_er[0].param + y * (uv_er[nr-1].param - uv_er[0].param); + double param_l = uv_el[0].param + y * (uv_el[nl-1].param - uv_el[0].param); + + gp_Pnt2d p0 = c2db->Value(param_b); + gp_Pnt2d p1 = c2dr->Value(param_r); + gp_Pnt2d p2 = c2dt->Value(param_t); + gp_Pnt2d p3 = c2dl->Value(param_l); + + 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(); + + //cout<<"x0="<ShapeToIndex( F ); + + int nb = quad->nbPts[0]; + int nr = quad->nbPts[1]; + int nt = quad->nbPts[2]; + int nl = quad->nbPts[3]; + int dh = abs(nb-nt); + int dv = abs(nr-nl); + + if( dh>=dv ) { + if( nt>nb ) { + // it is a base case => not shift quad but me be replacement is need + ShiftQuad(quad,0,WisF); + } + else { + // we have to shift quad on 2 + ShiftQuad(quad,2,WisF); + } + } + else { + if( nr>nl ) { + // we have to shift quad on 3 + ShiftQuad(quad,3,WisF); + } + else { + // we have to shift quad on 1 + ShiftQuad(quad,1,WisF); + } + } + + nb = quad->nbPts[0]; + nr = quad->nbPts[1]; + nt = quad->nbPts[2]; + nl = quad->nbPts[3]; + dh = abs(nb-nt); + dv = abs(nr-nl); + int nbh = Max(nb,nt); + int nbv = Max(nr,nl); + int addh = 0; + int addv = 0; + + // orientation of face and 3 main domain for future faces + // 0 top 1 + // 1------------1 + // | | | | + // | | | | + // | L | | R | + // left | | | | rigth + // | / \ | + // | / C \ | + // |/ \| + // 0------------0 + // 0 bottom 1 + + if(dh>dv) { + addv = (dh-dv)/2; + nbv = nbv + addv; + } + else { // dv>=dh + addh = (dv-dh)/2; + nbh = nbh + addh; + } + + Handle(Geom2d_Curve) c2d[4]; + for(i=0; i<4; i++) { + c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i], F, + quad->first[i], quad->last[i]); + } + + bool loadOk = true; + for(i=0; i<2; i++) { + quad->uv_edges[i] = LoadEdgePoints2(aMesh, F, quad->edge[i], false); + if(!quad->uv_edges[i]) loadOk = false; + } + for(i=2; i<4; i++) { + quad->uv_edges[i] = LoadEdgePoints2(aMesh, F, quad->edge[i], true); + if (!quad->uv_edges[i]) loadOk = false; + } + if (!loadOk) { + INFOS("StdMeshers_Quadrangle_2D::ComputeQuadPref - LoadEdgePoints failed"); + QuadDelete( quad ); + quad = 0; + return false; + } + + UVPtStruct* uv_eb = quad->uv_edges[0]; + UVPtStruct* uv_er = quad->uv_edges[1]; + UVPtStruct* uv_et = quad->uv_edges[2]; + UVPtStruct* uv_el = quad->uv_edges[3]; + + // arrays for normalized params + //cout<<"Dump B:"<X()<<","<Y()<<","<Z()<<")"<D0(uv_eb[0].param,a[0]); + c2d[0]->D0(uv_eb[nb-1].param,a[1]); + c2d[2]->D0(uv_et[nt-1].param,a[2]); + c2d[2]->D0(uv_et[0].param,a[3]); + //cout<<" a[0]("<0) { + // add top nodes + for(i=1; i<=dl; i++) + NodesL.SetValue(i+1,nl,uv_et[i].node); + // create and add needed nodes + TColgp_SequenceOfXY UVtmp; + for(i=1; i<=dl; i++) { + double x0 = npt.Value(i+1); + double x1 = x0; + // diagonal node + double y0 = npl.Value(i+1); + double y1 = npr.Value(i+1); + gp_XY UV = CalcUV(x0, x1, y0, y1, quad, a[0], a[1], a[2], a[3], + c2d[0], c2d[1], c2d[2], c2d[3]); + gp_Pnt P = S->Value(UV.X(),UV.Y()); + SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y()); + NodesL.SetValue(i+1,1,N); + if(UVL.Length()Value(UV.X(),UV.Y()); + SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y()); + NodesL.SetValue(i+1,j,N); + if( i==dl ) UVtmp.Append(UV); + } + } + for(i=1; i<=UVtmp.Length() && UVL.Length()X()<<","<Y()<<","<Z()<<")"; + // } + // cout<AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j), + NodesL.Value(i+1,j+1), NodesL.Value(i,j+1)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + else { + SMDS_MeshFace* F = + myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1), + NodesL.Value(i+1,j+1), NodesL.Value(i+1,j)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } + } + } + else { + // fill UVL using c2d + for(i=1; iD0(uv_el[i].param,p2d); + UVL.Append(p2d.XY()); + } + } + + // step2: create faces for right domain + StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr); + // add right nodes + for(j=1; j<=nr; j++) + NodesR.SetValue(1,j,uv_er[nr-j].node); + if(dr>0) { + // add top nodes + for(i=1; i<=dr; i++) + NodesR.SetValue(i+1,1,uv_et[nt-1-i].node); + // create and add needed nodes + TColgp_SequenceOfXY UVtmp; + for(i=1; i<=dr; i++) { + double x0 = npt.Value(nt-i); + double x1 = x0; + // diagonal node + double y0 = npl.Value(i+1); + double y1 = npr.Value(i+1); + gp_XY UV = CalcUV(x0, x1, y0, y1, quad, a[0], a[1], a[2], a[3], + c2d[0], c2d[1], c2d[2], c2d[3]); + gp_Pnt P = S->Value(UV.X(),UV.Y()); + SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y()); + NodesR.SetValue(i+1,nr,N); + if(UVR.Length()Value(UV.X(),UV.Y()); + SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y()); + NodesR.SetValue(i+1,j,N); + if( i==dr ) UVtmp.Prepend(UV); + } + } + for(i=1; i<=UVtmp.Length() && UVR.Length()AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j), + NodesR.Value(i+1,j+1), NodesR.Value(i,j+1)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + else { + SMDS_MeshFace* F = + myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1), + NodesR.Value(i+1,j+1), NodesR.Value(i+1,j)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } + } + } + else { + // fill UVR using c2d + for(i=1; iD0(uv_er[i].param,p2d); + UVR.Append(p2d.XY()); + } + } + + // step3: create faces for central domain + StdMeshers_Array2OfNode NodesC(1,nb,1,nbv); + // add first string using NodesL + for(i=1; i<=dl+1; i++) + NodesC.SetValue(1,i,NodesL(i,1)); + for(i=2; i<=nl; i++) + NodesC.SetValue(1,dl+i,NodesL(dl+1,i)); + // add last string using NodesR + for(i=1; i<=dr+1; i++) + NodesC.SetValue(nb,i,NodesR(i,nr)); + for(i=1; iD0(uv_eb[i-1].param,p2d); + } + // create and add needed nodes + // add linear layers + for(i=2; iValue(UV.X(),UV.Y()); + SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y()); + NodesC.SetValue(i,nbv-nnn+j,N); + } + } + // add diagonal layers + //cout<<"UVL.Length()="<Value(u,v); + SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(N, geomFaceID, u, v); + NodesC.SetValue(j,i+1,N); + } + } + // create faces + for(i=1; iAddFace(NodesC.Value(i,j), NodesC.Value(i+1,j), + NodesC.Value(i+1,j+1), NodesC.Value(i,j+1)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + else { + SMDS_MeshFace* F = + myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1), + NodesC.Value(i+1,j+1), NodesC.Value(i+1,j)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } + } + + QuadDelete(quad); + bool isOk = true; + return isOk; +} + + +//============================================================================= +/*! + * LoadEdgePoints2 + */ +//============================================================================= +UVPtStruct* StdMeshers_Quadrangle_2D::LoadEdgePoints2 (SMESH_Mesh & aMesh, + const TopoDS_Face& F, + const TopoDS_Edge& E, + bool IsReverse) +{ + //MESSAGE("StdMeshers_Quadrangle_2D::LoadEdgePoints"); + // --- 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(); + int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes(); + + if(!_quadraticMesh) { + 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; + } + } + else { + vector nodes(nbPoints+2); + nodes[0] = idFirst; + nodes[nbPoints+1] = idLast; + nbPoints = nbPoints/2; + int nn = 1; + while(ite->more()) { + const SMDS_MeshNode* node = ite->next(); + nodes[nn++] = node; + // check if node is medium + bool IsMedium = false; + SMDS_ElemIteratorPtr itn = node->GetInverseElementIterator(); + while (itn->more()) { + const SMDS_MeshElement* elem = itn->next(); + if ( elem->GetType() != SMDSAbs_Edge ) + continue; + if(elem->IsMediumNode(node)) { + IsMedium = true; + break; + } + } + if(IsMedium) + continue; + const SMDS_EdgePosition* epos = + static_cast(node->GetPosition().get()); + double param = epos->GetUParameter(); + params[param] = node; + } + } + + 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); + + const TopoDS_Wire& W = BRepTools::OuterWire(F); + bool FisF = (F.Orientation()==TopAbs_FORWARD); + bool WisF = (W.Orientation()==TopAbs_FORWARD); + bool isForward = (E.Orientation()==TopAbs_FORWARD); + //if(isForward) cout<<"E is FORWARD"<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()->GetElements(); +// while(iter->more()) { +// const SMDS_MeshElement* elem = iter->next(); +// SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator(); +// const SMDS_MeshNode* n1 = static_cast( nodeIt->next() ); +// const SMDS_MeshNode* n2 = static_cast( nodeIt->next() ); +// const SMDS_MeshNode* n3 = static_cast( nodeIt->next() ); +// NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 )); +// myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n3)); +// myNLinkNodeMap[link] = n3; +// } +// } + map params; SMDS_NodeIteratorPtr ite = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes(); + int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes(); - 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; + if(!_quadraticMesh) { + 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; + } + } + else { + nbPoints = nbPoints/2; + while(ite->more()) { + const SMDS_MeshNode* node = ite->next(); + // check if node is medium + bool IsMedium = false; + SMDS_ElemIteratorPtr itn = node->GetInverseElementIterator(); + while (itn->more()) { + const SMDS_MeshElement* elem = itn->next(); + if ( elem->GetType() != SMDSAbs_Edge ) + continue; + if(elem->IsMediumNode(node)) { + IsMedium = true; + break; + } + } + if(IsMedium) + continue; + 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()) - { + if (nbPoints != params.size()) { MESSAGE( "BAD NODE ON EDGE POSITIONS" ); return 0; }