X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Quadrangle_2D.cxx;h=0a7e3301754ee315929c1e865e4f646dca1dfd59;hb=4c65637b3ba0be986e1ce6e952689b2686475b2f;hp=1275988675e93267201dc26caaf5a346007792a8;hpb=f7aba4830d53719b963fdb7fccc98b760fdef2d1;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index 127598867..0a7e33017 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2014 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 @@ -6,7 +6,7 @@ // 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. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -42,8 +42,10 @@ #include "StdMeshers_QuadrangleParams.hxx" #include "StdMeshers_ViscousLayers2D.hxx" +#include #include #include +#include #include #include #include @@ -88,6 +90,7 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, myTrianglePreference(false), myTriaVertexID(-1), myNeedSmooth(false), + myCheckOri(false), myParams( NULL ), myQuadType(QUAD_STANDARD), myHelper( NULL ) @@ -213,15 +216,20 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, const TopoDS_Face& F = TopoDS::Face(aShape); aMesh.GetSubMesh( F ); + // do not initialize my fields before this as StdMeshers_ViscousLayers2D + // can call Compute() recursively + SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F ); + if ( !proxyMesh ) + return false; + + myProxyMesh = proxyMesh; + SMESH_MesherHelper helper (aMesh); myHelper = &helper; - myProxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F ); - if ( !myProxyMesh ) - return false; - _quadraticMesh = myHelper->IsQuadraticSubMesh(aShape); myNeedSmooth = false; + myCheckOri = false; FaceQuadStruct::Ptr quad = CheckNbEdges( aMesh, F, /*considerMesh=*/true ); if (!quad) @@ -288,6 +296,9 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, if ( res == COMPUTE_OK && myNeedSmooth ) smooth( quad ); + if ( res == COMPUTE_OK ) + res = check(); + return ( res == COMPUTE_OK ); } @@ -317,17 +328,25 @@ bool StdMeshers_Quadrangle_2D::computeTriangles(SMESH_Mesh& aMesh, if ( !setNormalizedGrid( quad )) return false; - if ( quad->nbNodeOut( QUAD_BOTTOM_SIDE )) - { - splitQuad( quad, 0, 1 ); - } if ( quad->nbNodeOut( QUAD_TOP_SIDE )) { splitQuad( quad, 0, quad->jSize-2 ); } + if ( quad->nbNodeOut( QUAD_BOTTOM_SIDE )) // this should not happen + { + splitQuad( quad, 0, 1 ); + } FaceQuadStruct::Ptr newQuad = myQuadList.back(); if ( quad != newQuad ) // split done { + { + FaceQuadStruct::Ptr botQuad = // a bottom part + ( quad->side[ QUAD_LEFT_SIDE ].from == 0 ) ? quad : newQuad; + if ( botQuad->nbNodeOut( QUAD_LEFT_SIDE ) > 0 ) + botQuad->side[ QUAD_LEFT_SIDE ].to += botQuad->nbNodeOut( QUAD_LEFT_SIDE ); + else if ( botQuad->nbNodeOut( QUAD_RIGHT_SIDE ) > 0 ) + botQuad->side[ QUAD_RIGHT_SIDE ].to += botQuad->nbNodeOut( QUAD_RIGHT_SIDE ); + } // make quad be a greatest one if ( quad->side[ QUAD_LEFT_SIDE ].NbPoints() == 2 || quad->side[ QUAD_RIGHT_SIDE ].NbPoints() == 2 ) @@ -455,7 +474,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, int nbright = (int) uv_e1.size(); int nbleft = (int) uv_e3.size(); - if (quad->nbNodeOut(0) && nbvertic == 2) + if (quad->nbNodeOut(0) && nbvertic == 2) // this should not occure { // Down edge is out // @@ -649,7 +668,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, } // right or left boundary quadrangles - if (quad->nbNodeOut( QUAD_RIGHT_SIDE ) && nbhoriz == 2) + if (quad->nbNodeOut( QUAD_RIGHT_SIDE ) && nbhoriz == 2) // this should not occure { int g = 0; // last processed node in the grid int stop = nbright - 1; @@ -725,9 +744,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh, // MESSAGE("left edge is out"); int g = nbvertic - 1; // last processed node in the grid int stop = 0; - i = nbleft - 1; - if (quad->side[3].from != stop ) stop++; - if (quad->side[3].to != i ) i--; + i = quad->side[ QUAD_LEFT_SIDE ].to-1; // nbleft - 1; for (; i > stop; i--) { const SMDS_MeshNode *a, *b, *c, *d; a = uv_e3[i].node; @@ -892,6 +909,38 @@ bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh, return true; } +//================================================================================ +/*! + * \brief Return true if the algorithm can mesh this shape + * \param [in] aShape - shape to check + * \param [in] toCheckAll - if true, this check returns OK if all shapes are OK, + * else, returns OK if at least one shape is OK + */ +//================================================================================ + +bool StdMeshers_Quadrangle_2D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll ) +{ + int nbFoundFaces = 0; + for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces ) + { + const TopoDS_Shape& aFace = exp.Current(); + int nbWire = SMESH_MesherHelper::Count( aFace, TopAbs_WIRE, false ); + if ( nbWire != 1 ) { + if ( toCheckAll ) return false; + continue; + } + + int nbNoDegenEdges = 0; + TopExp_Explorer eExp( aFace, TopAbs_EDGE ); + for ( ; eExp.More() && nbNoDegenEdges < 3; eExp.Next() ) { + if ( !SMESH_Algo::isDegenerated( TopoDS::Edge( eExp.Current() ))) + ++nbNoDegenEdges; + } + if ( toCheckAll && nbNoDegenEdges < 3 ) return false; + if ( !toCheckAll && nbNoDegenEdges >= 3 ) return true; + } + return ( toCheckAll && nbFoundFaces != 0 ); +} //================================================================================ /*! @@ -931,7 +980,7 @@ FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD ); const bool ignoreMediumNodes = _quadraticMesh; - // verify 1 wire only, with 4 edges + // verify 1 wire only list< TopoDS_Edge > edges; list< int > nbEdgesInWire; int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire); @@ -1304,38 +1353,30 @@ bool StdMeshers_Quadrangle_2D::setNormalizedGrid (FaceQuadStruct::Ptr quad) // min max 0 x0 1 // =down // + const FaceQuadStruct::Side & bSide = quad->side[0]; + const FaceQuadStruct::Side & rSide = quad->side[1]; + const FaceQuadStruct::Side & tSide = quad->side[2]; + const FaceQuadStruct::Side & lSide = quad->side[3]; - int nbhoriz = Min(quad->side[0].NbPoints(), quad->side[2].NbPoints()); - int nbvertic = Min(quad->side[1].NbPoints(), quad->side[3].NbPoints()); + int nbhoriz = Min( bSide.NbPoints(), tSide.NbPoints() ); + int nbvertic = Min( rSide.NbPoints(), lSide.NbPoints() ); if ( myQuadList.size() == 1 ) { // all sub-quads must have NO sides with nbNodeOut > 0 - quad->nbNodeOut(0) = Max( 0, quad->side[0].grid->NbPoints() - quad->side[2].grid->NbPoints()); - quad->nbNodeOut(1) = Max( 0, quad->side[1].grid->NbPoints() - quad->side[3].grid->NbPoints()); - quad->nbNodeOut(2) = Max( 0, quad->side[2].grid->NbPoints() - quad->side[0].grid->NbPoints()); - quad->nbNodeOut(3) = Max( 0, quad->side[3].grid->NbPoints() - quad->side[1].grid->NbPoints()); + quad->nbNodeOut(0) = Max( 0, bSide.grid->NbPoints() - tSide.grid->NbPoints() ); + quad->nbNodeOut(1) = Max( 0, rSide.grid->NbPoints() - lSide.grid->NbPoints() ); + quad->nbNodeOut(2) = Max( 0, tSide.grid->NbPoints() - bSide.grid->NbPoints() ); + quad->nbNodeOut(3) = Max( 0, lSide.grid->NbPoints() - rSide.grid->NbPoints() ); } - int from[4] = { - quad->side[0].from, - quad->side[1].from, - quad->side[2].from, - quad->side[3].from - }; - const vector& uv_e0_vec = quad->side[ 0 ].GetUVPtStruct(); - const vector& uv_e1_vec = quad->side[ 1 ].GetUVPtStruct(); - const vector& uv_e2_vec = quad->side[ 2 ].GetUVPtStruct(); - const vector& uv_e3_vec = quad->side[ 3 ].GetUVPtStruct(); - - if (uv_e0_vec.empty() || uv_e1_vec.empty() || uv_e2_vec.empty() || uv_e3_vec.empty()) + const vector& uv_e0 = bSide.GetUVPtStruct(); + const vector& uv_e1 = rSide.GetUVPtStruct(); + const vector& uv_e2 = tSide.GetUVPtStruct(); + const vector& uv_e3 = lSide.GetUVPtStruct(); + 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); - UVPtStruct* uv_e0 = (UVPtStruct*) & uv_e0_vec[0] + from[0]; - UVPtStruct* uv_e1 = (UVPtStruct*) & uv_e1_vec[0] + from[1]; - UVPtStruct* uv_e2 = (UVPtStruct*) & uv_e2_vec[0] + from[2]; - UVPtStruct* uv_e3 = (UVPtStruct*) & uv_e3_vec[0] + from[3]; - quad->uv_grid.resize( nbvertic * nbhoriz ); quad->iSize = nbhoriz; quad->jSize = nbvertic; @@ -1345,55 +1386,62 @@ bool StdMeshers_Quadrangle_2D::setNormalizedGrid (FaceQuadStruct::Ptr quad) // copy data of face boundary + FaceQuadStruct::SideIterator sideIter; + { // BOTTOM - const int j = 0; - const double x0 = uv_e0[ 0 ].normParam; - const double dx = uv_e0[ nbhoriz-1 ].normParam - uv_e0[ 0 ].normParam; - for (int i = 0; i < nbhoriz; i++) { // down - uv_e0[i].x = ( uv_e0[i].normParam - x0 ) / dx; - uv_e0[i].y = 0.; - uv_grid[ j * nbhoriz + i ] = uv_e0[i]; - quad->uv_box.Add( uv_e0[i].UV() ); + const int j = 0; + const double x0 = bSide.First().normParam; + const double dx = bSide.Last().normParam - bSide.First().normParam; + for ( sideIter.Init( bSide ); sideIter.More(); sideIter.Next() ) { + sideIter.UVPt().x = ( sideIter.UVPt().normParam - x0 ) / dx; + sideIter.UVPt().y = 0.; + uv_grid[ j * nbhoriz + sideIter.Count() ] = sideIter.UVPt(); + quad->uv_box.Add( sideIter.UVPt().UV() ); } } { // RIGHT - const int i = nbhoriz - 1; - const double y0 = uv_e1[ 0 ].normParam; - const double dy = uv_e1[ nbvertic-1 ].normParam - uv_e1[ 0 ].normParam; - int j = 0, nb = nbvertic; - if ( quad->UVPt( i, j ).node ) ++j; // avoid copying from a split emulated side - for ( ; j < nb; j++) { // right - uv_e1[j].x = 1.; - uv_e1[j].y = ( uv_e1[j].normParam - y0 ) / dy; - uv_grid[ j * nbhoriz + i ] = uv_e1[j]; - quad->uv_box.Add( uv_e1[j].UV() ); + const int i = nbhoriz - 1; + const double y0 = rSide.First().normParam; + const double dy = rSide.Last().normParam - rSide.First().normParam; + sideIter.Init( rSide ); + if ( quad->UVPt( i, sideIter.Count() ).node ) + sideIter.Next(); // avoid copying from a split emulated side + for ( ; sideIter.More(); sideIter.Next() ) { + sideIter.UVPt().x = 1.; + sideIter.UVPt().y = ( sideIter.UVPt().normParam - y0 ) / dy; + uv_grid[ sideIter.Count() * nbhoriz + i ] = sideIter.UVPt(); + quad->uv_box.Add( sideIter.UVPt().UV() ); } } { // TOP - const int j = nbvertic - 1; - const double x0 = uv_e2[ 0 ].normParam; - const double dx = uv_e2[ nbhoriz-1 ].normParam - uv_e2[ 0 ].normParam; + const int j = nbvertic - 1; + const double x0 = tSide.First().normParam; + const double dx = tSide.Last().normParam - tSide.First().normParam; int i = 0, nb = nbhoriz; + sideIter.Init( tSide ); if ( quad->UVPt( nb-1, j ).node ) --nb; // avoid copying from a split emulated side - for (; i < nb; i++) { // up - uv_e2[i].x = ( uv_e2[i].normParam - x0 ) / dx; - uv_e2[i].y = 1.; - uv_grid[ j * nbhoriz + i ] = uv_e2[i]; - quad->uv_box.Add( uv_e2[i].UV() ); + for ( ; i < nb; i++, sideIter.Next()) { + sideIter.UVPt().x = ( sideIter.UVPt().normParam - x0 ) / dx; + sideIter.UVPt().y = 1.; + uv_grid[ j * nbhoriz + i ] = sideIter.UVPt(); + quad->uv_box.Add( sideIter.UVPt().UV() ); } } { // LEFT const int i = 0; - const double y0 = uv_e3[ 0 ].normParam; - const double dy = uv_e3[ nbvertic-1 ].normParam - uv_e3[ 0 ].normParam; + const double y0 = lSide.First().normParam; + const double dy = lSide.Last().normParam - lSide.First().normParam; int j = 0, nb = nbvertic; - if ( quad->UVPt( i, j ).node ) ++j; // avoid copying from a split emulated side - if ( quad->UVPt( i, nb-1 ).node ) --nb; - for ( ; j < nb; j++) { // left - uv_e3[j].x = 0.; - uv_e3[j].y = ( uv_e3[j].normParam - y0 ) / dy; - uv_grid[ j * nbhoriz + i ] = uv_e3[j]; - quad->uv_box.Add( uv_e3[j].UV() ); + sideIter.Init( lSide ); + if ( quad->UVPt( i, j ).node ) + ++j, sideIter.Next(); // avoid copying from a split emulated side + if ( quad->UVPt( i, nb-1 ).node ) + --nb; + for ( ; j < nb; j++, sideIter.Next()) { + sideIter.UVPt().x = 0.; + sideIter.UVPt().y = ( sideIter.UVPt().normParam - y0 ) / dy; + uv_grid[ j * nbhoriz + i ] = sideIter.UVPt(); + quad->uv_box.Add( sideIter.UVPt().UV() ); } } @@ -1401,12 +1449,12 @@ bool StdMeshers_Quadrangle_2D::setNormalizedGrid (FaceQuadStruct::Ptr quad) for (int i = 1; i < nbhoriz-1; i++) { - const double x0 = uv_e0[i].x; - const double x1 = uv_e2[i].x; + const double x0 = quad->UVPt( i, 0 ).x; + const double x1 = quad->UVPt( i, nbvertic-1 ).x; for (int j = 1; j < nbvertic-1; j++) { - const double y0 = uv_e3[j].y; - const double y1 = uv_e1[j].y; + const double y0 = quad->UVPt( 0, j ).y; + const double y1 = quad->UVPt( nbhoriz-1, j ).y; // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0) double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); double y = y0 + x * (y1 - y0); @@ -1419,19 +1467,19 @@ bool StdMeshers_Quadrangle_2D::setNormalizedGrid (FaceQuadStruct::Ptr quad) // projection on 2d domain (u,v) - gp_UV a0 = uv_e0[0 ].UV(); - gp_UV a1 = uv_e0[nbhoriz-1].UV(); - gp_UV a2 = uv_e2[nbhoriz-1].UV(); - gp_UV a3 = uv_e2[0 ].UV(); + gp_UV a0 = quad->UVPt( 0, 0 ).UV(); + gp_UV a1 = quad->UVPt( nbhoriz-1, 0 ).UV(); + gp_UV a2 = quad->UVPt( nbhoriz-1, nbvertic-1 ).UV(); + gp_UV a3 = quad->UVPt( 0, nbvertic-1 ).UV(); for (int i = 1; i < nbhoriz-1; i++) { - gp_UV p0 = uv_e0[i].UV(); - gp_UV p2 = uv_e2[i].UV(); + gp_UV p0 = quad->UVPt( i, 0 ).UV(); + gp_UV p2 = quad->UVPt( i, nbvertic-1 ).UV(); for (int j = 1; j < nbvertic-1; j++) { - gp_UV p1 = uv_e1[j].UV(); - gp_UV p3 = uv_e3[j].UV(); + gp_UV p1 = quad->UVPt( nbhoriz-1, j ).UV(); + gp_UV p3 = quad->UVPt( 0, j ).UV(); int ij = j * nbhoriz + i; double x = uv_grid[ij].x; @@ -1451,31 +1499,55 @@ bool StdMeshers_Quadrangle_2D::setNormalizedGrid (FaceQuadStruct::Ptr quad) //purpose : auxilary function for computeQuadPref //======================================================================= -static void shiftQuad(FaceQuadStruct::Ptr& quad, const int num) +void StdMeshers_Quadrangle_2D::shiftQuad(FaceQuadStruct::Ptr& quad, const int num ) { - quad->shift( num, /*ori=*/true ); + quad->shift( num, /*ori=*/true, /*keepGrid=*/myQuadList.size() > 1 ); } //================================================================================ /*! - * \brief Rotate sides of a quad by nb + * \brief Rotate sides of a quad by given nb of quartes * \param nb - number of rotation quartes * \param ori - to keep orientation of sides as in an unit quad or not + * \param keepGrid - if \c true Side::grid is not changed, Side::from and Side::to + * are altered instead */ //================================================================================ -void FaceQuadStruct::shift( size_t nb, bool ori ) +void FaceQuadStruct::shift( size_t nb, bool ori, bool keepGrid ) { if ( nb == 0 ) return; - StdMeshers_FaceSidePtr sideArr[4] = { side[0], side[1], side[2], side[3] }; - for (int i = QUAD_BOTTOM_SIDE; i < NB_QUAD_SIDES; ++i) { + + vector< Side > newSides( side.size() ); + vector< Side* > sidePtrs( side.size() ); + for (int i = QUAD_BOTTOM_SIDE; i < NB_QUAD_SIDES; ++i) + { int id = (i + nb) % NB_QUAD_SIDES; - bool wasForward = (i < QUAD_TOP_SIDE); - bool newForward = (id < QUAD_TOP_SIDE); - if (ori && wasForward != newForward) - sideArr[ i ]->Reverse(); - side[ id ] = sideArr[ i ]; + if ( ori ) + { + bool wasForward = (i < QUAD_TOP_SIDE); + bool newForward = (id < QUAD_TOP_SIDE); + if ( wasForward != newForward ) + side[ i ].Reverse( keepGrid ); + } + newSides[ id ] = side[ i ]; + sidePtrs[ i ] = & side[ i ]; } + // make newSides refer newSides via Side::Contact's + for ( size_t i = 0; i < newSides.size(); ++i ) + { + FaceQuadStruct::Side& ns = newSides[ i ]; + for ( size_t iC = 0; iC < ns.contacts.size(); ++iC ) + { + FaceQuadStruct::Side* oSide = ns.contacts[iC].other_side; + vector< Side* >::iterator sIt = std::find( sidePtrs.begin(), sidePtrs.end(), oSide ); + if ( sIt != sidePtrs.end() ) + ns.contacts[iC].other_side = & newSides[ *sIt - sidePtrs[0] ]; + } + } + newSides.swap( side ); + + uv_grid.clear(); } //======================================================================= @@ -1532,20 +1604,17 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, const TopoDS_Face& aFace, FaceQuadStruct::Ptr quad) { - // Auxilary key in order to keep old variant - // of meshing after implementation new variant - // for bug 0016220 from Mantis. - bool OldVersion = (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED); + const bool OldVersion = (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED); + const bool WisF = true; SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); Handle(Geom_Surface) S = BRep_Tool::Surface(aFace); - bool WisF = true; - int i,j,geomFaceID = meshDS->ShapeToIndex(aFace); + int i,j, geomFaceID = meshDS->ShapeToIndex(aFace); - int nb = quad->side[0].grid->NbPoints(); - int nr = quad->side[1].grid->NbPoints(); - int nt = quad->side[2].grid->NbPoints(); - int nl = quad->side[3].grid->NbPoints(); + 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); @@ -1562,17 +1631,17 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, { // rotate the quad to have nt > nb [and nr > nl] if ( nb > nt ) - quad->shift( nr > nl ? 1 : 2, true ); + shiftQuad ( quad, nr > nl ? 1 : 2 ); else if ( nr > nl ) - quad->shift( nb == nt ? 1 : 0, true ); + shiftQuad( quad, nb == nt ? 1 : 0 ); else if ( nl > nr ) - quad->shift( 3, true ); + shiftQuad( quad, 3 ); } - nb = quad->side[0].grid->NbPoints(); - nr = quad->side[1].grid->NbPoints(); - nt = quad->side[2].grid->NbPoints(); - nl = quad->side[3].grid->NbPoints(); + 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); @@ -1607,13 +1676,27 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, // 0------------0 // 0 bottom 1 - 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); - if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl) - return error(COMPERR_BAD_INPUT_MESH); + const int bfrom = quad->side[0].from; + const int rfrom = quad->side[1].from; + const int tfrom = quad->side[2].from; + const int lfrom = quad->side[3].from; + { + const vector& uv_eb_vec = quad->side[0].GetUVPtStruct(true,0); + const vector& uv_er_vec = quad->side[1].GetUVPtStruct(false,1); + const vector& uv_et_vec = quad->side[2].GetUVPtStruct(true,1); + const vector& uv_el_vec = quad->side[3].GetUVPtStruct(false,0); + if (uv_eb_vec.empty() || + uv_er_vec.empty() || + uv_et_vec.empty() || + uv_el_vec.empty()) + return error(COMPERR_BAD_INPUT_MESH); + } + FaceQuadStruct::SideIterator uv_eb, uv_er, uv_et, uv_el; + uv_eb.Init( quad->side[0] ); + uv_er.Init( quad->side[1] ); + uv_et.Init( quad->side[2] ); + uv_el.Init( quad->side[3] ); gp_UV a0,a1,a2,a3, p0,p1,p2,p3, uv; double x,y; @@ -1625,7 +1708,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, if ( !myForcedPnts.empty() ) { - if ( dv != 0 && dh != 0 ) + if ( dv != 0 && dh != 0 ) // here myQuadList.size() == 1 { const int dmin = Min( dv, dh ); @@ -1659,10 +1742,11 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, x = uv_et[ dmin ].normParam; p0 = quad->side[0].grid->Value2d( x ).XY(); p2 = uv_et[ dmin ].UV(); - for ( int i = 1; i < nl; ++i ) + double y0 = uv_er[ dmin ].normParam; + for ( int i = 1; i < nl-1; ++i ) { - y = uv_er[ i + dmin ].normParam; - p1 = uv_er[ i + dmin ].UV(); + y = y0 + i / ( nl-1. ) * ( 1. - y0 ); + p1 = quad->side[1].grid->Value2d( y ).XY(); p3 = quad->side[3].grid->Value2d( y ).XY(); uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); pointsLCt[ i ].u = uv.X(); @@ -1679,10 +1763,11 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, y = uv_er[ dmin ].normParam; p1 = uv_er[ dmin ].UV(); p3 = quad->side[3].grid->Value2d( y ).XY(); + double x0 = uv_et[ dmin ].normParam; for ( int i = 1; i < nb-1; ++i ) { - x = uv_et[ i + dmin ].normParam; - p2 = uv_et[ i + dmin ].UV(); + x = x0 + i / ( nb-1. ) * ( 1. - x0 ); + p2 = quad->side[2].grid->Value2d( x ).XY(); p0 = quad->side[0].grid->Value2d( x ).XY(); uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); pointsCbCt[ i ].u = uv.X(); @@ -1691,7 +1776,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, sideCbCt = StdMeshers_FaceSide::New( pointsCbCt, aFace ); } // Make Cb quad - FaceQuadStruct* qCb = new FaceQuadStruct( quad->face ); + FaceQuadStruct* qCb = new FaceQuadStruct( quad->face, "Cb" ); myQuadList.push_back( FaceQuadStruct::Ptr( qCb )); qCb->side.resize(4); qCb->side[0] = quad->side[0]; @@ -1700,7 +1785,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, qCb->side[3] = sideLCb; qCb->side[1].to = dmin+1; // Make L quad - FaceQuadStruct* qL = new FaceQuadStruct( quad->face ); + FaceQuadStruct* qL = new FaceQuadStruct( quad->face, "L" ); myQuadList.push_back( FaceQuadStruct::Ptr( qL )); qL->side.resize(4); qL->side[0] = sideLCb; @@ -1715,6 +1800,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, qCt->side[1].from = dmin; qCt->side[2].from = dmin; qCt->uv_grid.clear(); + qCt->name = "Ct"; // Connect sides qCb->side[3].AddContact( dmin, & qCb->side[2], 0 ); @@ -1731,7 +1817,12 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, } // if ( dv != 0 && dh != 0 ) - // Case dv == 0 + const int db = quad->side[0].IsReversed() ? -1 : +1; + const int dr = quad->side[1].IsReversed() ? -1 : +1; + const int dt = quad->side[2].IsReversed() ? -1 : +1; + const int dl = quad->side[3].IsReversed() ? -1 : +1; + + // Case dv == 0, here possibly myQuadList.size() > 1 // // lw nb lw = dh/2 // +------------+ @@ -1744,33 +1835,33 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, // |/ \| // +------------+ const int lw = dh/2; // lateral width - const int bfrom = quad->side[0].from; - const int rfrom = quad->side[1].from; - const int tfrom = quad->side[2].from; - const int lfrom = quad->side[3].from; - - const double lL = quad->side[3].Length(); - const double lLwL = quad->side[2].Length( tfrom, tfrom + lw + 1 ); - const double yCbL = lLwL / ( lLwL + lL ); - - const double lR = quad->side[1].Length(); - const double lLwR = quad->side[2].Length( nt - lw - 1, nt ); - const double yCbR = lLwR / ( lLwR + lR ); + double yCbL, yCbR; + { + double lL = quad->side[3].Length(); + double lLwL = quad->side[2].Length( tfrom, + tfrom + ( lw ) * dt ); + yCbL = lLwL / ( lLwL + lL ); + + double lR = quad->side[1].Length(); + double lLwR = quad->side[2].Length( tfrom + ( lw + nb-1 ) * dt, + tfrom + ( lw + nb-1 + lw ) * dt); + yCbR = lLwR / ( lLwR + lR ); + } // Make sides separating domains Cb and L and R StdMeshers_FaceSidePtr sideLCb, sideRCb; UVPtStruct pTBL, pTBR; // points where 3 domains meat { vector pointsLCb( lw+1 ), pointsRCb( lw+1 ); - pointsLCb[0] = uv_eb[ 0 + bfrom ]; - pointsRCb[0] = uv_eb[ nb + bfrom ]; + pointsLCb[0] = uv_eb[ 0 ]; + pointsRCb[0] = uv_eb[ nb-1 ]; for ( int i = 1, i2 = nt-2; i <= lw; ++i, --i2 ) { x = quad->side[2].Param( i ); y = yCbL * i / lw; p0 = quad->side[0].Value2d( x ); p1 = quad->side[1].Value2d( y ); - p2 = uv_et[ i + tfrom ].UV(); + p2 = uv_et[ i ].UV(); p3 = quad->side[3].Value2d( y ); uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); pointsLCb[ i ].u = uv.X(); @@ -1779,8 +1870,10 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, x = quad->side[2].Param( i2 ); y = yCbR * i / lw; + p1 = quad->side[1].Value2d( y ); p0 = quad->side[0].Value2d( x ); - p2 = uv_et[ i2 + tfrom ].UV(); + p2 = uv_et[ i2 ].UV(); + p3 = quad->side[3].Value2d( y ); uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); pointsRCb[ i ].u = uv.X(); pointsRCb[ i ].v = uv.Y(); @@ -1796,29 +1889,29 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, { vector pointsLCt( nl ), pointsRCt( nl ); pointsLCt[0] = pTBL; - pointsLCt.back() = uv_et[ lw + tfrom ]; + pointsLCt.back() = uv_et[ lw ]; pointsRCt[0] = pTBR; - pointsRCt.back() = uv_et[ lw + nb - 1 + tfrom ]; + pointsRCt.back() = uv_et[ lw + nb - 1 ]; x = pTBL.x; p0 = quad->side[0].Value2d( x ); - p2 = uv_et[ lw + tfrom ].UV(); + p2 = uv_et[ lw ].UV(); int iR = lw + nb - 1; double xR = pTBR.x; gp_UV p0R = quad->side[0].Value2d( xR ); - gp_UV p2R = uv_et[ iR + tfrom ].UV(); - for ( int i = 1; i < nl; ++i ) + gp_UV p2R = uv_et[ iR ].UV(); + for ( int i = 1; i < nl-1; ++i ) { - y = yCbL + ( 1. - yCbL ) * i / nl; + y = yCbL + ( 1. - yCbL ) * i / (nl-1.); p1 = quad->side[1].Value2d( y ); p3 = quad->side[3].Value2d( y ); uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); pointsLCt[ i ].u = uv.X(); pointsLCt[ i ].v = uv.Y(); - y = yCbR + ( 1. - yCbR ) * i / nl; + y = yCbR + ( 1. - yCbR ) * i / (nl-1.); p1 = quad->side[1].Value2d( y ); p3 = quad->side[3].Value2d( y ); - uv = calcUV( xR,y, a0,a1,a2,a3, p0,p1,p2,p3 ); + uv = calcUV( xR,y, a0,a1,a2,a3, p0R,p1,p2R,p3 ); pointsRCt[ i ].u = uv.X(); pointsRCt[ i ].v = uv.Y(); } @@ -1836,8 +1929,8 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, for ( int i = 1; i < nb-1; ++i ) { x = quad->side[2].Param( i + lw ); - y = yCbL + ( yCbR - yCbL ) * i / nb; - p2 = uv_et[ i + lw + tfrom ].UV(); + y = yCbL + ( yCbR - yCbL ) * i / (nb-1.); + p2 = uv_et[ i + lw ].UV(); p0 = quad->side[0].Value2d( x ); uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 ); pointsCbCt[ i ].u = uv.X(); @@ -1846,7 +1939,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, sideCbCt = StdMeshers_FaceSide::New( pointsCbCt, aFace ); } // Make Cb quad - FaceQuadStruct* qCb = new FaceQuadStruct( quad->face ); + FaceQuadStruct* qCb = new FaceQuadStruct( quad->face, "Cb" ); myQuadList.push_back( FaceQuadStruct::Ptr( qCb )); qCb->side.resize(4); qCb->side[0] = quad->side[0]; @@ -1854,54 +1947,54 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, qCb->side[2] = sideCbCt; qCb->side[3] = sideLCb; // Make L quad - FaceQuadStruct* qL = new FaceQuadStruct( quad->face ); + FaceQuadStruct* qL = new FaceQuadStruct( quad->face, "L" ); myQuadList.push_back( FaceQuadStruct::Ptr( qL )); qL->side.resize(4); qL->side[0] = sideLCb; qL->side[1] = sideLCt; qL->side[2] = quad->side[2]; qL->side[3] = quad->side[3]; - qL->side[2].to = lw+1; + qL->side[2].to = ( lw + 1 ) * dt + tfrom; // Make R quad - FaceQuadStruct* qR = new FaceQuadStruct( quad->face ); + FaceQuadStruct* qR = new FaceQuadStruct( quad->face, "R" ); myQuadList.push_back( FaceQuadStruct::Ptr( qR )); qR->side.resize(4); qR->side[0] = sideRCb; qR->side[0].from = lw; qR->side[0].to = -1; + qR->side[0].di = -1; qR->side[1] = quad->side[1]; qR->side[2] = quad->side[2]; - qR->side[2].from = nb + lw + tfrom; + qR->side[2].from = ( lw + nb-1 ) * dt + tfrom; qR->side[3] = sideRCt; // Make Ct from the main quad FaceQuadStruct::Ptr qCt = quad; qCt->side[0] = sideCbCt; qCt->side[1] = sideRCt; - qCt->side[2].from = lw + tfrom; - qCt->side[2].to = nt - lw + tfrom; + qCt->side[2].from = ( lw ) * dt + tfrom; + qCt->side[2].to = ( lw + nb ) * dt + tfrom; qCt->side[3] = sideLCt; qCt->uv_grid.clear(); + qCt->name = "Ct"; // Connect sides qCb->side[3].AddContact( lw, & qCb->side[2], 0 ); qCb->side[3].AddContact( lw, & qCt->side[3], 0 ); - qCt->side[3].AddContact( 0, & qCt->side[0], 0 ); - qCt->side[0].AddContact( 0, & qL ->side[0], lw ); + qCt->side[3].AddContact( 0, & qCt->side[0], 0 ); + qCt->side[0].AddContact( 0, & qL ->side[0], lw ); qL ->side[0].AddContact( lw, & qL ->side[1], 0 ); qL ->side[0].AddContact( lw, & qCb->side[2], 0 ); // - qCb->side[1].AddContact( lw, & qCb->side[2], lw ); - qCb->side[1].AddContact( lw, & qCt->side[1], 0 ); - qCt->side[0].AddContact( lw, & qCt->side[1], 0 ); - qCt->side[0].AddContact( lw, & qR ->side[0], lw ); - qR ->side[3].AddContact( lw, & qR ->side[0], lw ); - qR ->side[3].AddContact( lw, & qCb->side[2], lw ); + qCb->side[1].AddContact( lw, & qCb->side[2], nb-1 ); + qCb->side[1].AddContact( lw, & qCt->side[1], 0 ); + qCt->side[0].AddContact( nb-1, & qCt->side[1], 0 ); + qCt->side[0].AddContact( nb-1, & qR ->side[0], lw ); + qR ->side[3].AddContact( 0, & qR ->side[0], lw ); + qR ->side[3].AddContact( 0, & qCb->side[2], nb-1 ); - if ( dh == dv ) - return computeQuadDominant( aMesh, aFace ); + return computeQuadDominant( aMesh, aFace ); - - } + } // if ( !myForcedPnts.empty() ) if ( dh > dv ) { addv = (dh-dv)/2; @@ -1998,12 +2091,6 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, NodesL.Value(i+1,j+1), NodesL.Value(i,j+1)); if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1), - NodesL.Value(i+1,j+1), NodesL.Value(i+1,j)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); - } } } } @@ -2061,12 +2148,6 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, NodesR.Value(i+1,j+1), NodesR.Value(i,j+1)); if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1), - NodesR.Value(i+1,j+1), NodesR.Value(i+1,j)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); - } } } } @@ -2140,12 +2221,6 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, NodesC.Value(i+1,j+1), NodesC.Value(i,j+1)); if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1), - NodesC.Value(i+1,j+1), NodesC.Value(i+1,j)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); - } } } } @@ -2178,12 +2253,6 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, NodesBRD.Value(i+1,j+1), NodesBRD.Value(i,j+1)); if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i,j+1), - NodesBRD.Value(i+1,j+1), NodesBRD.Value(i+1,j)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); - } } } int drl = abs(nr-nl); @@ -2290,12 +2359,6 @@ bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh, NodesC.Value(i+1,j+1), NodesC.Value(i,j+1)); if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1), - NodesC.Value(i+1,j+1), NodesC.Value(i+1,j)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); - } } } // end nrSetMeshElementOnShape(F, geomFaceID); } - else { - SMDS_MeshFace* F = - myHelper->AddFace(NodesLast.Value(i,1), NodesLast.Value(i,2), - NodesLast.Value(i+1,2), NodesLast.Value(i+1,2)); - if (F) meshDS->SetMeshElementOnShape(F, geomFaceID); - } } } // if ((drl+addv) > 0) @@ -2695,10 +2752,10 @@ bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, Handle(Geom_Surface) S = BRep_Tool::Surface(aFace); int i,j,geomFaceID = meshDS->ShapeToIndex(aFace); - int nb = quad->side[0].grid->NbPoints(); // bottom - int nr = quad->side[1].grid->NbPoints(); // right - int nt = quad->side[2].grid->NbPoints(); // top - int nl = quad->side[3].grid->NbPoints(); // left + int nb = quad->side[0].NbPoints(); // bottom + int nr = quad->side[1].NbPoints(); // right + int nt = quad->side[2].NbPoints(); // top + int nl = quad->side[3].NbPoints(); // left // Simple Reduce 10->8->6->4 (3 steps) Multiple Reduce 10->4 (1 step) // @@ -2787,10 +2844,10 @@ bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, } } - nb = quad->side[0].grid->NbPoints(); - nr = quad->side[1].grid->NbPoints(); - nt = quad->side[2].grid->NbPoints(); - nl = quad->side[3].grid->NbPoints(); + 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); @@ -3061,10 +3118,10 @@ bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh, } } - nb = quad->side[0].grid->NbPoints(); - nr = quad->side[1].grid->NbPoints(); - nt = quad->side[2].grid->NbPoints(); - nl = quad->side[3].grid->NbPoints(); + nb = quad->side[0].NbPoints(); + nr = quad->side[1].NbPoints(); + nt = quad->side[2].NbPoints(); + nl = quad->side[3].NbPoints(); // number of rows and columns int nrows = nr - 1; // and also == nl - 1 @@ -3631,6 +3688,18 @@ namespace // data for smoothing double d = v1 ^ v2; return d > 1e-100; } + //================================================================================ + /*! + * \brief Returns area of a triangle + */ + //================================================================================ + + double getArea( const gp_UV uv1, const gp_UV uv2, const gp_UV uv3 ) + { + gp_XY v1 = uv1 - uv2, v2 = uv3 - uv2; + double a = v2 ^ v1; + return a; + } } //================================================================================ @@ -3718,6 +3787,8 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad) // Get nodes to smooth + // TODO: do not smooth fixed nodes + typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap; TNo2SmooNoMap smooNoMap; @@ -3871,6 +3942,144 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad) } } +//================================================================================ +/*! + * \brief Checks validity of generated faces + */ +//================================================================================ + +bool StdMeshers_Quadrangle_2D::check() +{ + const bool isOK = true; + if ( !myCheckOri || myQuadList.empty() || !myQuadList.front() || !myHelper ) + return isOK; + + TopoDS_Face geomFace = TopoDS::Face( myHelper->GetSubShape() ); + SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); + SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace ); + bool toCheckUV; + if ( geomFace.Orientation() >= TopAbs_INTERNAL ) geomFace.Orientation( TopAbs_FORWARD ); + + // Get a reference orientation sign + + double okSign; + { + TError err; + TSideVector wireVec = + StdMeshers_FaceSide::GetFaceWires( geomFace, *myHelper->GetMesh(), true, err ); + StdMeshers_FaceSidePtr wire = wireVec[0]; + + // find a right angle VERTEX + int iVertex; + double maxAngle = -1e100; + for ( int i = 0; i < wire->NbEdges(); ++i ) + { + int iPrev = myHelper->WrapIndex( i-1, wire->NbEdges() ); + const TopoDS_Edge& e1 = wire->Edge( iPrev ); + const TopoDS_Edge& e2 = wire->Edge( i ); + double angle = myHelper->GetAngle( e1, e2, geomFace ); + if ( maxAngle < angle && angle < 0.9 * M_PI ) + { + maxAngle = angle; + iVertex = i; + } + } + if ( maxAngle < -2*M_PI ) return isOK; + + // get a sign of 2D area of a corner face + + int iPrev = myHelper->WrapIndex( iVertex-1, wire->NbEdges() ); + const TopoDS_Edge& e1 = wire->Edge( iPrev ); + const TopoDS_Edge& e2 = wire->Edge( iVertex ); + + gp_Vec2d v1, v2; gp_Pnt2d p; + double u[2]; + { + bool rev = ( e1.Orientation() == TopAbs_REVERSED ); + Handle(Geom2d_Curve) c = BRep_Tool::CurveOnSurface( e1, geomFace, u[0], u[1] ); + c->D1( u[ !rev ], p, v1 ); + if ( !rev ) + v1.Reverse(); + } + { + bool rev = ( e2.Orientation() == TopAbs_REVERSED ); + Handle(Geom2d_Curve) c = BRep_Tool::CurveOnSurface( e2, geomFace, u[0], u[1] ); + c->D1( u[ rev ], p, v2 ); + if ( rev ) + v2.Reverse(); + } + + okSign = v2 ^ v1; + + if ( maxAngle < 0 ) + okSign *= -1; + } + + // Look for incorrectly oriented faces + + std::list badFaces; + + const SMDS_MeshNode* nn [ 8 ]; // 8 is just for safety + gp_UV uv [ 8 ]; + SMDS_ElemIteratorPtr fIt = fSubMesh->GetElements(); + while ( fIt->more() ) // loop on faces bound to a FACE + { + const SMDS_MeshElement* f = fIt->next(); + + const int nbN = f->NbCornerNodes(); + for ( int i = 0; i < nbN; ++i ) + nn[ i ] = f->GetNode( i ); + + const SMDS_MeshNode* nInFace = 0; + if ( myHelper->HasSeam() ) + for ( int i = 0; i < nbN && !nInFace; ++i ) + if ( !myHelper->IsSeamShape( nn[i]->getshapeId() )) + nInFace = nn[i]; + + for ( int i = 0; i < nbN; ++i ) + uv[ i ] = myHelper->GetNodeUV( geomFace, nn[i], nInFace, &toCheckUV ); + + switch ( nbN ) { + case 4: + { + double sign1 = getArea( uv[0], uv[1], uv[2] ); + double sign2 = getArea( uv[0], uv[2], uv[3] ); + if ( sign1 * sign2 < 0 ) + { + sign2 = getArea( uv[1], uv[2], uv[3] ); + sign1 = getArea( uv[1], uv[3], uv[0] ); + if ( sign1 * sign2 < 0 ) + continue; // this should not happen + } + if ( sign1 * okSign < 0 ) + badFaces.push_back ( f ); + break; + } + case 3: + { + double sign = getArea( uv[0], uv[1], uv[2] ); + if ( sign * okSign < 0 ) + badFaces.push_back ( f ); + break; + } + default:; + } + } + + if ( !badFaces.empty() ) + { + SMESH_subMesh* fSM = myHelper->GetMesh()->GetSubMesh( geomFace ); + SMESH_ComputeErrorPtr& err = fSM->GetComputeError(); + err.reset ( new SMESH_ComputeError( COMPERR_ALGO_FAILED, + "Inverted elements generated")); + err->myBadElements.swap( badFaces ); + + return !isOK; + } + + return isOK; +} + /*//================================================================================ /*! * \brief Finds vertices at the most sharp face corners @@ -3933,7 +4142,8 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, TopoDS_Shape triaVertex = helper.GetMeshDS()->IndexToShape( myTriaVertexID ); if ( !triaVertex.IsNull() && triaVertex.ShapeType() == TopAbs_VERTEX && - helper.IsSubShape( triaVertex, theFace )) + helper.IsSubShape( triaVertex, theFace ) && + ( vertexByAngle.size() != 4 || vertexByAngle.begin()->first < 5 * M_PI/180. )) nbCorners = 3; else triaVertex.Nullify(); @@ -3984,6 +4194,9 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, isThereVariants = ( lostAngle * 1.1 >= lastAngle ); } + myCheckOri = ( vertexByAngle.size() > nbCorners || + vertexByAngle.begin()->first < 5.* M_PI/180 ); + // make theWire begin from a corner vertex or triaVertex if ( nbCorners == 3 ) while ( !triaVertex.IsSame( ( helper.IthVertex( 0, theWire.front() ))) || @@ -4124,7 +4337,7 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face& theFace, //================================================================================ FaceQuadStruct::Side::Side(StdMeshers_FaceSidePtr theGrid) - : grid(theGrid), nbNodeOut(0), from(0), to(theGrid ? theGrid->NbPoints() : 0 ) + : grid(theGrid), nbNodeOut(0), from(0), to(theGrid ? theGrid->NbPoints() : 0 ), di(1) { } @@ -4134,7 +4347,8 @@ FaceQuadStruct::Side::Side(StdMeshers_FaceSidePtr theGrid) */ //============================================================================= -FaceQuadStruct::FaceQuadStruct(const TopoDS_Face& F) : face( F ) +FaceQuadStruct::FaceQuadStruct(const TopoDS_Face& F, const std::string& theName) + : face( F ), name( theName ) { side.reserve(4); } @@ -4185,6 +4399,9 @@ bool StdMeshers_Quadrangle_2D::getEnforcedUV() surf->Bounds( u1,u2,v1,v2 ); GeomAPI_ProjectPointOnSurf project; project.Init(surf, u1,u2, v1,v2, tol ); + Bnd_Box bbox; + BRepBndLib::Add( face, bbox ); + double farTol = 0.01 * sqrt( bbox.SquareExtent() ); for ( size_t iP = 0; iP < points.size(); ++iP ) { @@ -4197,7 +4414,7 @@ bool StdMeshers_Quadrangle_2D::getEnforcedUV() << points[ iP ].X() << ", "<< points[ iP ].Y() << ", "<< points[ iP ].Z() << " )"); continue; } - if ( project.LowerDistance() > tol*1000 ) + if ( project.LowerDistance() > farTol ) { if ( isStrictCheck && iP < nbPoints ) return error @@ -4274,7 +4491,11 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() list< FaceQuadStruct::Ptr >::iterator quadIt = myQuadList.begin(); for ( ; quadIt != myQuadList.end(); ++quadIt ) for ( size_t iSide = 0; iSide < (*quadIt)->side.size(); ++iSide ) + { + if ( !setNormalizedGrid( *quadIt )) + return false; quadsBySide[ (*quadIt)->side[iSide] ].push_back( *quadIt ); + } SMESH_Mesh* mesh = myHelper->GetMesh(); SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); @@ -4289,9 +4510,9 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() for ( quadIt = myQuadList.begin(); quadIt != myQuadList.end(); ++quadIt ) { FaceQuadStruct::Ptr quad = *quadIt; - int i,j; - if ( !setNormalizedGrid( quad )) + if ( !setNormalizedGrid( *quadIt )) return false; + int i,j; if ( !quad->findCell( myForcedPnts[ iFP ], i, j )) continue; @@ -4302,25 +4523,25 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() for ( int dj = 0; dj < 2; ++dj ) { double dist2 = ( myForcedPnts[ iFP ].uv - quad->UVPt( i+di,j+dj ).UV() ).SquareModulus(); - ijByDist.insert( make_pair( dist2, make_pair( i+di,j+dj ))); + ijByDist.insert( make_pair( dist2, make_pair( di,dj ))); } // try all nodes starting from the closest one set< FaceQuadStruct::Ptr > changedQuads; multimap< double, pair< int, int > >::iterator d2ij = ijByDist.begin(); for ( ; !isNodeEnforced && d2ij != ijByDist.end(); ++d2ij ) { - i = d2ij->second.first; - j = d2ij->second.second; + int di = d2ij->second.first; + int dj = d2ij->second.second; // check if a node is at a side int iSide = -1; - if ( j == 0 ) + if ( dj== 0 && j == 0 ) iSide = QUAD_BOTTOM_SIDE; - else if ( j+1 == quad->jSize ) + else if ( dj == 1 && j+2 == quad->jSize ) iSide = QUAD_TOP_SIDE; - else if ( i == 0 ) + else if ( di == 0 && i == 0 ) iSide = QUAD_LEFT_SIDE; - else if ( i+1 == quad->iSize ) + else if ( di == 1 && i+2 == quad->iSize ) iSide = QUAD_RIGHT_SIDE; if ( iSide > -1 ) // ----- node is at a side @@ -4369,18 +4590,20 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() } else // ------------------ node is inside the quad { + i += di; + j += dj; // make a new side passing through IJ node and split the quad int indForced, iNewSide; if ( quad->iSize < quad->jSize ) // split vertically { quad->updateUV( myForcedPnts[ iFP ].uv, i, j, /*isVert=*/true ); - indForced = i; + indForced = j; iNewSide = splitQuad( quad, i, 0 ); } else { quad->updateUV( myForcedPnts[ iFP ].uv, i, j, /*isVert=*/false ); - indForced = j; + indForced = i; iNewSide = splitQuad( quad, 0, j ); } FaceQuadStruct::Ptr newQuad = myQuadList.back(); @@ -4390,7 +4613,10 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() quad->side[( iNewSide+2 ) % 4 ].forced_nodes.insert( indForced ); quadsBySide[ newSide ].push_back( quad ); - quadsBySide[ newSide ].push_back( newQuad ); + quadsBySide[ newQuad->side[0] ].push_back( newQuad ); + quadsBySide[ newQuad->side[1] ].push_back( newQuad ); + quadsBySide[ newQuad->side[2] ].push_back( newQuad ); + quadsBySide[ newQuad->side[3] ].push_back( newQuad ); isNodeEnforced = true; @@ -4403,6 +4629,9 @@ bool StdMeshers_Quadrangle_2D::addEnforcedNodes() for ( ; qIt != changedQuads.end(); ++qIt ) (*qIt)->uv_grid.clear(); + if ( isNodeEnforced ) + break; + } // loop on quads if ( !isNodeEnforced ) @@ -4509,6 +4738,7 @@ int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J) newQuad->side[ QUAD_BOTTOM_SIDE ].from = iBot; newQuad->side[ QUAD_TOP_SIDE ].from = iTop; + newQuad->name = ( TComm("Right of I=") << I ); quad->side[ QUAD_BOTTOM_SIDE ].to = iBot + 1; quad->side[ QUAD_TOP_SIDE ].to = iTop + 1; @@ -4516,7 +4746,7 @@ int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J) return QUAD_LEFT_SIDE; } - else if ( J > 0 ) //// split horizontally + else if ( J > 0 ) //// split horizontally, a new quad is below an old one { points.reserve( quad->iSize ); for ( int iP = 0; iP < quad->iSize; ++iP ) @@ -4547,6 +4777,7 @@ int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J) newQuad->side[ QUAD_RIGHT_SIDE ].to = iRgt+1; newQuad->side[ QUAD_LEFT_SIDE ].to = iLft+1; + newQuad->name = ( TComm("Below J=") << J ); quad->side[ QUAD_RIGHT_SIDE ].from = iRgt; quad->side[ QUAD_LEFT_SIDE ].from = iLft; @@ -4554,6 +4785,9 @@ int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J) return QUAD_TOP_SIDE; } + + myQuadList.pop_back(); + return -1; } //================================================================================ @@ -4590,10 +4824,10 @@ void StdMeshers_Quadrangle_2D::updateSideUV( FaceQuadStruct::Side& side, const int iTo = Max ( iForced, *iNext ) + 1; const int sideSize = iTo - iFrom; - vector points[4]; + vector points[4]; // side points of a temporary quad - // get from the quads grid points adjacent to the side - // to make two sides of another temporary quad + // from the quads get grid points adjacent to the side + // to make two sides of a temporary quad vector< FaceQuadStruct::Ptr > quads = quadsBySide.find( side )->second; // copy! for ( int is2nd = 0; is2nd < 2; ++is2nd ) { @@ -4607,7 +4841,8 @@ void StdMeshers_Quadrangle_2D::updateSideUV( FaceQuadStruct::Side& side, for ( size_t iQ = 0; iQ < quads.size(); ++iQ ) { FaceQuadStruct::Ptr q = quads[ iQ ]; - if ( !q ) continue; + if ( !q ) + continue; size_t iS; for ( iS = 0; iS < q->side.size(); ++iS ) if ( side.grid == q->side[ iS ].grid ) @@ -4619,10 +4854,12 @@ void StdMeshers_Quadrangle_2D::updateSideUV( FaceQuadStruct::Side& side, isOut = ( q->side[ iS ].to >= iCur || q->side[ iS ].from <= iCur ); if ( isOut ) continue; + if ( !setNormalizedGrid( q )) + continue; // found - copy points int i,j,di,dj,nb; - if ( iS % 2 ) // right ot left + if ( iS % 2 ) // right or left { i = ( iS == QUAD_LEFT_SIDE ) ? 1 : q->iSize-2; j = q->side[ iS ].ToQuadIndex( iCur ); @@ -4678,7 +4915,8 @@ void StdMeshers_Quadrangle_2D::updateSideUV( FaceQuadStruct::Side& side, points[T].push_back( points[R].back() ); // make the temporary quad - FaceQuadStruct::Ptr tmpQuad( new FaceQuadStruct( TopoDS::Face( myHelper->GetSubShape() ))); + FaceQuadStruct::Ptr tmpQuad + ( new FaceQuadStruct( TopoDS::Face( myHelper->GetSubShape() ), "tmpQuad")); tmpQuad->side.push_back( StdMeshers_FaceSide::New( points[B] )); // bottom tmpQuad->side.push_back( StdMeshers_FaceSide::New( points[R] )); // right tmpQuad->side.push_back( StdMeshers_FaceSide::New( points[T] )); @@ -4692,7 +4930,11 @@ void StdMeshers_Quadrangle_2D::updateSideUV( FaceQuadStruct::Side& side, // update UV of the side vector& sidePoints = (vector&) side.GetUVPtStruct(); for ( int i = iFrom; i < iTo; ++i ) - sidePoints[ i ] = tmpQuad->UVPt( 1, i-iFrom ); + { + const uvPtStruct& uvPt = tmpQuad->UVPt( 1, i-iFrom ); + sidePoints[ i ].u = uvPt.u; + sidePoints[ i ].v = uvPt.v; + } } //================================================================================ @@ -4717,7 +4959,6 @@ bool FaceQuadStruct::findCell( const gp_XY& UV, int & I, int & J ) y = Min( 1., Max( 0., y )); // precise the position - //int i, j; normPa2IJ( x,y, I,J ); if ( !isNear( UV, I,J )) { @@ -4764,11 +5005,11 @@ void FaceQuadStruct::normPa2IJ(double X, double Y, int & I, int & J ) oldI = I, oldJ = J; while ( X <= UVPt( I,J ).x && I != 0 ) --I; - while ( X > UVPt( I+1,J ).x && I+1 < iSize ) + while ( X > UVPt( I+1,J ).x && I+2 < iSize ) ++I; while ( Y <= UVPt( I,J ).y && J != 0 ) --J; - while ( Y > UVPt( I,J+1 ).y && J+1 < jSize ) + while ( Y > UVPt( I,J+1 ).y && J+2 < jSize ) ++J; } while ( oldI != I || oldJ != J ); } @@ -4798,9 +5039,9 @@ bool FaceQuadStruct::isNear( const gp_XY& UV, int & I, int & J, int nbLoops ) return true; if ( I > 0 && bcI < 0. ) --I; - if ( I+1 < iSize && bcI > 1. ) ++I; + if ( I+2 < iSize && bcI > 1. ) ++I; if ( J > 0 && bcJ < 0. ) --J; - if ( J+1 < jSize && bcJ > 1. ) ++J; + if ( J+2 < jSize && bcJ > 1. ) ++J; uv1 = UVPt( I+1,J+1).UV(); if ( I != oldI || J != oldJ ) @@ -4813,9 +5054,9 @@ bool FaceQuadStruct::isNear( const gp_XY& UV, int & I, int & J, int nbLoops ) return true; if ( I > 0 && bcI > 1. ) --I; - if ( I+1 < iSize && bcI < 0. ) ++I; + if ( I+2 < iSize && bcI < 0. ) ++I; if ( J > 0 && bcJ > 1. ) --J; - if ( J+1 < jSize && bcJ < 0. ) ++J; + if ( J+2 < jSize && bcJ < 0. ) ++J; if ( I == oldI && J == oldJ ) return false; @@ -4840,7 +5081,7 @@ bool FaceQuadStruct::isNear( const gp_XY& UV, int & I, int & J, int nbLoops ) //================================================================================ /*! - * \brief Checks if a given UV is equal to a given frid point + * \brief Checks if a given UV is equal to a given grid point */ //================================================================================ @@ -4987,6 +5228,7 @@ FaceQuadStruct::Side& FaceQuadStruct::Side::operator=(const Side& otherSide) grid = otherSide.grid; from = otherSide.from; to = otherSide.to; + di = otherSide.di; forced_nodes = otherSide.forced_nodes; contacts = otherSide.contacts; nbNodeOut = otherSide.nbNodeOut; @@ -5012,7 +5254,7 @@ FaceQuadStruct::Side& FaceQuadStruct::Side::operator=(const Side& otherSide) int FaceQuadStruct::Side::ToSideIndex( int quadNodeIndex ) const { - return ( from > to ) ? ( from - quadNodeIndex ) : ( quadNodeIndex + from ); + return from + di * quadNodeIndex; } //================================================================================ @@ -5023,7 +5265,31 @@ int FaceQuadStruct::Side::ToSideIndex( int quadNodeIndex ) const int FaceQuadStruct::Side::ToQuadIndex( int sideNodeIndex ) const { - return ( from > to ) ? ( from - sideNodeIndex ) : ( sideNodeIndex - from ); + return ( sideNodeIndex - from ) * di; +} + +//================================================================================ +/*! + * \brief Reverse the side + */ +//================================================================================ + +bool FaceQuadStruct::Side::Reverse(bool keepGrid) +{ + if ( grid ) + { + if ( keepGrid ) + { + from -= di; + to -= di; + std::swap( from, to ); + di *= -1; + } + else + { + grid->Reverse(); + } + } } //================================================================================ @@ -5086,8 +5352,8 @@ void FaceQuadStruct::Side::AddContact( int ip, Side* side, int iop ) double FaceQuadStruct::Side::Param( int i ) const { const vector& points = GetUVPtStruct(); - return (( points[ from + i ].normParam - points[ from ].normParam ) / - ( points[ to - 1 ].normParam - points[ from ].normParam )); + return (( points[ from + i * di ].normParam - points[ from ].normParam ) / + ( points[ to - 1 * di ].normParam - points[ from ].normParam )); } //================================================================================ @@ -5100,7 +5366,7 @@ gp_XY FaceQuadStruct::Side::Value2d( double x ) const { const vector& points = GetUVPtStruct(); double u = ( points[ from ].normParam + - x * ( points[ to-1 ].normParam - points[ from ].normParam )); + x * ( points[ to-di ].normParam - points[ from ].normParam )); return grid->Value2d( u ).XY(); } @@ -5112,8 +5378,19 @@ gp_XY FaceQuadStruct::Side::Value2d( double x ) const double FaceQuadStruct::Side::Length(int theFrom, int theTo) const { + if ( IsReversed() != ( theTo < theFrom )) + std::swap( theTo, theFrom ); + const vector& points = GetUVPtStruct(); - double r = ( points[ Max( to, theTo )-1 ].normParam - - points[ Max( from, theFrom ) ].normParam ); + double r; + if ( theFrom == theTo && theTo == -1 ) + r = Abs( First().normParam - + Last ().normParam ); + else if ( IsReversed() ) + r = Abs( points[ Max( to, theTo+1 ) ].normParam - + points[ Min( from, theFrom ) ].normParam ); + else + r = Abs( points[ Min( to, theTo-1 ) ].normParam - + points[ Max( from, theFrom ) ].normParam ); return r * grid->Length(); }