X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=124b61e1e66e93fdabe4468d1de23cc8247b1e2c;hp=9487a6148b003600eb140c0fdc2178a5b975f21e;hb=4cd2499bddcd3da3ec8900fe825bc98669b789b5;hpb=b6fe67a04f6b63b1d1b09c333ca140fa2be79ff8 diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 9487a6148..124b61e1e 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -1,4 +1,4 @@ -// 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 @@ -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,10 @@ #include #include #include -#include #include #include -#include +#include #include //#define __myDEBUG @@ -70,7 +82,7 @@ using namespace std; //================================================================================ -namespace VISCOUS +namespace VISCOUS_3D { typedef int TGeomID; @@ -107,13 +119,15 @@ 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 + class _ShrinkShapeListener : SMESH_subMeshEventListener { - _SrinkShapeListener(): SMESH_subMeshEventListener(/*isDeletable=*/false) {} - static SMESH_subMeshEventListener* Get() { static _SrinkShapeListener l; return &l; } + _ShrinkShapeListener() + : SMESH_subMeshEventListener(/*isDeletable=*/false, + "StdMeshers_ViscousLayers::_ShrinkShapeListener") {} public: + static SMESH_subMeshEventListener* Get() { static _ShrinkShapeListener l; return &l; } virtual void ProcessEvent(const int event, const int eventType, SMESH_subMesh* solidSM, @@ -125,23 +139,6 @@ namespace VISCOUS SMESH_subMeshEventListener::ProcessEvent(event,eventType,solidSM,data,hyp); } } - static void ToClearSubMeshWithSolid( SMESH_subMesh* sm, - const TopoDS_Shape& solid) - { - SMESH_subMesh* solidSM = sm->GetFather()->GetSubMesh( solid ); - SMESH_subMeshEventListenerData* data = solidSM->GetEventListenerData( Get()); - if ( data ) - { - if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sm ) == - data->mySubMeshes.end()) - data->mySubMeshes.push_back( sm ); - } - else - { - data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sm ); - sm->SetEventListener( Get(), data, /*whereToListenTo=*/solidSM ); - } - } }; //-------------------------------------------------------------------------------- /*! @@ -151,7 +148,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, @@ -189,6 +188,32 @@ namespace VISCOUS } }; + //================================================================================ + /*! + * \brief sets a sub-mesh event listener to clear sub-meshes of sub-shapes of + * the main shape when sub-mesh of the main shape is cleared, + * for example to clear sub-meshes of FACEs when sub-mesh of a SOLID + * is cleared + */ + //================================================================================ + + void ToClearSubWithMain( SMESH_subMesh* sub, const TopoDS_Shape& main) + { + SMESH_subMesh* mainSM = sub->GetFather()->GetSubMesh( main ); + SMESH_subMeshEventListenerData* data = + mainSM->GetEventListenerData( _ShrinkShapeListener::Get()); + if ( data ) + { + if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sub ) == + data->mySubMeshes.end()) + data->mySubMeshes.push_back( sub ); + } + else + { + data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sub ); + sub->SetEventListener( _ShrinkShapeListener::Get(), data, /*whereToListenTo=*/mainSM ); + } + } //-------------------------------------------------------------------------------- /*! * \brief Simplex (triangle or tetrahedron) based on 1 (tria) or 2 (tet) nodes of @@ -215,17 +240,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; + } }; //-------------------------------------------------------------------------------- /*! @@ -235,6 +265,7 @@ namespace VISCOUS { double _r; // radius double _k; // factor to correct node smoothed position + double _h2lenRatio; // avgNormProj / (2*avgDist) public: static _Curvature* New( double avgNormProj, double avgDist ) { @@ -245,10 +276,12 @@ namespace VISCOUS c->_r = avgDist * avgDist / avgNormProj; c->_k = avgDist * avgDist / c->_r / c->_r; c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive + c->_h2lenRatio = avgNormProj / ( avgDist + avgDist ); } return c; } double lenDelta(double len) const { return _k * ( _r + len ); } + double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; } }; struct _LayerEdge; //-------------------------------------------------------------------------------- @@ -268,6 +301,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] ); + } }; //-------------------------------------------------------------------------------- /*! @@ -362,9 +399,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; @@ -377,6 +413,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); }; //-------------------------------------------------------------------------------- /*! @@ -392,6 +435,7 @@ namespace VISCOUS Handle(Geom_Surface)& surface, SMESH_MesherHelper& helper, const double refSign, + bool isCentroidal, bool set3D); }; //-------------------------------------------------------------------------------- @@ -425,7 +469,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, @@ -433,13 +478,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 ); @@ -487,7 +539,8 @@ namespace VISCOUS virtual SMDSAbs_ElementType GetType() const { return SMDSAbs_Face; } virtual vtkIdType GetVtkType() const { return -1; } virtual SMDSAbs_EntityType GetEntityType() const { return SMDSEntity_Last; } - virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const + virtual SMDSAbs_GeometryType GetGeomType() const { return SMDSGeom_TRIANGLE; } +virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));} }; //-------------------------------------------------------------------------------- @@ -506,7 +559,7 @@ namespace VISCOUS _nn[3]=_le2->_nodes[0]; } }; -} // namespace VISCOUS +} // namespace VISCOUS_3D //================================================================================ // StdMeshers_ViscousLayers hypothesis @@ -518,10 +571,15 @@ StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH _name = StdMeshers_ViscousLayers::GetHypType(); _param_algo_dim = -3; // auxiliary hyp used by 3D algos } // -------------------------------------------------------------------------------- -void StdMeshers_ViscousLayers::SetIgnoreFaces(const std::vector& faceIds) +void StdMeshers_ViscousLayers::SetBndShapesToIgnore(const std::vector& faceIds) +{ + if ( faceIds != _ignoreBndShapeIds ) + _ignoreBndShapeIds = faceIds, NotifySubMeshesHypothesisModification(); +} // -------------------------------------------------------------------------------- +bool StdMeshers_ViscousLayers::IsIgnoredShape(const int shapeID) const { - if ( faceIds != _ignoreFaceIds ) - _ignoreFaceIds = faceIds, NotifySubMeshesHypothesisModification(); + return ( find( _ignoreBndShapeIds.begin(), _ignoreBndShapeIds.end(), shapeID ) + != _ignoreBndShapeIds.end() ); } // -------------------------------------------------------------------------------- void StdMeshers_ViscousLayers::SetTotalThickness(double thickness) { @@ -543,7 +601,7 @@ StdMeshers_ViscousLayers::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theShape, const bool toMakeN2NMap) const { - using namespace VISCOUS; + using namespace VISCOUS_3D; _ViscousBuilder bulder; SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape ); if ( err && !err->IsOK() ) @@ -579,17 +637,17 @@ std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save) save << " " << _nbLayers << " " << _thickness << " " << _stretchFactor - << " " << _ignoreFaceIds.size(); - for ( unsigned i = 0; i < _ignoreFaceIds.size(); ++i ) - save << " " << _ignoreFaceIds[i]; + << " " << _ignoreBndShapeIds.size(); + for ( unsigned i = 0; i < _ignoreBndShapeIds.size(); ++i ) + save << " " << _ignoreBndShapeIds[i]; return save; } // -------------------------------------------------------------------------------- std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load) { int nbFaces, faceID; load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces; - while ( _ignoreFaceIds.size() < nbFaces && load >> faceID ) - _ignoreFaceIds.push_back( faceID ); + while ( _ignoreBndShapeIds.size() < nbFaces && load >> faceID ) + _ignoreBndShapeIds.push_back( faceID ); return load; } // -------------------------------------------------------------------------------- bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh* theMesh, @@ -669,38 +727,102 @@ 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); if ( !( ok = ( edges.size() > 0 ))) return dir; // get average dir of edges going fromV - gp_Vec edgeDir; - for ( unsigned i = 0; i < edges.size(); ++i ) - { - edgeDir = getEdgeDir( edges[i], fromV ); - double size2 = edgeDir.SquareMagnitude(); - if ( size2 > numeric_limits::min() ) - edgeDir /= sqrt( size2 ); - else - ok = false; - dir += edgeDir.XYZ(); - } + gp_XYZ edgeDir; + //if ( edges.size() > 1 ) + for ( unsigned i = 0; i < edges.size(); ++i ) + { + edgeDir = getEdgeDir( edges[i], fromV ); + double size2 = edgeDir.SquareModulus(); + if ( size2 > numeric_limits::min() ) + edgeDir /= sqrt( size2 ); + else + ok = false; + dir += edgeDir; + } gp_XYZ fromEdgeDir = getFaceDir( F, edges[0], node, helper, ok ); - if ( edges.size() == 1 || dir.SquareModulus() < 1e-10) + if ( edges.size() == 1 ) dir = fromEdgeDir; + else if ( dir.SquareModulus() < 0.1 ) // ~< 20 degrees + dir = fromEdgeDir + getFaceDir( F, edges[1], node, helper, ok ); else if ( dir * fromEdgeDir < 0 ) dir *= -1; if ( ok ) { //dir /= edges.size(); if ( cosin ) { - double angle = edgeDir.Angle( dir ); + double angle = gp_Vec( edgeDir ).Angle( dir ); *cosin = cos( angle ); } } 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 @@ -712,35 +834,42 @@ 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 ))"< shapesToSmooth; SMESH_MesherHelper helper( *_mesh ); - bool ok; + bool ok = true; for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS ) { @@ -1559,8 +1715,8 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, } case SMDS_TOP_VERTEX: { TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS())); - gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK); - double angle = inFaceDir.Angle( edge._normal ); // [0,PI] + gp_XYZ inFaceDir = getFaceDir( F, V, node, helper, normOK); + double angle = gp_Vec( inFaceDir).Angle( edge._normal ); // [0,PI] edge._cosin = cos( angle ); //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl; break; @@ -1792,14 +1948,14 @@ void _LayerEdge::SetCosin( double cosin ) void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node, vector<_Simplex>& simplices, const set& ingnoreShapes, - const _SolidData* dataToCheckOri) + const _SolidData* dataToCheckOri, + const bool toSort) { - SMESH_MeshEditor editor( _mesh ); SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face); while ( fIt->more() ) { const SMDS_MeshElement* f = fIt->next(); - const TGeomID shapeInd = editor.FindShape( f ); + const TGeomID shapeInd = f->getshapeId(); if ( ingnoreShapes.count( shapeInd )) continue; const int nbNodes = f->NbCornerNodes(); int srcInd = f->GetNodeIndex( node ); @@ -1809,7 +1965,25 @@ void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node, std::swap( nPrev, nNext ); simplices.push_back( _Simplex( nPrev, nNext )); } - simplices.resize( simplices.size() ); + + if ( toSort ) + { + vector<_Simplex> sortedSimplices( simplices.size() ); + sortedSimplices[0] = simplices[0]; + int nbFound = 0; + for ( size_t i = 1; i < simplices.size(); ++i ) + { + for ( size_t j = 1; j < simplices.size(); ++j ) + if ( sortedSimplices[i-1]._nNext == simplices[j]._nPrev ) + { + sortedSimplices[i] = simplices[j]; + nbFound++; + break; + } + } + if ( nbFound == simplices.size() - 1 ) + simplices.swap( sortedSimplices ); + } } //================================================================================ @@ -1900,7 +2074,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) { if ( data._edges[i]->IsOnEdge() ) continue; data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon ); - if ( geomSize > intersecDist ) + if ( geomSize > intersecDist && intersecDist > 0 ) geomSize = intersecDist; } if ( data._stepSize > 0.3 * geomSize ) @@ -1994,7 +2168,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) //================================================================================ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, - int nbSteps, + const int nbSteps, double & distToIntersection) { if ( data._endEdgeToSmooth.empty() ) @@ -2028,20 +2202,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 )) + { + // 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 ")<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 * uMidl < 0. ) + uLast += ( uMidl > 0 ? +2. : -2. ) * M_PI; + const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() ); + + gp_Ax2d axis( center, vec0 ); + gp_Circ2d circ( axis, radius ); + 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 @@ -2235,7 +2675,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++; @@ -2424,7 +2864,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher& searcher, } if ( intFound ) { - if ( dist < segLen*(1.01)) + if ( dist < segLen*(1.01) && dist > -(_len-segLen) ) segmentIntersected = true; if ( distance > dist ) distance = dist, iFace = j; @@ -2627,7 +3067,8 @@ bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface, double lenDelta = 0; if ( _curvature ) { - lenDelta = _curvature->lenDelta( _len ); + //lenDelta = _curvature->lenDelta( _len ); + lenDelta = _curvature->lenDeltaByDist( dist01 ); newPos.ChangeCoord() += _normal * lenDelta; } @@ -2957,6 +3398,15 @@ bool _ViscousBuilder::refine(_SolidData& data) } } + if ( !getMeshDS()->IsEmbeddedMode() ) + // Log node movement + for ( unsigned i = 0; i < data._edges.size(); ++i ) + { + _LayerEdge& edge = *data._edges[i]; + SMESH_TNodeXYZ p ( edge._nodes.back() ); + getMeshDS()->MoveNode( p._node, p.X(), p.Y(), p.Z() ); + } + // TODO: make quadratic prisms and polyhedrons(?) helper.SetElementsOnShape(true); @@ -3036,14 +3486,16 @@ 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 } } } SMESH_MesherHelper helper( *_mesh ); + helper.ToFixNodeParameters( true ); // EDGE's to shrink - map< int, _Shrinker1D > e2shrMap; + map< TGeomID, _Shrinker1D > e2shrMap; // loop on FACES to srink mesh on map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin(); @@ -3052,7 +3504,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); @@ -3065,8 +3516,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(); @@ -3080,12 +3531,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; } @@ -3133,6 +3587,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() ); { @@ -3142,7 +3599,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(); @@ -3161,6 +3620,7 @@ bool _ViscousBuilder::shrink() _Shrinker1D& srinker = e2shrMap[ edgeIndex ]; eShri1D.insert( & srinker ); srinker.AddEdge( edge, helper ); + VISCOUS_3D::ToClearSubWithMain( _mesh->GetSubMesh( edge->_sWOL ), data._solid ); // restore params of nodes on EGDE if the EDGE has been already // srinked while srinking another FACE srinker.RestoreParams(); @@ -3185,8 +3645,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(); @@ -3207,7 +3665,9 @@ 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=*/isConcaveFace); } if ( badNb < oldBadNb ) nbNoImpSteps = 0; @@ -3219,18 +3679,44 @@ 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 - _SrinkShapeListener::ToClearSubMeshWithSolid( sm, data._solid ); + // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh + VISCOUS_3D::ToClearSubWithMain( sm, data._solid ); + + if ( !getMeshDS()->IsEmbeddedMode() ) + // Log node movement + for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) + { + SMESH_TNodeXYZ p ( nodesToSmooth[i]._node ); + getMeshDS()->MoveNode( nodesToSmooth[i]._node, p.X(), p.Y(), p.Z() ); + } - }// loop on FACES to srink mesh on + } // loop on FACES to srink mesh on // Replace source nodes by target nodes in shrinked mesh edges @@ -3416,6 +3902,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() ); @@ -3619,7 +4267,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper ) GeomAdaptor_Curve aCurve(C, f,l); const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l); - int nbExpectNodes = eSubMesh->NbNodes() - e->_nodes.size(); + int nbExpectNodes = eSubMesh->NbNodes(); _initU .reserve( nbExpectNodes ); _normPar.reserve( nbExpectNodes ); _nodes .reserve( nbExpectNodes ); @@ -3727,6 +4375,7 @@ void _Shrinker1D::RestoreParams() } _done = false; } + //================================================================================ /*! * \brief Replace source nodes by target nodes in shrinked mesh edges @@ -3830,6 +4479,8 @@ bool _ViscousBuilder::addBoundaryElements() F = e2f->second.Oriented( TopAbs_FORWARD ); reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED ); if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED ) + reverse = !reverse, F.Reverse(); + if ( helper.IsReversedSubMesh( TopoDS::Face(F) )) reverse = !reverse; } else @@ -3861,12 +4512,28 @@ bool _ViscousBuilder::addBoundaryElements() vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes; vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes; if ( isOnFace ) - for ( unsigned z = 1; z < nn1.size(); ++z ) + for ( size_t z = 1; z < nn1.size(); ++z ) sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] )); else - for ( unsigned z = 1; z < nn1.size(); ++z ) + for ( size_t z = 1; z < nn1.size(); ++z ) sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z])); } + + // Make edges + for ( int isFirst = 0; isFirst < 2; ++isFirst ) + { + _LayerEdge* edge = isFirst ? ledges.front() : ledges.back(); + if ( !edge->_sWOL.IsNull() && edge->_sWOL.ShapeType() == TopAbs_EDGE ) + { + vector< const SMDS_MeshNode*>& nn = edge->_nodes; + if ( nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() ) + continue; + helper.SetSubShape( edge->_sWOL ); + helper.SetElementsOnShape( true ); + for ( size_t z = 1; z < nn.size(); ++z ) + helper.AddEdge( nn[z-1], nn[z] ); + } + } } }