X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=bd55a02c400d31af28369303aacf61f07aaaf58f;hp=d806c877863982bccc64551834bbcf3372781816;hb=55510e5ccd337f72bb1caa755b131cf0916d91b3;hpb=58187d6f689f07303e5a147e853155e8c5f058d8 diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index d806c8778..bd55a02c4 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -1,60 +1,74 @@ -// SMESH SMESH : implementaion of SMESH idl descriptions +// Copyright (C) 2007-2008 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 +// +// 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. // -// Copyright (C) 2003 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 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 -// -// See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org +// 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 // +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // +// SMESH SMESH : implementaion of SMESH idl descriptions // File : StdMeshers_Quadrangle_2D.cxx // Moved here from SMESH_Quadrangle_2D.cxx // Author : Paul RASCLE, EDF // Module : SMESH -// $Header$ - -using namespace std; +// #include "StdMeshers_Quadrangle_2D.hxx" + +#include "StdMeshers_FaceSide.hxx" + +#include "StdMeshers_QuadrangleParams.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 -#include #include - +#include #include -#include -#include -#include -#include - +#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; +DEFINE_BASECOLLECTION (StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr) +DEFINE_ARRAY2(StdMeshers_Array2OfNode, + StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr) +#endif + +using namespace std; + +typedef gp_XY gp_UV; +typedef SMESH_Comment TComm; //============================================================================= /*! @@ -62,13 +76,17 @@ using namespace std; */ //============================================================================= -StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMESH_Gen* 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); + _compatibleHypothesis.push_back("QuadrangleParams"); + _compatibleHypothesis.push_back("QuadranglePreference"); + _compatibleHypothesis.push_back("TrianglePreference"); + myTool = 0; } //============================================================================= @@ -89,16 +107,75 @@ 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 + const list &hyps = + GetUsedHypothesis(aMesh, aShape, false); + const SMESHDS_Hypothesis *theHyp = 0; + + if( hyps.size() == 1 ) { + myTriaVertexID = -1; + theHyp = hyps.front(); + if(strcmp("QuadrangleParams", theHyp->GetName()) == 0) { + const StdMeshers_QuadrangleParams* theHyp1 = + (const StdMeshers_QuadrangleParams*)theHyp; + myTriaVertexID = theHyp1->GetTriaVertex(); + myQuadranglePreference= false; + myTrianglePreference= false; + } + if(strcmp("QuadranglePreference", theHyp->GetName()) == 0) { + myQuadranglePreference= true; + myTrianglePreference= false; + myTriaVertexID = -1; + } + else if(strcmp("TrianglePreference", theHyp->GetName()) == 0){ + myQuadranglePreference= false; + myTrianglePreference= true; + myTriaVertexID = -1; + } + } + + else if( hyps.size() > 1 ) { + theHyp = hyps.front(); + if(strcmp("QuadrangleParams", theHyp->GetName()) == 0) { + const StdMeshers_QuadrangleParams* theHyp1 = + (const StdMeshers_QuadrangleParams*)theHyp; + myTriaVertexID = theHyp1->GetTriaVertex(); + theHyp = hyps.back(); + if(strcmp("QuadranglePreference", theHyp->GetName()) == 0) { + myQuadranglePreference= true; + myTrianglePreference= false; + } + else if(strcmp("TrianglePreference", theHyp->GetName()) == 0){ + myQuadranglePreference= false; + myTrianglePreference= true; + } + } + else { + if(strcmp("QuadranglePreference", theHyp->GetName()) == 0) { + myQuadranglePreference= true; + myTrianglePreference= false; + } + else if(strcmp("TrianglePreference", theHyp->GetName()) == 0){ + myQuadranglePreference= false; + myTrianglePreference= true; + } + const StdMeshers_QuadrangleParams* theHyp2 = + (const StdMeshers_QuadrangleParams*)hyps.back(); + myTriaVertexID = theHyp2->GetTriaVertex(); + } + } + + else { + myQuadranglePreference = false; + myTrianglePreference = false; + myTriaVertexID = -1; + } return isOk; } @@ -110,46 +187,60 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis //============================================================================= bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape) throw (SALOME_Exception) + const TopoDS_Shape& aShape)// throw (SALOME_Exception) { - Unexpect aCatch(SalomeException); - //MESSAGE("StdMeshers_Quadrangle_2D::Compute"); + // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside + //Unexpect aCatchSalomeException); + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); aMesh.GetSubMesh(aShape); - FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape); + SMESH_MesherHelper helper(aMesh); + myTool = &helper; + + _quadraticMesh = myTool->IsQuadraticSubMesh(aShape); + + FaceQuadStruct *quad = CheckNbEdges( aMesh, aShape ); + std::auto_ptr quadDeleter( quad ); // to delete quad at exit from Compute() if (!quad) return false; + if(myQuadranglePreference) { + int n1 = quad->side[0]->NbPoints(); + int n2 = quad->side[1]->NbPoints(); + int n3 = quad->side[2]->NbPoints(); + int n4 = quad->side[3]->NbPoints(); + 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); + return ok; + } + } + + // set normalized grid on unit square in parametric domain + + if (!SetNormalizedGrid(aMesh, aShape, 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 nbdown = quad->side[0]->NbPoints(); + int nbup = quad->side[2]->NbPoints(); - int nbright = quad->nbPts[1]; - int nbleft = quad->nbPts[3]; -// bool isRightOut = (nbright > nbleft); -// bool isLeftOut = (nbright < nbleft); + int nbright = quad->side[1]->NbPoints(); + int nbleft = quad->side[3]->NbPoints(); int nbhoriz = Min(nbdown, nbup); int nbvertic = Min(nbright, nbleft); - //int nbVertices = nbhoriz * nbvertic; - //int nbQuad = (nbhoriz - 1) * (nbvertic - 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); // internal mesh nodes - int i, j, faceID = meshDS->ShapeToIndex( F ); + 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; @@ -157,11 +248,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, faceID, u, v); + meshDS->SetNodeOnFace(node, geomFaceID, u, v); quad->uv_grid[ij].node = node; } } - + // mesh faces // [2] @@ -175,18 +266,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; @@ -194,20 +284,25 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, b = quad->uv_grid[j * nbhoriz + i + 1].node; c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node; d = quad->uv_grid[(j + 1) * nbhoriz + i].node; - SMDS_MeshFace * face = meshDS->AddFace(a, b, c, d); - meshDS->SetMeshElementOnShape(face, faceID); + SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + if(face) { + 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]; - UVPtStruct *uv_e3 = quad->uv_edges[3]; + const vector& uv_e0 = quad->side[0]->GetUVPtStruct(true,0 ); + const vector& uv_e1 = quad->side[1]->GetUVPtStruct(false,1); + const vector& uv_e2 = quad->side[2]->GetUVPtStruct(true,1 ); + const vector& uv_e3 = quad->side[3]->GetUVPtStruct(false,0); + + if ( uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty() ) + return error( COMPERR_BAD_INPUT_MESH ); double eps = Precision::Confusion(); // Boundary quadrangles - + if (quad->isEdgeOut[0]) { // Down edge is out // @@ -219,14 +314,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++) { @@ -234,18 +329,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 @@ -265,15 +361,23 @@ 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 = myTool->AddFace(a, b, c); + if(face) 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, faceID); + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + + if(!myTrianglePreference){ + SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + if(face) meshDS->SetMeshElementOnShape(face, geomFaceID); + } + else { + SplitQuad(meshDS, geomFaceID, a, b, c, d); + } // if node d is not at position g - make additional triangles if (near - 1 > g) { @@ -283,8 +387,8 @@ 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, faceID); + SMDS_MeshFace* face = myTool->AddFace(a, c, d); + if(face) meshDS->SetMeshElementOnShape(face, geomFaceID); } } g = near; @@ -345,15 +449,22 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = meshDS->AddFace(a, b, c); - meshDS->SetMeshElementOnShape(face, faceID); - } else { // make quadrangle + SMDS_MeshFace* face = myTool->AddFace(a, b, c); + if(face) 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, faceID); + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + if(!myTrianglePreference){ + SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + if(face) meshDS->SetMeshElementOnShape(face, geomFaceID); + } + else { + SplitQuad(meshDS, geomFaceID, a, b, c, d); + } if (near + 1 < g) { // if d not is at g - make additional triangles for (int k = near + 1; k < g; k++) { @@ -362,8 +473,8 @@ 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, faceID); + SMDS_MeshFace* face = myTool->AddFace(a, c, d); + if(face) meshDS->SetMeshElementOnShape(face, geomFaceID); } } g = near; @@ -410,15 +521,23 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = meshDS->AddFace(a, b, c); - meshDS->SetMeshElementOnShape(face, faceID); - } else { // make quadrangle + SMDS_MeshFace* face = myTool->AddFace(a, b, c); + if(face) 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, faceID); + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + + if(!myTrianglePreference){ + SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + if(face) meshDS->SetMeshElementOnShape(face, geomFaceID); + } + else { + SplitQuad(meshDS, geomFaceID, a, b, c, d); + } if (near - 1 > g) { // if d not is at g - make additional triangles for (int k = near - 1; k > g; k--) { @@ -427,8 +546,8 @@ 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, faceID); + SMDS_MeshFace* face = myTool->AddFace(a, c, d); + if(face) meshDS->SetMeshElementOnShape(face, geomFaceID); } } g = near; @@ -472,15 +591,22 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } if (near == g) { // make triangle - SMDS_MeshFace* face = meshDS->AddFace(a, b, c); - meshDS->SetMeshElementOnShape(face, faceID); - } else { // make quadrangle + SMDS_MeshFace* face = myTool->AddFace(a, b, c); + if(face) 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, faceID); + //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d); + if(!myTrianglePreference){ + SMDS_MeshFace* face = myTool->AddFace(a, b, c, d); + if(face) meshDS->SetMeshElementOnShape(face, geomFaceID); + } + else { + SplitQuad(meshDS, geomFaceID, a, b, c, d); + } if (near + 1 < g) { // if d not is at g - make additional triangles for (int k = near + 1; k < g; k++) { @@ -489,8 +615,8 @@ 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, faceID); + SMDS_MeshFace* face = myTool->AddFace(a, c, d); + if(face) meshDS->SetMeshElementOnShape(face, geomFaceID); } } g = near; @@ -499,89 +625,484 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, } } - QuadDelete(quad); bool isOk = true; return isOk; } + //============================================================================= /*! - * + * Evaluate */ //============================================================================= -FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute - (SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) throw(SALOME_Exception) +bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + MapShapeNbElems& aResMap) + { - Unexpect aCatch(SalomeException); -// MESSAGE("StdMeshers_Quadrangle_2D::CheckAnd2Dcompute"); + aMesh.GetSubMesh(aShape); + + std::vector aNbNodes(4); + bool IsQuadratic = false; + if( !CheckNbEdgesForEvaluate( aMesh, aShape, aResMap, aNbNodes, IsQuadratic ) ) { + std::vector aResVec(17); + for(int i=0; i<17; i++) aResVec[i] = 0; + SMESH_subMesh * sm = aMesh.GetSubMesh(aShape); + aResMap.insert(std::make_pair(sm,aResVec)); + SMESH_ComputeErrorPtr& smError = sm->GetComputeError(); + smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this)); + return false; + } + + if(myQuadranglePreference) { + int n1 = aNbNodes[0]; + int n2 = aNbNodes[1]; + int n3 = aNbNodes[2]; + int n4 = aNbNodes[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 + return EvaluateQuadPref(aMesh, aShape, aNbNodes, aResMap, IsQuadratic); + //return true; + } + } + + int nbdown = aNbNodes[0]; + int nbup = aNbNodes[2]; - //SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape); + int nbright = aNbNodes[1]; + int nbleft = aNbNodes[3]; + + int nbhoriz = Min(nbdown, nbup); + int nbvertic = Min(nbright, nbleft); + + int dh = Max(nbdown, nbup) - nbhoriz; + int dv = Max(nbright, nbleft) - nbvertic; + + //int kdh = 0; + //if(dh>0) kdh = 1; + //int kdv = 0; + //if(dv>0) kdv = 1; + + int nbNodes = (nbhoriz-2)*(nbvertic-2); + //int nbFaces3 = dh + dv + kdh*(nbvertic-1)*2 + kdv*(nbhoriz-1)*2; + int nbFaces3 = dh + dv; + //if( kdh==1 && kdv==1 ) nbFaces3 -= 2; + //if( dh>0 && dv>0 ) nbFaces3 -= 2; + //int nbFaces4 = (nbhoriz-1-kdh)*(nbvertic-1-kdv); + int nbFaces4 = (nbhoriz-1)*(nbvertic-1); + + std::vector aVec(17); + for(int i=0; i<17; i++) aVec[i] = 0; + if(IsQuadratic) { + aVec[4] = nbFaces3; + aVec[6] = nbFaces4; + int nbbndedges = nbdown + nbup + nbright + nbleft -4; + int nbintedges = ( nbFaces4*4 + nbFaces3*3 - nbbndedges ) / 2; + aVec[0] = nbNodes + nbintedges; + if( aNbNodes.size()==5 ) { + aVec[4] = nbFaces3 + aNbNodes[3] -1; + aVec[6] = nbFaces4 - aNbNodes[3] +1; + } + } + else { + aVec[0] = nbNodes; + aVec[3] = nbFaces3; + aVec[5] = nbFaces4; + if( aNbNodes.size()==5 ) { + aVec[3] = nbFaces3 + aNbNodes[3] - 1; + aVec[5] = nbFaces4 - aNbNodes[3] + 1; + } + } + SMESH_subMesh * sm = aMesh.GetSubMesh(aShape); + aResMap.insert(std::make_pair(sm,aVec)); + + return true; +} + + +//================================================================================ +/*! + * \brief Return true if only two given edges meat at their common vertex + */ +//================================================================================ - // const TopoDS_Face& FF = TopoDS::Face(aShape); - // bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD); - // TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD)); +static bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1, + const TopoDS_Edge& e2, + SMESH_Mesh & mesh) +{ + TopoDS_Vertex v; + if ( !TopExp::CommonVertex( e1, e2, v )) + return false; + TopTools_ListIteratorOfListOfShape ancestIt( mesh.GetAncestors( v )); + for ( ; ancestIt.More() ; ancestIt.Next() ) + if ( ancestIt.Value().ShapeType() == TopAbs_EDGE ) + if ( !e1.IsSame( ancestIt.Value() ) && !e2.IsSame( ancestIt.Value() )) + return false; + return true; +} + +//============================================================================= +/*! + * + */ +//============================================================================= + +FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMesh, + const TopoDS_Shape & aShape) + //throw(SALOME_Exception) +{ const TopoDS_Face & F = TopoDS::Face(aShape); - //bool faceIsForward = (F.Orientation() == TopAbs_FORWARD); + const bool ignoreMediumNodes = _quadraticMesh; // verify 1 wire only, with 4 edges - - if (NumberOfWires(F) != 1) - { - MESSAGE("only 1 wire by face (quadrangles)"); + TopoDS_Vertex V; + list< TopoDS_Edge > edges; + list< int > nbEdgesInWire; + int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire); + if (nbWire != 1) { + error(COMPERR_BAD_SHAPE, TComm("Wrong number of wires: ") << nbWire); 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; + FaceQuadStruct* quad = new FaceQuadStruct; 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 + quad->side.reserve(nbEdgesInWire.front()); + + int nbSides = 0; + list< TopoDS_Edge >::iterator edgeIt = edges.begin(); + if ( nbEdgesInWire.front() == 3 ) { // exactly 3 edges + if(myTriaVertexID>0) { + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); + TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID)); + if(!V.IsNull()) { + TopoDS_Edge E1,E2,E3; + for(; edgeIt != edges.end(); ++edgeIt) { + TopoDS_Edge E = TopoDS::Edge(*edgeIt); + TopoDS_Vertex VF, VL; + TopExp::Vertices(E, VF, VL, true); + if( VF.IsSame(V) ) + E1 = E; + else if( VL.IsSame(V) ) + E3 = E; + else + E2 = E; + } + quad->side.reserve(4); + 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)); + std::vector UVPSleft = quad->side[0]->GetUVPtStruct(true,0); + std::vector UVPStop = quad->side[1]->GetUVPtStruct(false,1); + std::vector UVPSright = quad->side[2]->GetUVPtStruct(true,1); + const SMDS_MeshNode* aNode = UVPSleft[0].node; + gp_Pnt2d aPnt2d( UVPSleft[0].u, UVPSleft[0].v ); + StdMeshers_FaceSide* VertFS = + new StdMeshers_FaceSide(aNode, aPnt2d, quad->side[1]); + quad->side.push_back(VertFS); + return quad; + } } - nbEdges++; + return 0; + } + else if ( nbEdgesInWire.front() == 4 ) { // exactly 4 edges + for ( ; edgeIt != edges.end(); ++edgeIt, nbSides++ ) + quad->side.push_back( new StdMeshers_FaceSide(F, *edgeIt, &aMesh, + nbSides 4 ) { // more than 4 edges - try to unite some + list< TopoDS_Edge > sideEdges; + while ( !edges.empty()) { + sideEdges.clear(); + sideEdges.splice( sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end() + bool sameSide = true; + while ( !edges.empty() && sameSide ) { + sameSide = SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() ); + if ( sameSide ) + sideEdges.splice( sideEdges.end(), edges, edges.begin()); + } + if ( nbSides == 0 ) { // go backward from the first edge + sameSide = true; + while ( !edges.empty() && sameSide ) { + sameSide = SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() ); + if ( sameSide ) + sideEdges.splice( sideEdges.begin(), edges, --edges.end()); + } + } + quad->side.push_back( new StdMeshers_FaceSide(F, sideEdges, &aMesh, + nbSidesside.clear(); + quad->side.reserve(nbEdgesInWire.front()); + nbSides = 0; + + SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire); + while ( !edges.empty()) { + sideEdges.clear(); + sideEdges.splice( sideEdges.end(), edges, edges.begin()); + bool sameSide = true; + while ( !edges.empty() && sameSide ) { + sameSide = + SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() ) && + twoEdgesMeatAtVertex( sideEdges.back(), edges.front(), aMesh ); + if ( sameSide ) + sideEdges.splice( sideEdges.end(), edges, edges.begin()); + } + if ( nbSides == 0 ) { // go backward from the first edge + sameSide = true; + while ( !edges.empty() && sameSide ) { + sameSide = + SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() ) && + twoEdgesMeatAtVertex( sideEdges.front(), edges.back(), aMesh ); + if ( sameSide ) + sideEdges.splice( sideEdges.begin(), edges, --edges.end()); + } + } + quad->side.push_back( new StdMeshers_FaceSide(F, sideEdges, &aMesh, + nbSidesside[i]->NbEdges(); ++e ) + MESSAGE ( myTool->GetMeshDS()->ShapeToIndex( quad->side[i]->Edge( e )) << " " ); + MESSAGE ( ")\n" ); + } + //cout << endl; +#endif + if ( !nbSides ) + nbSides = nbEdgesInWire.front(); + error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides); + delete quad; + quad = 0; } - if (nbEdges != 4) - { - MESSAGE("face must have 4 edges /quadrangles"); - QuadDelete(quad); - return 0; - //throw SALOME_Exception(LOCALIZED("face must have 4 edges /quadrangles")); + return quad; +} + + +//============================================================================= +/*! + * + */ +//============================================================================= + +bool StdMeshers_Quadrangle_2D::CheckNbEdgesForEvaluate(SMESH_Mesh& aMesh, + const TopoDS_Shape & aShape, + MapShapeNbElems& aResMap, + std::vector& aNbNodes, + bool& IsQuadratic) + +{ + const TopoDS_Face & F = TopoDS::Face(aShape); + + // verify 1 wire only, with 4 edges + TopoDS_Vertex V; + list< TopoDS_Edge > edges; + list< int > nbEdgesInWire; + int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire); + if (nbWire != 1) { + return false; } -// 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")); -// } + aNbNodes.resize(4); - // set normalized grid on unit square in parametric domain + int nbSides = 0; + list< TopoDS_Edge >::iterator edgeIt = edges.begin(); + SMESH_subMesh * sm = aMesh.GetSubMesh( *edgeIt ); + MapShapeNbElemsItr anIt = aResMap.find(sm); + if(anIt==aResMap.end()) { + return false; + } + std::vector aVec = (*anIt).second; + IsQuadratic = (aVec[2] > aVec[1]); + if ( nbEdgesInWire.front() == 3 ) { // exactly 3 edges + if(myTriaVertexID>0) { + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); + TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID)); + if(!V.IsNull()) { + TopoDS_Edge E1,E2,E3; + for(; edgeIt != edges.end(); ++edgeIt) { + TopoDS_Edge E = TopoDS::Edge(*edgeIt); + TopoDS_Vertex VF, VL; + TopExp::Vertices(E, VF, VL, true); + if( VF.IsSame(V) ) + E1 = E; + else if( VL.IsSame(V) ) + E3 = E; + else + E2 = E; + } + SMESH_subMesh * sm = aMesh.GetSubMesh(E1); + MapShapeNbElemsItr anIt = aResMap.find(sm); + if(anIt==aResMap.end()) return false; + std::vector aVec = (*anIt).second; + if(IsQuadratic) + aNbNodes[0] = (aVec[0]-1)/2 + 2; + else + aNbNodes[0] = aVec[0] + 2; + sm = aMesh.GetSubMesh(E2); + anIt = aResMap.find(sm); + if(anIt==aResMap.end()) return false; + aVec = (*anIt).second; + if(IsQuadratic) + aNbNodes[1] = (aVec[0]-1)/2 + 2; + else + aNbNodes[1] = aVec[0] + 2; + sm = aMesh.GetSubMesh(E3); + anIt = aResMap.find(sm); + if(anIt==aResMap.end()) return false; + aVec = (*anIt).second; + if(IsQuadratic) + aNbNodes[2] = (aVec[0]-1)/2 + 2; + else + aNbNodes[2] = aVec[0] + 2; + aNbNodes[3] = aNbNodes[1]; + aNbNodes.resize(5); + nbSides = 4; + } + } + } + if ( nbEdgesInWire.front() == 4 ) { // exactly 4 edges + for(; edgeIt != edges.end(); edgeIt++) { + SMESH_subMesh * sm = aMesh.GetSubMesh( *edgeIt ); + MapShapeNbElemsItr anIt = aResMap.find(sm); + if(anIt==aResMap.end()) { + return false; + } + std::vector aVec = (*anIt).second; + if(IsQuadratic) + aNbNodes[nbSides] = (aVec[0]-1)/2 + 2; + else + aNbNodes[nbSides] = aVec[0] + 2; + nbSides++; + } + } + else if ( nbEdgesInWire.front() > 4 ) { // more than 4 edges - try to unite some + list< TopoDS_Edge > sideEdges; + while ( !edges.empty()) { + sideEdges.clear(); + sideEdges.splice( sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end() + bool sameSide = true; + while ( !edges.empty() && sameSide ) { + sameSide = SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() ); + if ( sameSide ) + sideEdges.splice( sideEdges.end(), edges, edges.begin()); + } + if ( nbSides == 0 ) { // go backward from the first edge + sameSide = true; + while ( !edges.empty() && sameSide ) { + sameSide = SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() ); + if ( sameSide ) + sideEdges.splice( sideEdges.begin(), edges, --edges.end()); + } + } + list::iterator ite = sideEdges.begin(); + aNbNodes[nbSides] = 1; + for(; ite!=sideEdges.end(); ite++) { + SMESH_subMesh * sm = aMesh.GetSubMesh( *ite ); + MapShapeNbElemsItr anIt = aResMap.find(sm); + if(anIt==aResMap.end()) { + return false; + } + std::vector aVec = (*anIt).second; + if(IsQuadratic) + aNbNodes[nbSides] += (aVec[0]-1)/2 + 1; + else + aNbNodes[nbSides] += aVec[0] + 1; + } + ++nbSides; + } + // issue 20222. Try to unite only edges shared by two same faces + if (nbSides < 4) { + nbSides = 0; + SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire); + while ( !edges.empty()) { + sideEdges.clear(); + sideEdges.splice( sideEdges.end(), edges, edges.begin()); + bool sameSide = true; + while ( !edges.empty() && sameSide ) { + sameSide = + SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() ) && + twoEdgesMeatAtVertex( sideEdges.back(), edges.front(), aMesh ); + if ( sameSide ) + sideEdges.splice( sideEdges.end(), edges, edges.begin()); + } + if ( nbSides == 0 ) { // go backward from the first edge + sameSide = true; + while ( !edges.empty() && sameSide ) { + sameSide = + SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() ) && + twoEdgesMeatAtVertex( sideEdges.front(), edges.back(), aMesh ); + if ( sameSide ) + sideEdges.splice( sideEdges.begin(), edges, --edges.end()); + } + } + list::iterator ite = sideEdges.begin(); + aNbNodes[nbSides] = 1; + for(; ite!=sideEdges.end(); ite++) { + SMESH_subMesh * sm = aMesh.GetSubMesh( *ite ); + MapShapeNbElemsItr anIt = aResMap.find(sm); + if(anIt==aResMap.end()) { + return false; + } + std::vector aVec = (*anIt).second; + if(IsQuadratic) + aNbNodes[nbSides] += (aVec[0]-1)/2 + 1; + else + aNbNodes[nbSides] += aVec[0] + 1; + } + ++nbSides; + } + } + } + if (nbSides != 4) { + if ( !nbSides ) + nbSides = nbEdgesInWire.front(); + error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides); + return false; + } + + return true; +} + + +//============================================================================= +/*! + * CheckAnd2Dcompute + */ +//============================================================================= - SetNormalizedGrid(aMesh, F, quad); +FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute + (SMESH_Mesh & aMesh, + const TopoDS_Shape & aShape, + const bool CreateQuadratic) //throw(SALOME_Exception) +{ + _quadraticMesh = CreateQuadratic; + + FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape); + + if(!quad) return 0; + + // set normalized grid on unit square in parametric domain + bool stat = SetNormalizedGrid(aMesh, aShape, quad); + if(!stat) { + if(!quad) + delete quad; + quad = 0; + } return quad; } @@ -592,20 +1113,23 @@ FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute */ //============================================================================= -void StdMeshers_Quadrangle_2D::QuadDelete (FaceQuadStruct * quad) +faceQuadStruct::~faceQuadStruct() { - //MESSAGE("StdMeshers_Quadrangle_2D::QuadDelete"); - if (quad) + for (int i = 0; i < side.size(); i++) { + if (side[i]) delete side[i]; + } + if (uv_grid) delete [] uv_grid; +} + +namespace { + inline const vector& GetUVPtStructIn(FaceQuadStruct* quad, int i, int nbSeg) { - 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; + bool isXConst = ( i == BOTTOM_SIDE || i == TOP_SIDE ); + double constValue = ( i == BOTTOM_SIDE || i == LEFT_SIDE ) ? 0 : 1; + return + quad->isEdgeOut[i] ? + quad->side[i]->SimulateUVPtStruct(nbSeg,isXConst,constValue) : + quad->side[i]->GetUVPtStruct(isXConst,constValue); } } @@ -615,18 +1139,17 @@ void StdMeshers_Quadrangle_2D::QuadDelete (FaceQuadStruct * quad) */ //============================================================================= -void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, +bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, const TopoDS_Shape& aShape, - FaceQuadStruct* quad) throw (SALOME_Exception) + 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é // MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid"); - const TopoDS_Face& F = TopoDS::Face(aShape); +// const TopoDS_Face& F = TopoDS::Face(aShape); // 1 --- find orientation of the 4 edges, by test on extrema @@ -642,102 +1165,26 @@ void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, // =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]); - if (!quad->uv_edges[i]) loadOk = false; - // 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]); - if (!quad->uv_edges[i]) loadOk = false; - // !quad->isEdgeForward[i]); - } - - if (!loadOk) - { -// MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid - LoadEdgePoints failed"); - QuadDelete( quad ); - quad = 0; - return; - } // 3 --- 2D normalized values on unit square [0..1][0..1] -// int nbdown = quad->nbPts[0]; -// int nbright = quad->nbPts[1]; - int nbhoriz = Min(quad->nbPts[0], quad->nbPts[2]); - int nbvertic = Min(quad->nbPts[1], quad->nbPts[3]); -// MESSAGE("nbhoriz, nbvertic = " << nbhoriz << nbvertic); + int nbhoriz = Min(quad->side[0]->NbPoints(), quad->side[2]->NbPoints()); + int nbvertic = Min(quad->side[1]->NbPoints(), quad->side[3]->NbPoints()); + + quad->isEdgeOut[0] = (quad->side[0]->NbPoints() > quad->side[2]->NbPoints()); + quad->isEdgeOut[1] = (quad->side[1]->NbPoints() > quad->side[3]->NbPoints()); + quad->isEdgeOut[2] = (quad->side[2]->NbPoints() > quad->side[0]->NbPoints()); + quad->isEdgeOut[3] = (quad->side[3]->NbPoints() > quad->side[1]->NbPoints()); - 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]); + UVPtStruct *uv_grid = quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz]; - quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz]; + const vector& uv_e0 = GetUVPtStructIn( quad, 0, nbhoriz - 1 ); + const vector& uv_e1 = GetUVPtStructIn( quad, 1, nbvertic - 1 ); + const vector& uv_e2 = GetUVPtStructIn( quad, 2, nbhoriz - 1 ); + const vector& uv_e3 = GetUVPtStructIn( quad, 3, nbvertic - 1 ); - 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]; + if ( uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty() ) + //return error( "Can't find nodes on sides"); + return error( COMPERR_BAD_INPUT_MESH ); // nodes Id on "in" edges if (! quad->isEdgeOut[0]) { @@ -769,26 +1216,9 @@ void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, } } - // 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++) - { + 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 @@ -807,305 +1237,856 @@ void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, } // 4 --- projection on 2d domain (u,v) - gp_Pnt2d a0 = pf[0]; - gp_Pnt2d a1 = pf[1]; - gp_Pnt2d a2 = pf[2]; - gp_Pnt2d a3 = pf[3]; + 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++) - { - for (int j = 0; j < nbvertic; j++) - { + 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].param + x * (uv_e0[nbhoriz - 1].param - uv_e0[0].param); // sud - double param_2 = uv_e2[0].param + x * (uv_e2[nbhoriz - 1].param - uv_e2[0].param); // nord - double param_1 = uv_e1[0].param + y * (uv_e1[nbvertic - 1].param - uv_e1[0].param); // est - double param_3 = uv_e3[0].param + y * (uv_e3[nbvertic - 1].param - uv_e3[0].param); // ouest + 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 "<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); + 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(); - 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(); + gp_UV uv = (1 - y) * p0 + x * p1 + y * p2 + (1 - x) * p3; + uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3; - 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 = uv.X(); + uv_grid[ij].v = uv.Y(); + } + } + return true; +} - uv_grid[ij].u = u; - uv_grid[ij].v = v; +//======================================================================= +//function : ShiftQuad +//purpose : auxilary function for ComputeQuadPref +//======================================================================= - //MESSAGE("-uv- "<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 ); + bool newForward = ( id < TOP_SIDE ); + if ( wasForward != newForward ) + side[ i ]->Reverse(); + quad->side[ id ] = side[ i ]; } } -//============================================================================= -/*! - * LoadEdgePoints - */ -//============================================================================= -UVPtStruct* StdMeshers_Quadrangle_2D::LoadEdgePoints (SMESH_Mesh & aMesh, - const TopoDS_Face& F, - const TopoDS_Edge& E, - double first, double last) -// bool isForward) +//======================================================================= +//function : CalcUV +//purpose : auxilary function for ComputeQuadPref +//======================================================================= + +static gp_UV CalcUV(double x0, double x1, double y0, double y1, + FaceQuadStruct* quad, + const gp_UV& a0, const gp_UV& a1, + const gp_UV& a2, const gp_UV& a3) { - //MESSAGE("StdMeshers_Quadrangle_2D::LoadEdgePoints"); + const vector& uv_eb = quad->side[0]->GetUVPtStruct(true,0 ); + const vector& uv_er = quad->side[1]->GetUVPtStruct(false,1); + const vector& uv_et = quad->side[2]->GetUVPtStruct(true,1 ); + const vector& uv_el = quad->side[3]->GetUVPtStruct(false,0); - //SMDS_Mesh* meshDS = aMesh.GetMeshDS(); + double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); + double y = y0 + x * (y1 - y0); - // --- IDNodes of first and last Vertex + 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); - TopoDS_Vertex VFirst, VLast; - TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l + 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(); - 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(); + gp_UV uv = p0 * (1 - y) + p1 * x + p2 * y + p3 * (1 - x); - 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(); + uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3; - // --- edge internal IDNodes (relies on good order storage, not checked) + return uv; +} - map params; - SMDS_NodeIteratorPtr ite = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes(); +//======================================================================= +//function : CalcUV2 +//purpose : auxilary function for ComputeQuadPref +//======================================================================= - 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; - } +static gp_UV CalcUV2(double x, double y, + FaceQuadStruct* quad, + const gp_UV& a0, const gp_UV& a1, + const gp_UV& a2, const gp_UV& a3) +{ + const vector& uv_eb = quad->side[0]->GetUVPtStruct(true,0 ); + const vector& uv_er = quad->side[1]->GetUVPtStruct(false,1); + const vector& uv_et = quad->side[2]->GetUVPtStruct(true,1 ); + const vector& uv_el = quad->side[3]->GetUVPtStruct(false,0); - 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 x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); + //double y = y0 + x * (y1 - y0); - double f, l; - Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l); + 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); - 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 "<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(); - ASSERT(paramin != paramax); - for (int i = 0; i < nbPoints + 2; i++) - { - uvslf[i].normParam = (uvslf[i].param - paramin) / (paramax - paramin); - //SCRUTE(uvslf[i].normParam); - } + gp_UV uv = p0 * (1 - y) + p1 * x + p2 * y + p3 * (1 - x); + + uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3; - return uvslf; + return uv; } -//============================================================================= + +//======================================================================= /*! - * MakeEdgePoints + * Create only quandrangle faces */ -//============================================================================= -UVPtStruct* StdMeshers_Quadrangle_2D::MakeEdgePoints (SMESH_Mesh & aMesh, - const TopoDS_Face& F, - const TopoDS_Edge& E, - double first, double last, - int nb_segm) +//======================================================================= + +bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh, + const TopoDS_Shape& aShape, + FaceQuadStruct* quad) { -// MESSAGE("StdMeshers_Quadrangle_2D::MakeEdgePoints"); - - UVPtStruct* uvslf = new UVPtStruct[nb_segm + 1]; - list 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); - } + // Auxilary key in order to keep old variant + // of meshing after implementation new variant + // for bug 0016220 from Mantis. + bool OldVersion = false; + + SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + const TopoDS_Face& F = TopoDS::Face(aShape); + Handle(Geom_Surface) S = BRep_Tool::Surface(F); +// const TopoDS_Wire& W = BRepTools::OuterWire(F); + bool WisF = true; +// if(W.Orientation()==TopAbs_FORWARD) +// WisF = true; + //if(WisF) cout<<"W is FORWARD"<ShapeToIndex( F ); + + int nb = quad->side[0]->NbPoints(); + int nr = quad->side[1]->NbPoints(); + int nt = quad->side[2]->NbPoints(); + int nl = quad->side[3]->NbPoints(); + 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); } - catch (Standard_Failure) { - return 0; + else { + // we have to shift quad on 2 + ShiftQuad(quad,2,WisF); } } - 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); + else { + if( nr>nl ) { + // we have to shift quad on 1 + ShiftQuad(quad,1,WisF); + } + else { + // we have to shift quad on 3 + ShiftQuad(quad,3,WisF); } } - 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++; + nb = quad->side[0]->NbPoints(); + nr = quad->side[1]->NbPoints(); + nt = quad->side[2]->NbPoints(); + nl = quad->side[3]->NbPoints(); + dh = abs(nb-nt); + dv = abs(nr-nl); + int nbh = Max(nb,nt); + int nbv = Max(nr,nl); + int addh = 0; + int addv = 0; + + // ----------- Old version --------------- + // 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 + + // ----------- New version --------------- + // orientation of face and 3 main domain for future faces + // 0 top 1 + // 1------------1 + // | |____| | + // | / \ | + // | / C \ | + // left |/________\| rigth + // | | + // | | + // | | + // 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; + } + + const vector& uv_eb = quad->side[0]->GetUVPtStruct(true,0 ); + const vector& uv_er = quad->side[1]->GetUVPtStruct(false,1); + const vector& uv_et = quad->side[2]->GetUVPtStruct(true,1 ); + const vector& uv_el = quad->side[3]->GetUVPtStruct(false,0); + + // arrays for normalized params + //cout<<"Dump B:"<X()<<","<Y()<<","<Z()<<")"<::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++; + // insert to left + dl = nbv - nl; + dpr = (npl.Value(2) - npl.Value(1))/(dl+1); + for(i=1; i<=dl; i++) { + npl.InsertAfter(1,npl.Value(2)-dpr); + } + } + //cout<<"npb:"; + //for(i=1; i<=npb.Length(); i++) { + // cout<<" "<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_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3); + 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)); + if(F) 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)); + if(F) meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } + } + } + else { + // fill UVL using c2d + for(i=1; i0) { + // 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_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3); + 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)); + if(F) 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)); + if(F) meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } + } + } + else { + // fill UVR using c2d + for(i=1; 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)); + if(F) 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)); + if(F) meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } } } - return uvslf; -} + else { // New version (!OldVersion) + // step1: create faces for bottom rectangle domain + StdMeshers_Array2OfNode NodesBRD(1,nb,1,nnn-1); + // fill UVL and UVR using c2d + for(j=0; jValue(u,v); + SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(N, geomFaceID, u, v); + NodesBRD.SetValue(j,i+1,N); + } + } + int nbf=0; + for(j=1; jAddFace(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), + NodesBRD.Value(i+1,j+1), NodesBRD.Value(i+1,j)); + if(F) meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } + } + int drl = abs(nr-nl); + // create faces for region C + StdMeshers_Array2OfNode NodesC(1,nb,1,drl+1+addv); + // add nodes from previous region + for(j=1; j<=nb; j++) { + NodesC.SetValue(j,1,NodesBRD.Value(j,nnn-1)); + } + if( (drl+addv) > 0 ) { + int n1,n2; + if(nr>nl) { + n1 = 1; + n2 = drl + 1; + TColgp_SequenceOfXY UVtmp; + double drparam = npr.Value(nr) - npr.Value(nnn-1); + double dlparam = npl.Value(nnn) - npl.Value(nnn-1); + double y0,y1; + for(i=1; i<=drl; i++) { + // add existed nodes from right edge + NodesC.SetValue(nb,i+1,uv_er[nnn+i-2].node); + //double dtparam = npt.Value(i+1); + y1 = npr.Value(nnn+i-1); // param on right edge + double dpar = (y1 - npr.Value(nnn-1))/drparam; + y0 = npl.Value(nnn-1) + dpar*dlparam; // param on left edge + double dy = y1 - y0; + for(j=1; jValue(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(j,i+1,N); + } + } + double dy0 = (1-y0)/(addv+1); + double dy1 = (1-y1)/(addv+1); + for(i=1; i<=addv; i++) { + double yy0 = y0 + dy0*i; + double yy1 = y1 + dy1*i; + double dyy = yy1 - yy0; + for(j=1; j<=nb; j++) { + double x = npt.Value(i+1+drl) + + npb.Value(j) * ( npt.Value(nt-i) - npt.Value(i+1+drl) ); + double y = yy0 + dyy*x; + gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3); + 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()); + NodesC.SetValue(j,i+drl+1,N); + } + } + } + else { // nrValue(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(j,i+1,N); + } + } + double dy0 = (1-y0)/(addv+1); + double dy1 = (1-y1)/(addv+1); + for(i=1; i<=addv; i++) { + double yy0 = y0 + dy0*i; + double yy1 = y1 + dy1*i; + double dyy = yy1 - yy0; + for(j=1; j<=nb; j++) { + double x = npt.Value(i+1) + + npb.Value(j) * ( npt.Value(nt-i-drl) - npt.Value(i+1) ); + double y = yy0 + dyy*x; + gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3); + 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()); + NodesC.SetValue(j,i+drl+1,N); + } + } + } + // create faces + for(j=1; j<=drl+addv; j++) { + for(i=1; iAddFace(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), + NodesC.Value(i+1,j+1), NodesC.Value(i+1,j)); + if(F) meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } + } // end nr=n2; i--) { + nnn++; + NodesLast.SetValue(nnn,1,NodesC.Value(nb,i)); + } + for(i=1; iAddFace(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), + NodesLast.Value(i+1,2), NodesLast.Value(i+1,2)); + if(F) meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } + } // if( (drl+addv) > 0 ) -ostream & StdMeshers_Quadrangle_2D::SaveTo(ostream & save) -{ - return save; + } // end new version implementation + + bool isOk = true; + return isOk; } -//============================================================================= + +//======================================================================= /*! - * + * Evaluate only quandrangle faces */ -//============================================================================= +//======================================================================= -istream & StdMeshers_Quadrangle_2D::LoadFrom(istream & load) +bool StdMeshers_Quadrangle_2D::EvaluateQuadPref(SMESH_Mesh & aMesh, + const TopoDS_Shape& aShape, + std::vector& aNbNodes, + MapShapeNbElems& aResMap, + bool IsQuadratic) { - return load; -} + // Auxilary key in order to keep old variant + // of meshing after implementation new variant + // for bug 0016220 from Mantis. + bool OldVersion = false; -//============================================================================= -/*! - * - */ -//============================================================================= + const TopoDS_Face& F = TopoDS::Face(aShape); + Handle(Geom_Surface) S = BRep_Tool::Surface(F); -ostream & operator <<(ostream & save, StdMeshers_Quadrangle_2D & hyp) -{ - return hyp.SaveTo( save ); + int nb = aNbNodes[0]; + int nr = aNbNodes[1]; + int nt = aNbNodes[2]; + int nl = aNbNodes[3]; + int dh = abs(nb-nt); + int dv = abs(nr-nl); + + if( dh>=dv ) { + if( nt>nb ) { + // it is a base case => not shift + } + else { + // we have to shift on 2 + nb = aNbNodes[2]; + nr = aNbNodes[3]; + nt = aNbNodes[0]; + nl = aNbNodes[1]; + } + } + else { + if( nr>nl ) { + // we have to shift quad on 1 + nb = aNbNodes[3]; + nr = aNbNodes[0]; + nt = aNbNodes[1]; + nl = aNbNodes[2]; + } + else { + // we have to shift quad on 3 + nb = aNbNodes[1]; + nr = aNbNodes[2]; + nt = aNbNodes[3]; + nl = aNbNodes[0]; + } + } + + dh = abs(nb-nt); + dv = abs(nr-nl); + int nbh = Max(nb,nt); + int nbv = Max(nr,nl); + int addh = 0; + int addv = 0; + + if(dh>dv) { + addv = (dh-dv)/2; + nbv = nbv + addv; + } + else { // dv>=dh + addh = (dv-dh)/2; + nbh = nbh + addh; + } + + int dl,dr; + if(OldVersion) { + // add some params to right and left after the first param + // insert to right + dr = nbv - nr; + // insert to left + dl = nbv - nl; + } + + int nnn = Min(nr,nl); + + int nbNodes = 0; + int nbFaces = 0; + if(OldVersion) { + // step1: create faces for left domain + if(dl>0) { + nbNodes += dl*(nl-1); + nbFaces += dl*(nl-1); + } + // step2: create faces for right domain + if(dr>0) { + nbNodes += dr*(nr-1); + nbFaces += dr*(nr-1); + } + // step3: create faces for central domain + nbNodes += (nb-2)*(nnn-1) + (nbv-nnn-1)*(nb-2); + nbFaces += (nb-1)*(nbv-1); + } + else { // New version (!OldVersion) + nbNodes += (nnn-2)*(nb-2); + nbFaces += (nnn-2)*(nb-1); + int drl = abs(nr-nl); + nbNodes += drl*(nb-1) + addv*nb; + nbFaces += (drl+addv)*(nb-1) + (nt-1); + } // end new version implementation + + std::vector aVec(17); + for(int i=0; i<17; i++) aVec[i] = 0; + if(IsQuadratic) { + aVec[6] = nbFaces; + aVec[0] = nbNodes + nbFaces*4; + if( aNbNodes.size()==5 ) { + aVec[4] = aNbNodes[3] - 1; + aVec[6] = nbFaces - aNbNodes[3] + 1; + } + } + else { + aVec[0] = nbNodes; + aVec[5] = nbFaces; + if( aNbNodes.size()==5 ) { + aVec[3] = aNbNodes[3] - 1; + aVec[5] = nbFaces - aNbNodes[3] + 1; + } + } + SMESH_subMesh * sm = aMesh.GetSubMesh(aShape); + aResMap.insert(std::make_pair(sm,aVec)); + + return true; } + //============================================================================= -/*! - * +/*! Split quadrangle in to 2 triangles by smallest diagonal + * */ //============================================================================= - -istream & operator >>(istream & load, StdMeshers_Quadrangle_2D & hyp) +void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS, + int theFaceID, + const SMDS_MeshNode* theNode1, + const SMDS_MeshNode* theNode2, + const SMDS_MeshNode* theNode3, + const SMDS_MeshNode* theNode4) { - return hyp.LoadFrom( load ); + gp_Pnt a(theNode1->X(),theNode1->Y(),theNode1->Z()); + gp_Pnt b(theNode2->X(),theNode2->Y(),theNode2->Z()); + gp_Pnt c(theNode3->X(),theNode3->Y(),theNode3->Z()); + 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); + if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID ); + face = myTool->AddFace(theNode2, theNode3, theNode4); + if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID ); + + } + else{ + face = myTool->AddFace(theNode1, theNode2 ,theNode3); + if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID ); + face = myTool->AddFace(theNode1, theNode3, theNode4); + if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID ); + } } + +