X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=41cff96229d65cd7f1468843ba25942a3506ffb2;hp=73536f8eb82acd8d75eaa3fe8f0689f663d98659;hb=d4a710ce52f6e76786a7b3845e2f7975dc9a00b1;hpb=e201ae8a2a7d37a84c195803ee97237be92f2f9b diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 73536f8eb..41cff9622 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-2012 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 @@ -32,21 +32,34 @@ #include "SMESHDS_Hypothesis.hxx" #include "SMESH_Algo.hxx" #include "SMESH_ComputeError.hxx" +#include "SMESH_ControlsDef.hxx" #include "SMESH_Gen.hxx" #include "SMESH_Group.hxx" #include "SMESH_Mesh.hxx" #include "SMESH_MesherHelper.hxx" +#include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" #include "SMESH_subMeshEventListener.hxx" -#include "SMESH_ProxyMesh.hxx" #include "utilities.h" +#include #include +#include +#include +#include #include +#include +#include +#include #include +#include #include +#include +#include #include +#include +#include #include #include #include @@ -58,11 +71,11 @@ #include #include #include -#include #include #include -#include +#include +#include //#define __myDEBUG @@ -106,11 +119,13 @@ namespace VISCOUS //-------------------------------------------------------------------------------- /*! * \brief Listener of events of 3D sub-meshes computed with viscous layers. - * It is used to clear an inferior dim sub-mesh modified by viscous layers + * It is used to clear an inferior dim sub-meshes modified by viscous layers */ class _SrinkShapeListener : SMESH_subMeshEventListener { - _SrinkShapeListener(): SMESH_subMeshEventListener(/*isDeletable=*/false) {} + _SrinkShapeListener() + : SMESH_subMeshEventListener(/*isDeletable=*/false, + "StdMeshers_ViscousLayers::_SrinkShapeListener") {} static SMESH_subMeshEventListener* Get() { static _SrinkShapeListener l; return &l; } public: virtual void ProcessEvent(const int event, @@ -150,7 +165,9 @@ namespace VISCOUS */ class _ViscousListener : SMESH_subMeshEventListener { - _ViscousListener(): SMESH_subMeshEventListener(/*isDeletable=*/false) {} + _ViscousListener(): + SMESH_subMeshEventListener(/*isDeletable=*/false, + "StdMeshers_ViscousLayers::_ViscousListener") {} static SMESH_subMeshEventListener* Get() { static _ViscousListener l; return &l; } public: virtual void ProcessEvent(const int event, @@ -214,17 +231,22 @@ namespace VISCOUS - M[0][2]*M[1][1]*M[2][0]); return determinant > 1e-100; } - bool IsForward(const gp_XY& tgtUV, - const TopoDS_Face& face, - SMESH_MesherHelper& helper, - const double refSign) const + bool IsForward(const gp_XY& tgtUV, + const SMDS_MeshNode* smoothedNode, + const TopoDS_Face& face, + SMESH_MesherHelper& helper, + const double refSign) const { - gp_XY prevUV = helper.GetNodeUV( face, _nPrev ); - gp_XY nextUV = helper.GetNodeUV( face, _nNext ); + gp_XY prevUV = helper.GetNodeUV( face, _nPrev, smoothedNode ); + gp_XY nextUV = helper.GetNodeUV( face, _nNext, smoothedNode ); gp_Vec2d v1( tgtUV, prevUV ), v2( tgtUV, nextUV ); double d = v1 ^ v2; return d*refSign > 1e-100; } + bool IsNeighbour(const _Simplex& other) const + { + return _nPrev == other._nNext || _nNext == other._nPrev; + } }; //-------------------------------------------------------------------------------- /*! @@ -267,6 +289,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 +347,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 +387,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 +401,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); }; //-------------------------------------------------------------------------------- /*! @@ -383,6 +423,7 @@ namespace VISCOUS Handle(Geom_Surface)& surface, SMESH_MesherHelper& helper, const double refSign, + bool isCentroidal, bool set3D); }; //-------------------------------------------------------------------------------- @@ -416,7 +457,8 @@ namespace VISCOUS _SolidData& data); void getSimplices( const SMDS_MeshNode* node, vector<_Simplex>& simplices, const set& ingnoreShapes, - const _SolidData* dataToCheckOri = 0); + const _SolidData* dataToCheckOri = 0, + const bool toSort = false); bool sortEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom); void limitStepSize( _SolidData& data, @@ -424,13 +466,20 @@ 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(); bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F, SMESH_MesherHelper& helper, const SMESHDS_SubMesh* faceSubMesh ); + void fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper); bool addBoundaryElements(); bool error( const string& text, int solidID=-1 ); @@ -660,7 +709,7 @@ namespace while ( eIt->more()) { const TopoDS_Edge* e = static_cast( eIt->next() ); - if ( helper.IsSubShape( *e, F ) && BRep_Tool::Curve( *e, loc,f,l)) + if ( helper.IsSubShape( *e, F ) && !BRep_Tool::Curve( *e, loc,f,l).IsNull() ) edges.push_back( *e ); } gp_XYZ dir(0,0,0); @@ -692,6 +741,67 @@ namespace } return dir; } + //================================================================================ + /*! + * \brief Returns true if a FACE is bound by a concave EDGE + */ + //================================================================================ + + bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper ) + { + gp_Vec2d drv1, drv2; + gp_Pnt2d p; + TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE ); + for ( ; eExp.More(); eExp.Next() ) + { + const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() ); + if ( BRep_Tool::Degenerated( E )) continue; + // check if 2D curve is concave + BRepAdaptor_Curve2d curve( E, F ); + const int nbIntervals = curve.NbIntervals( GeomAbs_C2 ); + TColStd_Array1OfReal intervals(1, nbIntervals + 1 ); + curve.Intervals( intervals, GeomAbs_C2 ); + bool isConvex = true; + for ( int i = 1; i <= nbIntervals && isConvex; ++i ) + { + double u1 = intervals( i ); + double u2 = intervals( i+1 ); + curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 ); + double cross = drv2 ^ drv1; + if ( E.Orientation() == TopAbs_REVERSED ) + cross = -cross; + isConvex = ( cross < 1e-9 ); + } + // check if concavity is strong enough to care about it + //const double maxAngle = 5 * Standard_PI180; + if ( !isConvex ) + { + //cout << "Concave FACE " << helper.GetMeshDS()->ShapeToIndex( F ) << endl; + return true; + // map< double, const SMDS_MeshNode* > u2nodes; + // if ( !SMESH_Algo::GetSortedNodesOnEdge( helper.GetMeshDS(), E, + // /*ignoreMedium=*/true, u2nodes)) + // continue; + // map< double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin(); + // gp_Pnt2d uvPrev = helper.GetNodeUV( F, u2n->second ); + // double uPrev = u2n->first; + // for ( ++u2n; u2n != u2nodes.end(); ++u2n ) + // { + // gp_Pnt2d uv = helper.GetNodeUV( F, u2n->second ); + // gp_Vec2d segmentDir( uvPrev, uv ); + // curve.D1( uPrev, p, drv1 ); + // try { + // if ( fabs( segmentDir.Angle( drv1 )) > maxAngle ) + // return true; + // } + // catch ( ... ) {} + // uvPrev = uv; + // uPrev = u2n->first; + // } + } + } + return false; + } //-------------------------------------------------------------------------------- // DEBUG. Dump intermediate node positions into a python script #ifdef __myDEBUG @@ -703,31 +813,38 @@ namespace py = new ofstream(fname); *py << "from smesh import *" << endl << "meshSO = GetCurrentStudy().FindObjectID('0:1:2:3')" << endl - << "mesh = Mesh( meshSO.GetObject()._narrow( SMESH.SMESH_Mesh ))"<SmoothOnEdge(surface, F, helper); + + // try a simple solution on an analytic EDGE + if ( !smoothAnalyticEdge( data, iBeg, iEnd, surface, F, helper )) + { + // smooth on EDGE's + int step = 0; + do { + moved = false; + for ( int i = iBeg; i < iEnd; ++i ) + { + moved |= data._edges[i]->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 +2247,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 +2279,268 @@ 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() ); + dumpMove( tgtNode ); + } + } + 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]); + if ( data._edges[iFrom]->_2neibors->_nodes[0] == + data._edges[iTo-1]->_2neibors->_nodes[1] ) // closed edge + { + int iPeriodic = helper.GetPeriodicIndex(); + if ( iPeriodic == 1 || iPeriodic == 2 ) + { + uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic ))); + if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic )) + std::swap( uv0, uv1 ); + } + } + const gp_XY rangeUV = uv1 - uv0; + for ( int i = iFrom; i < iTo; ++i ) + { + double r = len[i-iFrom] / len.back(); + gp_XY newUV = uv0 + r * rangeUV; + 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() ); + dumpMove( tgtNode ); + + 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.*M_PI; // 0.0 - 2*PI + if ( uMidl < 0 ) uMidl += 2.*M_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() ); + dumpMove( tgtNode ); + + 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 +2612,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 +2623,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 +2639,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 +2652,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 +2899,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 +2920,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 +2950,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 +3033,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 +3055,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 +3101,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 +3122,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 +3135,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 +3163,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() ); @@ -3017,6 +3453,7 @@ bool _ViscousBuilder::shrink() SMDS_ElemIteratorPtr fIt = smDS->GetElements(); while ( fIt->more() ) proxySub->AddElement( fIt->next() ); + // as a result 3D algo will use elements from proxySub and not from smDS } } } @@ -3033,7 +3470,6 @@ bool _ViscousBuilder::shrink() _SolidData& data = *f2sd->second; TNode2Edge& n2eMap = data._n2eMap; const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first )); - const bool reverse = ( data._reversedFaceIds.count( f2sd->first )); Handle(Geom_Surface) surface = BRep_Tool::Surface(F); @@ -3046,8 +3482,8 @@ bool _ViscousBuilder::shrink() // Prepare data for shrinking // =========================== - // Collect nodes to smooth as src nodes are not yet replaced by tgt ones - // and thus all nodes on FACE connected to 2d elements are to be smoothed + // Collect nodes to smooth, as src nodes are not yet replaced by tgt ones + // and thus all nodes on a FACE connected to 2d elements are to be smoothed vector < const SMDS_MeshNode* > smoothNodes; { SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); @@ -3061,12 +3497,15 @@ bool _ViscousBuilder::shrink() // Find out face orientation double refSign = 1; const set ignoreShapes; + bool isOkUV; if ( !smoothNodes.empty() ) { - gp_XY uv = helper.GetNodeUV( F, smoothNodes[0] ); vector<_Simplex> simplices; getSimplices( smoothNodes[0], simplices, ignoreShapes ); - if ( simplices[0].IsForward(uv, F, helper,refSign) != (!reverse)) + helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of silpmex nodes + helper.GetNodeUV( F, simplices[0]._nNext, 0, &isOkUV ); + gp_XY uv = helper.GetNodeUV( F, smoothNodes[0], 0, &isOkUV ); + if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper,refSign) ) refSign = -1; } @@ -3114,6 +3553,9 @@ bool _ViscousBuilder::shrink() } } + // find out if a FACE is concave + const bool isConcaveFace = isConcave( F, helper ); + // Create _SmoothNode's on face F vector< _SmoothNode > nodesToSmooth( smoothNodes.size() ); { @@ -3123,7 +3565,9 @@ bool _ViscousBuilder::shrink() const SMDS_MeshNode* n = smoothNodes[i]; nodesToSmooth[ i ]._node = n; // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices - getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes ); + getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, isConcaveFace ); + // fix up incorrect uv of nodes on the FACE + helper.GetNodeUV( F, n, 0, &isOkUV); dumpMove( n ); } dumpFunctionEnd(); @@ -3166,8 +3610,6 @@ bool _ViscousBuilder::shrink() shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper ); } dumpFunctionEnd(); - if ( !shrinked ) - break; // Move nodes on EDGE's set< _Shrinker1D* >::iterator shr = eShri1D.begin(); @@ -3188,7 +3630,8 @@ bool _ViscousBuilder::shrink() moved = false; for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) { - moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,/*set3D=*/false ); + moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, + /*isCentroidal=*/isConcaveFace,/*set3D=*/false ); } if ( badNb < oldBadNb ) nbNoImpSteps = 0; @@ -3200,18 +3643,38 @@ 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 ); + } + } + if ( !highQuality && isConcaveFace ) + fixBadFaces( F, helper ); // fix narrow faces by swaping diagonals + for ( int st = highQuality ? 10 : 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 ); + nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, + /*isCentroidal=*/isConcaveFace,/*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 + } // loop on FACES to srink mesh on // Replace source nodes by target nodes in shrinked mesh edges @@ -3397,6 +3860,124 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, // return true; } +//================================================================================ +/*! + * \brief Try to fix triangles with high aspect ratio by swaping diagonals + */ +//================================================================================ + +void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper) +{ + SMESH::Controls::AspectRatio qualifier; + SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3); + const double maxAspectRatio = 4.; + + // find bad triangles + + vector< const SMDS_MeshElement* > badTrias; + vector< double > badAspects; + SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F ); + SMDS_ElemIteratorPtr fIt = sm->GetElements(); + while ( fIt->more() ) + { + const SMDS_MeshElement * f = fIt->next(); + if ( f->NbCornerNodes() != 3 ) continue; + for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = SMESH_TNodeXYZ( f->GetNode(iP)); + double aspect = qualifier.GetValue( points ); + if ( aspect > maxAspectRatio ) + { + badTrias.push_back( f ); + badAspects.push_back( aspect ); + } + } + if ( badTrias.empty() ) + return; + + // find couples of faces to swap diagonal + + typedef pair < const SMDS_MeshElement* , const SMDS_MeshElement* > T2Trias; + vector< T2Trias > triaCouples; + + TIDSortedElemSet involvedFaces, emptySet; + for ( size_t iTia = 0; iTia < badTrias.size(); ++iTia ) + { + T2Trias trias [3]; + double aspRatio [3]; + int i1, i2, i3; + + involvedFaces.insert( badTrias[iTia] ); + for ( int iP = 0; iP < 3; ++iP ) + points(iP+1) = SMESH_TNodeXYZ( badTrias[iTia]->GetNode(iP)); + + // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping + int bestCouple = -1; + for ( int iSide = 0; iSide < 3; ++iSide ) + { + const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide ); + const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 ); + trias [iSide].first = badTrias[iTia]; + trias [iSide].second = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, involvedFaces, + & i1, & i2 ); + if ( ! trias[iSide].second || trias[iSide].second->NbCornerNodes() != 3 ) + continue; + + // aspect ratio of an adjacent tria + for ( int iP = 0; iP < 3; ++iP ) + points2(iP+1) = SMESH_TNodeXYZ( trias[iSide].second->GetNode(iP)); + double aspectInit = qualifier.GetValue( points2 ); + + // arrange nodes as after diag-swaping + if ( helper.WrapIndex( i1+1, 3 ) == i2 ) + i3 = helper.WrapIndex( i1-1, 3 ); + else + i3 = helper.WrapIndex( i1+1, 3 ); + points1 = points; + points1( 1+ iSide ) = points2( 1+ i3 ); + points2( 1+ i2 ) = points1( 1+ ( iSide+2 ) % 3 ); + + // aspect ratio after diag-swaping + aspRatio[ iSide ] = qualifier.GetValue( points1 ) + qualifier.GetValue( points2 ); + if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] ) + continue; + + if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] ) + bestCouple = iSide; + } + + if ( bestCouple >= 0 ) + { + triaCouples.push_back( trias[bestCouple] ); + involvedFaces.insert ( trias[bestCouple].second ); + } + else + { + involvedFaces.erase( badTrias[iTia] ); + } + } + if ( triaCouples.empty() ) + return; + + // swap diagonals + + SMESH_MeshEditor editor( helper.GetMesh() ); + dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")< uv( _simplices.size() ); + for ( size_t i = 0; i < _simplices.size(); ++i ) + uv[i] = helper.GetNodeUV( face, _simplices[i]._nPrev, _node ); + // compute new UV for the node gp_XY newPos (0,0); - for ( unsigned i = 0; i < _simplices.size(); ++i ) - newPos += helper.GetNodeUV( face, _simplices[i]._nPrev ); - newPos /= _simplices.size(); + if ( isCentroidal && _simplices.size() > 3 ) + { + // average centers of diagonals wieghted with their reciprocal lengths + if ( _simplices.size() == 4 ) + { + double w1 = 1. / ( uv[2]-uv[0] ).SquareModulus(); + double w2 = 1. / ( uv[3]-uv[1] ).SquareModulus(); + newPos = ( w1 * ( uv[2]+uv[0] ) + w2 * ( uv[3]+uv[1] )) / ( w1+w2 ) / 2; + } + else + { + double sumWeight = 0; + int nb = _simplices.size() == 4 ? 2 : _simplices.size(); + for ( int i = 0; i < nb; ++i ) + { + int iFrom = i + 2; + int iTo = i + _simplices.size() - 1; + for ( int j = iFrom; j < iTo; ++j ) + { + int i2 = SMESH_MesherHelper::WrapIndex( j, _simplices.size() ); + double w = 1. / ( uv[i]-uv[i2] ).SquareModulus(); + sumWeight += w; + newPos += w * ( uv[i]+uv[i2] ); + } + } + newPos /= 2 * sumWeight; + } + } + else + { + // Laplacian smooth + isCentroidal = false; + for ( size_t i = 0; i < _simplices.size(); ++i ) + newPos += uv[i]; + newPos /= _simplices.size(); + } // count quality metrics (orientation) of triangles around the node int nbOkBefore = 0; gp_XY tgtUV = helper.GetNodeUV( face, _node ); for ( unsigned i = 0; i < _simplices.size(); ++i ) - nbOkBefore += _simplices[i].IsForward( tgtUV, face, helper, refSign ); + nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign ); int nbOkAfter = 0; for ( unsigned i = 0; i < _simplices.size(); ++i ) - nbOkAfter += _simplices[i].IsForward( newPos, face, helper, refSign ); + nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign ); if ( nbOkAfter < nbOkBefore ) + { + // if ( isCentroidal ) + // return Smooth( badNb, surface, helper, refSign, !isCentroidal, set3D ); + badNb += _simplices.size() - nbOkBefore; return false; + } SMDS_FacePosition* pos = static_cast( _node->GetPosition() ); pos->SetUParameter( newPos.X() ); @@ -3597,7 +4222,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 +4276,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] ); @@ -3708,6 +4333,7 @@ void _Shrinker1D::RestoreParams() } _done = false; } + //================================================================================ /*! * \brief Replace source nodes by target nodes in shrinked mesh edges @@ -3801,12 +4427,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 +4450,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 +4463,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;