X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=73ceef8747f306cb2ec491517a68a1cde760dad0;hb=7eb6cc5064939f87836d9c55161514f365c16ec4;hp=73536f8eb82acd8d75eaa3fe8f0689f663d98659;hpb=e201ae8a2a7d37a84c195803ee97237be92f2f9b;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 73536f8eb..73ceef874 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -1,20 +1,20 @@ -// Copyright (C) 2007-2011 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2011 CEA/DEN, EDF R&D, OPEN CASCADE // -// This library is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 2.1 of the License. +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. // -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. // -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // -// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // // File : StdMeshers_ViscousLayers.cxx @@ -43,9 +43,18 @@ #include "utilities.h" #include +#include +#include +#include #include +#include +#include +#include #include +#include #include +#include +#include #include #include #include @@ -63,6 +72,7 @@ #include #include #include +#include //#define __myDEBUG @@ -267,6 +277,10 @@ namespace VISCOUS gp_XYZ* _plnNorm; _2NearEdges() { _nodes[0]=_nodes[1]=0; _plnNorm = 0; } + void reverse() { + std::swap( _nodes[0], _nodes[1] ); + std::swap( _wgt[0], _wgt[1] ); + } }; //-------------------------------------------------------------------------------- /*! @@ -321,6 +335,14 @@ namespace VISCOUS void Copy( _LayerEdge& other, SMESH_MesherHelper& helper ); void SetCosin( double cosin ); }; + struct _LayerEdgeCmp + { + bool operator () (const _LayerEdge* e1, const _LayerEdge* e2) const + { + const bool cmpNodes = ( e1 && e2 && e1->_nodes.size() && e2->_nodes.size() ); + return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 ); + } + }; //-------------------------------------------------------------------------------- typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge; @@ -353,9 +375,8 @@ namespace VISCOUS // FACE's WOL, srink on which is forbiden due to algo on the adjacent SOLID set< TGeomID > _noShrinkFaces; - // end index in _edges of _LayerEdge's based on EDGE (map key) to - // FACE (maybe NULL) they are inflated along - //map< int, TopoDS_Face > _endEdge2Face; + // to + map< TGeomID,Handle(Geom_Curve)> _edge2curve; // end indices in _edges of _LayerEdge on one shape to smooth vector< int > _endEdgeToSmooth; @@ -368,6 +389,13 @@ namespace VISCOUS const StdMeshers_ViscousLayers* h=0, _MeshOfSolid* m=0) :_solid(s), _hyp(h), _proxyMesh(m) {} ~_SolidData(); + + Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge& E, + const int iFrom, + const int iTo, + Handle(Geom_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper); }; //-------------------------------------------------------------------------------- /*! @@ -424,7 +452,13 @@ namespace VISCOUS const double cosin); void limitStepSize( _SolidData& data, const double minSize); bool inflate(_SolidData& data); - bool smoothAndCheck(_SolidData& data, int nbSteps, double & distToIntersection); + bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection); + bool smoothAnalyticEdge( _SolidData& data, + const int iFrom, + const int iTo, + Handle(Geom_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper); bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper ); bool refine(_SolidData& data); bool shrink(); @@ -1279,7 +1313,7 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE || curN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) { - double dist = SMESH_MeshEditor::TNodeXYZ( face->GetNode(i)).Distance( nextN ); + double dist = SMESH_TNodeXYZ( face->GetNode(i)).Distance( nextN ); if ( dist < minSize ) minSize = dist, iN = i; } @@ -1308,7 +1342,7 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize) if ( data._stepSizeNodes[0] ) { double dist = - SMESH_MeshEditor::TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]); + SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]); data._stepSizeCoeff = data._stepSize / dist; } } @@ -1329,7 +1363,7 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, list< TGeomID > shapesToSmooth; SMESH_MesherHelper helper( *_mesh ); - bool ok; + bool ok = true; for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS ) { @@ -1595,7 +1629,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, } else { - edge._pos.push_back( SMESH_MeshEditor::TNodeXYZ( node )); + edge._pos.push_back( SMESH_TNodeXYZ( node )); if ( posType == SMDS_TOP_FACE ) { @@ -1603,7 +1637,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, double avgNormProj = 0, avgLen = 0; for ( unsigned i = 0; i < edge._simplices.size(); ++i ) { - gp_XYZ vec = edge._pos.back() - SMESH_MeshEditor::TNodeXYZ( edge._simplices[i]._nPrev ); + gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev ); avgNormProj += edge._normal * vec; avgLen += vec.Modulus(); } @@ -1693,9 +1727,9 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, if ( _nodes[0]->GetPosition()->GetTypeOfPosition() != SMDS_TOP_EDGE ) return; - gp_XYZ pos = SMESH_MeshEditor::TNodeXYZ( _nodes[0] ); - gp_XYZ vec1 = pos - SMESH_MeshEditor::TNodeXYZ( n1 ); - gp_XYZ vec2 = pos - SMESH_MeshEditor::TNodeXYZ( n2 ); + gp_XYZ pos = SMESH_TNodeXYZ( _nodes[0] ); + gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 ); + gp_XYZ vec2 = pos - SMESH_TNodeXYZ( n2 ); // Set _curvature @@ -1836,7 +1870,7 @@ void _ViscousBuilder::makeGroupOfLE() for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j ) { _LayerEdge& edge = *_sdVec[i]._edges[j]; - SMESH_MeshEditor::TNodeXYZ nXYZ( edge._nodes[0] ); + SMESH_TNodeXYZ nXYZ( edge._nodes[0] ); nXYZ += edge._normal * _sdVec[i]._stepSize; dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <GetID() << ", mesh.AddNode( " << nXYZ.X()<<","<< nXYZ.Y()<<","<< nXYZ.Z()<<")])"); @@ -1969,7 +2003,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) limitStepSize( data, 0.25 * distToIntersection ); if ( data._stepSizeNodes[0] ) data._stepSize = data._stepSizeCoeff * - SMESH_MeshEditor::TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]); + SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]); } if (nbSteps == 0 ) @@ -1985,7 +2019,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) //================================================================================ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, - int nbSteps, + const int nbSteps, double & distToIntersection) { if ( data._endEdgeToSmooth.empty() ) @@ -2019,21 +2053,25 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, TGeomID sInd = data._edges[ iBeg ]->_nodes[0]->getshapeId(); if ( data._edges[ iBeg ]->IsOnEdge() ) - { - dumpFunction(SMESH_Comment("smooth")<SmoothOnEdge(surface, F, helper); + { // try a simple solution on an analytic EDGE + if ( !smoothAnalyticEdge( data, iBeg, iEnd, surface, F, helper )) + { + dumpFunction(SMESH_Comment("smooth")<SmoothOnEdge(surface, F, helper); + } + dumpCmd( SMESH_Comment("# end step ")<_nodes.back() ); + SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() ); for ( unsigned j = 0; j < edge->_simplices.size(); ++j ) if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ )) { @@ -2083,14 +2121,23 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, distToIntersection = Precision::Infinite(); double dist; - const SMDS_MeshElement* intFace = 0, *closestFace = 0; + const SMDS_MeshElement* intFace = 0; +#ifdef __myDEBUG + const SMDS_MeshElement* closestFace = 0; int iLE = 0; +#endif for ( unsigned i = 0; i < data._edges.size(); ++i ) { if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace )) return false; if ( distToIntersection > dist ) - distToIntersection = dist, closestFace = intFace, iLE = i; + { + distToIntersection = dist; +#ifdef __myDEBUG + iLE = i; + closestFace = intFace; +#endif + } } #ifdef __myDEBUG if ( closestFace ) @@ -2106,6 +2153,253 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, return true; } +//================================================================================ +/*! + * \brief Return a curve of the EDGE to be used for smoothing and arrange + * _LayerEdge's to be in a consequent order + */ +//================================================================================ + +Handle(Geom_Curve) _SolidData::CurveForSmooth( const TopoDS_Edge& E, + const int iFrom, + const int iTo, + Handle(Geom_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper) +{ + TGeomID eIndex = helper.GetMeshDS()->ShapeToIndex( E ); + + map< TGeomID, Handle(Geom_Curve)>::iterator i2curve = _edge2curve.find( eIndex ); + + if ( i2curve == _edge2curve.end() ) + { + // sort _LayerEdge's by position on the EDGE + { + map< double, _LayerEdge* > u2edge; + for ( int i = iFrom; i < iTo; ++i ) + u2edge.insert( make_pair( helper.GetNodeU( E, _edges[i]->_nodes[0] ), _edges[i] )); + + ASSERT( u2edge.size() == iTo - iFrom ); + map< double, _LayerEdge* >::iterator u2e = u2edge.begin(); + for ( int i = iFrom; i < iTo; ++i, ++u2e ) + _edges[i] = u2e->second; + + // set _2neibors according to the new order + for ( int i = iFrom; i < iTo-1; ++i ) + if ( _edges[i]->_2neibors->_nodes[1] != _edges[i+1]->_nodes.back() ) + _edges[i]->_2neibors->reverse(); + if ( u2edge.size() > 1 && + _edges[iTo-1]->_2neibors->_nodes[0] != _edges[iTo-2]->_nodes.back() ) + _edges[iTo-1]->_2neibors->reverse(); + } + + SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( eIndex ); + + TopLoc_Location loc; double f,l; + + Handle(Geom_Line) line; + Handle(Geom_Circle) circle; + bool isLine, isCirc; + if ( F.IsNull() ) // 3D case + { + // check if the EDGE is a line + Handle(Geom_Curve) curve = BRep_Tool::Curve( E, loc, f, l); + if ( curve->IsKind( STANDARD_TYPE( Geom_TrimmedCurve ))) + curve = Handle(Geom_TrimmedCurve)::DownCast( curve )->BasisCurve(); + + line = Handle(Geom_Line)::DownCast( curve ); + circle = Handle(Geom_Circle)::DownCast( curve ); + isLine = (!line.IsNull()); + isCirc = (!circle.IsNull()); + + if ( !isLine && !isCirc ) // Check if the EDGE is close to a line + { + Bnd_B3d bndBox; + SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); + while ( nIt->more() ) + bndBox.Add( SMESH_TNodeXYZ( nIt->next() )); + gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin(); + + SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->_nodes[0] ); + SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->_nodes[1] ); + const double lineTol = 1e-2 * ( p0 - p1 ).Modulus(); + for ( int i = 0; i < 3 && !isLine; ++i ) + isLine = ( size.Coord( i+1 ) <= lineTol ); + } + if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle + { + // TODO + } + } + else // 2D case + { + // check if the EDGE is a line + Handle(Geom2d_Curve) curve = BRep_Tool::CurveOnSurface( E, F, f, l); + if ( curve->IsKind( STANDARD_TYPE( Geom2d_TrimmedCurve ))) + curve = Handle(Geom2d_TrimmedCurve)::DownCast( curve )->BasisCurve(); + + Handle(Geom2d_Line) line2d = Handle(Geom2d_Line)::DownCast( curve ); + Handle(Geom2d_Circle) circle2d = Handle(Geom2d_Circle)::DownCast( curve ); + isLine = (!line2d.IsNull()); + isCirc = (!circle2d.IsNull()); + + if ( !isLine && !isCirc) // Check if the EDGE is close to a line + { + Bnd_B2d bndBox; + SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); + while ( nIt->more() ) + bndBox.Add( helper.GetNodeUV( F, nIt->next() )); + gp_XY size = bndBox.CornerMax() - bndBox.CornerMin(); + + const double lineTol = 1e-2 * sqrt( bndBox.SquareExtent() ); + for ( int i = 0; i < 2 && !isLine; ++i ) + isLine = ( size.Coord( i+1 ) <= lineTol ); + } + if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle + { + // TODO + } + if ( isLine ) + { + line = new Geom_Line( gp::OX() ); // only type does matter + } + else if ( isCirc ) + { + gp_Pnt2d p = circle2d->Location(); + gp_Ax2 ax( gp_Pnt( p.X(), p.Y(), 0), gp::DX()); + circle = new Geom_Circle( ax, 1.); // only center position does matter + } + } + + Handle(Geom_Curve)& res = _edge2curve[ eIndex ]; + if ( isLine ) + res = line; + else if ( isCirc ) + res = circle; + + return res; + } + return i2curve->second; +} + +//================================================================================ +/*! + * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE + */ +//================================================================================ + +bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data, + const int iFrom, + const int iTo, + Handle(Geom_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper) +{ + TopoDS_Shape S = helper.GetSubShapeByNode( data._edges[ iFrom ]->_nodes[0], + helper.GetMeshDS()); + TopoDS_Edge E = TopoDS::Edge( S ); + + Handle(Geom_Curve) curve = data.CurveForSmooth( E, iFrom, iTo, surface, F, helper ); + if ( curve.IsNull() ) return false; + + // compute a relative length of segments + vector< double > len( iTo-iFrom+1 ); + { + double curLen, prevLen = len[0] = 1.0; + for ( int i = iFrom; i < iTo; ++i ) + { + curLen = prevLen * data._edges[i]->_2neibors->_wgt[0] / data._edges[i]->_2neibors->_wgt[1]; + len[i-iFrom+1] = len[i-iFrom] + curLen; + prevLen = curLen; + } + } + + if ( curve->IsKind( STANDARD_TYPE( Geom_Line ))) + { + if ( F.IsNull() ) // 3D + { + SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->_nodes[0]); + SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->_nodes[1]); + for ( int i = iFrom; i < iTo; ++i ) + { + double r = len[i-iFrom] / len.back(); + gp_XYZ newPos = p0 * ( 1. - r ) + p1 * r; + data._edges[i]->_pos.back() = newPos; + SMDS_MeshNode* tgtNode = const_cast( data._edges[i]->_nodes.back() ); + tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); + } + } + else + { + gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]); + gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]); + for ( int i = iFrom; i < iTo; ++i ) + { + double r = len[i-iFrom] / len.back(); + gp_XY newUV = uv0 * ( 1. - r ) + uv1 * r; + data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); + + gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() ); + SMDS_MeshNode* tgtNode = const_cast( data._edges[i]->_nodes.back() ); + tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); + + SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); + pos->SetUParameter( newUV.X() ); + pos->SetVParameter( newUV.Y() ); + } + } + return true; + } + + if ( curve->IsKind( STANDARD_TYPE( Geom_Circle ))) + { + Handle(Geom_Circle) circle = Handle(Geom_Circle)::DownCast( curve ); + gp_Pnt center3D = circle->Location(); + + if ( F.IsNull() ) // 3D + { + return false; // TODO ??? + } + else // 2D + { + const gp_XY center( center3D.X(), center3D.Y() ); + + gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]); + gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back()); + gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]); + gp_Vec2d vec0( center, uv0 ); + gp_Vec2d vecM( center, uvM); + gp_Vec2d vec1( center, uv1 ); + double uLast = vec0.Angle( vec1 ); // -PI - +PI + double uMidl = vec0.Angle( vecM ); + if ( uLast < 0 ) uLast += 2*PI; // 0.0 - 2*PI + if ( uMidl < 0 ) uMidl += 2*PI; + const bool sense = ( uMidl < uLast ); + const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() ); + + gp_Ax2d axis( center, vec0 ); + gp_Circ2d circ ( axis, radius, sense ); + for ( int i = iFrom; i < iTo; ++i ) + { + double newU = uLast * len[i-iFrom] / len.back(); + gp_Pnt2d newUV = ElCLib::Value( newU, circ ); + data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); + + gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() ); + SMDS_MeshNode* tgtNode = const_cast( data._edges[i]->_nodes.back() ); + tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); + + SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); + pos->SetUParameter( newUV.X() ); + pos->SetVParameter( newUV.Y() ); + } + } + return true; + } + + return false; +} + //================================================================================ /*! * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with @@ -2177,7 +2471,8 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, // 1) Find intersections double dist; const SMDS_MeshElement* face; - map< _LayerEdge*, set< _LayerEdge* > > edge2CloseEdge; + typedef map< _LayerEdge*, set< _LayerEdge*, _LayerEdgeCmp >, _LayerEdgeCmp > TLEdge2LEdgeSet; + TLEdge2LEdgeSet edge2CloseEdge; const double eps = data._epsilon * data._epsilon; for ( unsigned i = 0; i < data._edges.size(); ++i ) @@ -2187,7 +2482,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, if ( edge->FindIntersection( *searcher, dist, eps, &face )) { const TmpMeshFaceOnEdge* f = (const TmpMeshFaceOnEdge*) face; - set< _LayerEdge* > & ee = edge2CloseEdge[ edge ]; + set< _LayerEdge*, _LayerEdgeCmp > & ee = edge2CloseEdge[ edge ]; ee.insert( f->_le1 ); ee.insert( f->_le2 ); if ( f->_le1->IsOnEdge() && f->_le1->_sWOL.IsNull() ) @@ -2203,12 +2498,12 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, { dumpFunction(SMESH_Comment("updateNormals")< >::iterator e2ee = edge2CloseEdge.begin(); + TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin(); for ( ; e2ee != edge2CloseEdge.end(); ++e2ee ) { _LayerEdge* edge1 = e2ee->first; _LayerEdge* edge2 = 0; - set< _LayerEdge* >& ee = e2ee->second; + set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second; // find EDGEs the edges reside TopoDS_Edge E1, E2; @@ -2216,7 +2511,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, if ( S.ShapeType() != TopAbs_EDGE ) continue; // TODO: find EDGE by VERTEX E1 = TopoDS::Edge( S ); - set< _LayerEdge* >::iterator eIt = ee.begin(); + set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin(); while ( E2.IsNull() && eIt != ee.end()) { _LayerEdge* e2 = *eIt++; @@ -2463,7 +2758,7 @@ gp_Ax1 _LayerEdge::LastSegment(double& segLen) const gp_Ax1 segDir; if ( iPrev < 0 ) { - segDir.SetLocation( SMESH_MeshEditor::TNodeXYZ( _nodes[0] )); + segDir.SetLocation( SMESH_TNodeXYZ( _nodes[0] )); segDir.SetDirection( _normal ); segLen = 0; } @@ -2484,7 +2779,7 @@ gp_Ax1 _LayerEdge::LastSegment(double& segLen) const Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc ); pPrev = surface->Value( pPrev.X(), pPrev.Y() ).Transformed( loc ); } - dir = SMESH_MeshEditor::TNodeXYZ( _nodes.back() ) - pPrev.XYZ(); + dir = SMESH_TNodeXYZ( _nodes.back() ) - pPrev.XYZ(); } segDir.SetLocation( pPrev ); segDir.SetDirection( dir ); @@ -2514,9 +2809,9 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment, gp_XYZ orig = lastSegment.Location().XYZ(); gp_XYZ dir = lastSegment.Direction().XYZ(); - SMESH_MeshEditor::TNodeXYZ vert0( n0 ); - SMESH_MeshEditor::TNodeXYZ vert1( n1 ); - SMESH_MeshEditor::TNodeXYZ vert2( n2 ); + SMESH_TNodeXYZ vert0( n0 ); + SMESH_TNodeXYZ vert1( n1 ); + SMESH_TNodeXYZ vert2( n2 ); /* calculate distance from vert0 to ray origin */ gp_XYZ tvec = orig - vert0; @@ -2597,11 +2892,11 @@ bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface, ASSERT( IsOnEdge() ); SMDS_MeshNode* tgtNode = const_cast( _nodes.back() ); - SMESH_MeshEditor::TNodeXYZ oldPos( tgtNode ); + SMESH_TNodeXYZ oldPos( tgtNode ); double dist01, distNewOld; - SMESH_MeshEditor::TNodeXYZ p0( _2neibors->_nodes[0]); - SMESH_MeshEditor::TNodeXYZ p1( _2neibors->_nodes[1]); + SMESH_TNodeXYZ p0( _2neibors->_nodes[0]); + SMESH_TNodeXYZ p1( _2neibors->_nodes[1]); dist01 = p0.Distance( _2neibors->_nodes[1] ); gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1]; @@ -2619,7 +2914,7 @@ bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface, if ( _2neibors->_plnNorm ) { // put newPos on the plane defined by source node and _plnNorm - gp_XYZ new2src = SMESH_MeshEditor::TNodeXYZ( _nodes[0] ) - newPos.XYZ(); + gp_XYZ new2src = SMESH_TNodeXYZ( _nodes[0] ) - newPos.XYZ(); double new2srcProj = (*_2neibors->_plnNorm) * new2src; newPos.ChangeCoord() += (*_2neibors->_plnNorm) * new2srcProj; } @@ -2665,7 +2960,7 @@ bool _LayerEdge::Smooth(int& badNb) // compute new position for the last _pos gp_XYZ newPos (0,0,0); for ( unsigned i = 0; i < _simplices.size(); ++i ) - newPos += SMESH_MeshEditor::TNodeXYZ( _simplices[i]._nPrev ); + newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev ); newPos /= _simplices.size(); if ( _curvature ) @@ -2686,7 +2981,7 @@ bool _LayerEdge::Smooth(int& badNb) // } // count quality metrics (orientation) of tetras around _tgtNode int nbOkBefore = 0; - SMESH_MeshEditor::TNodeXYZ tgtXYZ( _nodes.back() ); + SMESH_TNodeXYZ tgtXYZ( _nodes.back() ); for ( unsigned i = 0; i < _simplices.size(); ++i ) nbOkBefore += _simplices[i].IsForward( _nodes[0], &tgtXYZ ); @@ -2699,7 +2994,7 @@ bool _LayerEdge::Smooth(int& badNb) SMDS_MeshNode* n = const_cast< SMDS_MeshNode* >( _nodes.back() ); - _len -= prevPos.Distance(SMESH_MeshEditor::TNodeXYZ( n )); + _len -= prevPos.Distance(SMESH_TNodeXYZ( n )); _len += prevPos.Distance(newPos); n->setXYZ( newPos.X(), newPos.Y(), newPos.Z()); @@ -2727,7 +3022,7 @@ void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper ) } SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() ); - SMESH_MeshEditor::TNodeXYZ oldXYZ( n ); + SMESH_TNodeXYZ oldXYZ( n ); gp_XYZ nXYZ = oldXYZ + _normal * ( len - _len ) * _lenFactor; n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() ); @@ -3200,15 +3495,32 @@ bool _ViscousBuilder::shrink() if ( badNb > 0 ) return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first ); } - // No wrongly shaped faces remain; final smooth. Set node XYZ - for ( int st = 3; st; --st ) + // No wrongly shaped faces remain; final smooth. Set node XYZ. + // First, find out a needed quality of smoothing (high for quadrangles only) + bool highQuality; + { + const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles(); + if ( hasTria != hasQuad ) + { + highQuality = hasQuad; + } + else + { + set nbNodesSet; + SMDS_ElemIteratorPtr fIt = smDS->GetElements(); + while ( fIt->more() && nbNodesSet.size() < 2 ) + nbNodesSet.insert( fIt->next()->NbCornerNodes() ); + highQuality = ( *nbNodesSet.begin() == 4 ); + } + } + for ( int st = highQuality ? 8 : 3; st; --st ) { dumpFunction(SMESH_Comment("shrinkFace")<first<<"_st"<<++smooStep); // debug for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,/*set3D=*/st==1 ); dumpFunctionEnd(); } - // Set event listener to clear FACE sub-mesh together with SOLID sub-mesh + // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh _SrinkShapeListener::ToClearSubMeshWithSolid( sm, data._solid ); }// loop on FACES to srink mesh on @@ -3597,7 +3909,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper ) return; TopLoc_Location loc; Handle(Geom_Curve) C = BRep_Tool::Curve(E, loc, f,l); - GeomAdaptor_Curve aCurve(C); + GeomAdaptor_Curve aCurve(C, f,l); const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l); int nbExpectNodes = eSubMesh->NbNodes() - e->_nodes.size(); @@ -3651,7 +3963,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper) if ( set3D || _done ) { Handle(Geom_Curve) C = BRep_Tool::Curve(E, f,l); - GeomAdaptor_Curve aCurve(C); + GeomAdaptor_Curve aCurve(C, f,l); if ( _edges[0] ) f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] ); @@ -3801,12 +4113,12 @@ bool _ViscousBuilder::addBoundaryElements() // Find out orientation and type of face to create - bool reverse = false, tria = false, isOnFace; + bool reverse = false, isOnFace; map< TGeomID, TopoDS_Shape >::iterator e2f = data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E )); TopoDS_Shape F; - if ( isOnFace = ( e2f != data._shrinkShape2Shape.end() )) + if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() ))) { F = e2f->second.Oriented( TopAbs_FORWARD ); reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED ); @@ -3824,8 +4136,6 @@ bool _ViscousBuilder::addBoundaryElements() !_ignoreShapeIds.count( e2f->first )) F = *pF; } - - tria = true; } // Find the sub-mesh to add new faces SMESHDS_SubMesh* sm = 0; @@ -3839,8 +4149,6 @@ bool _ViscousBuilder::addBoundaryElements() // Make faces const int dj1 = reverse ? 0 : 1; const int dj2 = reverse ? 1 : 0; - vector newFaces; - newFaces.reserve(( ledges.size() - 1 ) * (ledges[0]->_nodes.size() - 1 )); for ( unsigned j = 1; j < ledges.size(); ++j ) { vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes;