X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=7b5cbba4128d26e3d019b88e266ad560bac33378;hp=6042aa4220d4c28d1f7fc046d9eaaa815c314e03;hb=10191484fe88a27e962b8e4b57e09d390d8705c7;hpb=e74c29b7867ca230102fc580d831dda6b81b3b4c diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 6042aa422..7b5cbba41 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -95,7 +95,7 @@ #include #ifdef _DEBUG_ -//#define __myDEBUG +#define __myDEBUG //#define __NOT_INVALIDATE_BAD_SMOOTH //#define __NODES_AT_POS #endif @@ -114,7 +114,7 @@ namespace VISCOUS_3D enum UIndex { U_TGT = 1, U_SRC, LEN_TGT }; const double theMinSmoothCosin = 0.1; - const double theSmoothThickToElemSizeRatio = 0.3; + const double theSmoothThickToElemSizeRatio = 0.6; const double theMinSmoothTriaAngle = 30; const double theMinSmoothQuadAngle = 45; @@ -427,20 +427,21 @@ namespace VISCOUS_3D enum EFlags { TO_SMOOTH = 0x0000001, MOVED = 0x0000002, // set by _neibors[i]->SetNewLength() - SMOOTHED = 0x0000004, // set by this->Smooth() + SMOOTHED = 0x0000004, // set by _LayerEdge::Smooth() DIFFICULT = 0x0000008, // near concave VERTEX ON_CONCAVE_FACE = 0x0000010, BLOCKED = 0x0000020, // not to inflate any more INTERSECTED = 0x0000040, // close intersection with a face found NORMAL_UPDATED = 0x0000080, - MARKED = 0x0000100, // local usage - MULTI_NORMAL = 0x0000200, // a normal is invisible by some of surrounding faces - NEAR_BOUNDARY = 0x0000400, // is near FACE boundary forcing smooth - SMOOTHED_C1 = 0x0000800, // is on _eosC1 - DISTORTED = 0x0001000, // was bad before smoothing - RISKY_SWOL = 0x0002000, // SWOL is parallel to a source FACE - SHRUNK = 0x0004000, // target node reached a tgt position while shrink() - UNUSED_FLAG = 0x0100000 // to add use flags after + UPD_NORMAL_CONV = 0x0000100, // to update normal on boundary of concave FACE + MARKED = 0x0000200, // local usage + MULTI_NORMAL = 0x0000400, // a normal is invisible by some of surrounding faces + NEAR_BOUNDARY = 0x0000800, // is near FACE boundary forcing smooth + SMOOTHED_C1 = 0x0001000, // is on _eosC1 + DISTORTED = 0x0002000, // was bad before smoothing + RISKY_SWOL = 0x0004000, // SWOL is parallel to a source FACE + SHRUNK = 0x0008000, // target node reached a tgt position while shrink() + UNUSED_FLAG = 0x0100000 // to add user flags after }; bool Is ( int flag ) const { return _flags & flag; } void Set ( int flag ) { _flags |= flag; } @@ -497,11 +498,14 @@ namespace VISCOUS_3D const gp_XYZ& PrevPos() const { return _pos[ _pos.size() - 2 ]; } gp_XYZ PrevCheckPos( _EdgesOnShape* eos=0 ) const; gp_Ax1 LastSegment(double& segLen, _EdgesOnShape& eos) const; - gp_XY LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const; + gp_XY LastUV( const TopoDS_Face& F, _EdgesOnShape& eos, int which=-1 ) const; bool IsOnEdge() const { return _2neibors; } + bool IsOnFace() const { return ( _nodes[0]->GetPosition()->GetDim() == 2 ); } + int BaseShapeDim() const { return _nodes[0]->GetPosition()->GetDim(); } gp_XYZ Copy( _LayerEdge& other, _EdgesOnShape& eos, SMESH_MesherHelper& helper ); void SetCosin( double cosin ); void SetNormal( const gp_XYZ& n ) { _normal = n; } + void SetMaxLen( double l ) { _maxLen = l; } int NbSteps() const { return _pos.size() - 1; } // nb inlation steps bool IsNeiborOnEdge( const _LayerEdge* edge ) const; void SetSmooLen( double len ) { // set _len at which smoothing is needed @@ -664,6 +668,8 @@ namespace VISCOUS_3D _SolidData* _data; // parent SOLID + _LayerEdge* operator[](size_t i) const { return (_LayerEdge*) _edges[i]; } + size_t size() const { return _edges.size(); } TopAbs_ShapeEnum ShapeType() const { return _shape.IsNull() ? TopAbs_SHAPE : _shape.ShapeType(); } TopAbs_ShapeEnum SWOLType() const @@ -678,7 +684,7 @@ namespace VISCOUS_3D //-------------------------------------------------------------------------------- /*! - * \brief Convex FACE whose radius of curvature is less than the thickness of + * \brief Convex FACE whose radius of curvature is less than the thickness of * layers. It is used to detect distortion of prisms based on a convex * FACE and to update normals to enable further increasing the thickness */ @@ -692,7 +698,14 @@ namespace VISCOUS_3D // map a sub-shape to _SolidData::_edgesOnShape map< TGeomID, _EdgesOnShape* > _subIdToEOS; + bool _isTooCurved; bool _normalsFixed; + bool _normalsFixedOnBorders; // used in putOnOffsetSurface() + + double GetMaxCurvature( _SolidData& data, + _EdgesOnShape& eof, + BRepLProp_SLProps& surfProp, + SMESH_MesherHelper& helper); bool GetCenterOfCurvature( _LayerEdge* ledge, BRepLProp_SLProps& surfProp, @@ -757,8 +770,6 @@ namespace VISCOUS_3D int _nbShapesToSmooth; - //map< TGeomID,Handle(Geom_Curve)> _edge2curve; - vector< _CollisionEdges > _collisionEdges; set< TGeomID > _concaveFaces; @@ -836,6 +847,8 @@ namespace VISCOUS_3D void Append( const gp_Pnt& center, _LayerEdge* ledge ) { + if ( ledge->Is( _LayerEdge::MULTI_NORMAL )) + return; if ( _curvaCenters.size() > 0 ) _segLength2.push_back( center.SquareDistance( _curvaCenters.back() )); _curvaCenters.push_back( center ); @@ -870,6 +883,7 @@ namespace VISCOUS_3D const gp_XY& uvToFix, const double refSign ); }; + struct PyDump; //-------------------------------------------------------------------------------- /*! * \brief Builder of viscous layers @@ -942,12 +956,15 @@ namespace VISCOUS_3D void makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper& ); void putOnOffsetSurface( _EdgesOnShape& eos, int infStep, vector< _EdgesOnShape* >& eosC1, - int smooStep=0, bool moveAll=false ); + int smooStep=0, int moveAll=false ); void findCollisionEdges( _SolidData& data, SMESH_MesherHelper& helper ); + void findEdgesToUpdateNormalNearConvexFace( _ConvexFace & convFace, + _SolidData& data, + SMESH_MesherHelper& helper ); void limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& helper ); void limitMaxLenByCurvature( _LayerEdge* e1, _LayerEdge* e2, _EdgesOnShape& eos1, _EdgesOnShape& eos2, - SMESH_MesherHelper& helper ); + const bool isSmoothable ); bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper, int stepNb, double stepSize ); bool updateNormalsOfConvexFaces( _SolidData& data, SMESH_MesherHelper& helper, @@ -987,6 +1004,7 @@ namespace VISCOUS_3D TopTools_MapOfShape _shrinkedFaces; int _tmpFaceID; + PyDump* _pyDump; }; //-------------------------------------------------------------------------------- /*! @@ -1034,6 +1052,7 @@ namespace VISCOUS_3D size_t _iSeg[2]; // index of segment where extreme tgt node is projected _EdgesOnShape& _eos; double _curveLen; // length of the EDGE + std::pair _eToSmooth[2]; // indices of _LayerEdge's in _eos static Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge& E, _EdgesOnShape& eos, @@ -1047,36 +1066,31 @@ namespace VISCOUS_3D bool Perform(_SolidData& data, Handle(ShapeAnalysis_Surface)& surface, const TopoDS_Face& F, - SMESH_MesherHelper& helper ) - { - if ( _leParams.empty() || ( !isAnalytic() && _offPoints.empty() )) - prepare( data ); + SMESH_MesherHelper& helper ); - if ( isAnalytic() ) - return smoothAnalyticEdge( data, surface, F, helper ); - else - return smoothComplexEdge ( data, surface, F, helper ); - } void prepare(_SolidData& data ); + void findEdgesToSmooth(); + + bool isToSmooth( int iE ); + bool smoothAnalyticEdge( _SolidData& data, Handle(ShapeAnalysis_Surface)& surface, const TopoDS_Face& F, SMESH_MesherHelper& helper); - bool smoothComplexEdge( _SolidData& data, Handle(ShapeAnalysis_Surface)& surface, const TopoDS_Face& F, SMESH_MesherHelper& helper); - gp_XYZ getNormalNormal( const gp_XYZ & normal, const gp_XYZ& edgeDir); - _LayerEdge* getLEdgeOnV( bool is2nd ) { return _eos._edges[ is2nd ? _eos._edges.size()-1 : 0 ]->_2neibors->_edges[ is2nd ]; } bool isAnalytic() const { return !_anaCurve.IsNull(); } + + void offPointsToPython() const; // debug }; //-------------------------------------------------------------------------------- /*! @@ -1198,8 +1212,8 @@ namespace VISCOUS_3D //================================================================================ // StdMeshers_ViscousLayers hypothesis // -StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen) - :SMESH_Hypothesis(hypId, studyId, gen), +StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, SMESH_Gen* gen) + :SMESH_Hypothesis(hypId, gen), _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1), _method( SURF_OFFSET_SMOOTH ) { @@ -1596,8 +1610,8 @@ namespace VISCOUS_3D // look for two neighbor not in-FACE nodes of face for ( int i = 0; i < 2; ++i ) { - if ( nNext[i]->GetPosition()->GetDim() != 2 && - nNext[i]->GetID() < nodeOnEdge->GetID() ) + if (( nNext[i]->GetPosition()->GetDim() != 2 ) && + ( nodeOnEdge->GetPosition()->GetDim() == 0 || nNext[i]->GetID() < nodeOnEdge->GetID() )) { // look for an in-FACE node for ( int iN = 0; iN < nbN; ++iN ) @@ -1663,18 +1677,19 @@ namespace VISCOUS_3D // HOWTO use: run python commands written in a console to see // construction steps of viscous layers #ifdef __myDEBUG - ofstream* py; - int theNbPyFunc; - struct PyDump { + ostream* py; + int theNbPyFunc; + struct PyDump + { PyDump(SMESH_Mesh& m) { int tag = 3 + m.GetId(); const char* fname = "/tmp/viscous.py"; cout << "execfile('"< ostream & operator<<( const T &anything ) { return *this ; } + }; + void Pause() { py = &_mystream; } + void Resume() { py = _pyStream; } + MyStream _mystream; + ostream* _pyStream; }; #define dumpFunction(f) { _dumpFunction(f, __LINE__);} #define dumpMove(n) { _dumpMove(n, __LINE__);} @@ -1710,7 +1733,7 @@ namespace VISCOUS_3D #else - struct PyDump { PyDump(SMESH_Mesh&) {} void Finish() {} }; + struct PyDump { PyDump(SMESH_Mesh&) {} void Finish() {} void Pause() {} void Resume() {} }; #define dumpFunction(f) f #define dumpMove(n) #define dumpMoveComm(n,txt) @@ -1852,6 +1875,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh, return SMESH_ComputeErrorPtr(); // everything already computed PyDump debugDump( theMesh ); + _pyDump = &debugDump; // TODO: ignore already computed SOLIDs if ( !findSolidsWithLayers()) @@ -1865,15 +1889,19 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh, size_t iSD = 0; for ( iSD = 0; iSD < _sdVec.size(); ++iSD ) // find next SOLID to compute if ( _sdVec[iSD]._before.IsEmpty() && + !_sdVec[iSD]._solid.IsNull() && _sdVec[iSD]._n2eMap.empty() ) break; if ( ! makeLayer(_sdVec[iSD]) ) // create _LayerEdge's return _error; - if ( _sdVec[iSD]._n2eMap.size() == 0 ) + if ( _sdVec[iSD]._n2eMap.size() == 0 ) // no layers in a SOLID + { + _sdVec[iSD]._solid.Nullify(); continue; - + } + if ( ! inflate(_sdVec[iSD]) ) // increase length of _LayerEdge's return _error; @@ -2317,7 +2345,7 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith) for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() ) { TGeomID vID = getMeshDS()->ShapeToIndex( vIt.Value() ); - bool noShrinkV = false; + bool noShrinkV = false, noShrinkIfAdjMeshed = false; if ( iSolid < _sdVec.size() ) { @@ -2327,7 +2355,8 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith) i2S = _sdVec[i ]._shrinkShape2Shape.find( vID ); i2SAdj = _sdVec[iSolid]._shrinkShape2Shape.find( vID ); if ( i2SAdj == _sdVec[iSolid]._shrinkShape2Shape.end() ) - noShrinkV = ( i2S->second.ShapeType() == TopAbs_EDGE || isStructured ); + noShrinkV = (( isStructured ) || + ( noShrinkIfAdjMeshed = i2S->second.ShapeType() == TopAbs_EDGE )); else noShrinkV = ( ! i2S->second.IsSame( i2SAdj->second )); } @@ -2339,8 +2368,31 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith) else { // the adjacent SOLID has NO layers at all - noShrinkV = ( isStructured || - _sdVec[i]._shrinkShape2Shape[ vID ].ShapeType() == TopAbs_EDGE ); + if ( isStructured ) + { + noShrinkV = true; + } + else + { + noShrinkV = noShrinkIfAdjMeshed = + ( _sdVec[i]._shrinkShape2Shape[ vID ].ShapeType() == TopAbs_EDGE ); + } + } + + if ( noShrinkV && noShrinkIfAdjMeshed ) + { + // noShrinkV if FACEs in the adjacent SOLID are meshed + PShapeIteratorPtr fIt = helper.GetAncestors( _sdVec[i]._shrinkShape2Shape[ vID ], + *_mesh, TopAbs_FACE, &solid ); + while ( fIt->more() ) + { + const TopoDS_Shape* f = fIt->next(); + if ( !f->IsSame( fWOL )) + { + noShrinkV = ! _mesh->GetSubMesh( *f )->IsEmpty(); + break; + } + } } if ( noShrinkV ) _sdVec[i]._noShrinkShapes.insert( vID ); @@ -2773,8 +2825,6 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) { SMESH_MesherHelper helper( *_mesh ); - const int nbTestPnt = 5; // on a FACE sub-shape - BRepLProp_SLProps surfProp( 2, 1e-6 ); data._convexFaces.clear(); @@ -2786,58 +2836,27 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) continue; TopoDS_Face F = TopoDS::Face( eof._shape ); - SMESH_subMesh * sm = eof._subMesh; const TGeomID faceID = eof._shapeID; BRepAdaptor_Surface surface( F, false ); surfProp.SetSurface( surface ); - bool isTooCurved = false; - _ConvexFace cnvFace; - const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. ); - SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true); - while ( smIt->more() ) - { - sm = smIt->next(); - const TGeomID subID = sm->GetId(); - // find _LayerEdge's of a sub-shape - _EdgesOnShape* eos; - if (( eos = data.GetShapeEdges( subID ))) - cnvFace._subIdToEOS.insert( make_pair( subID, eos )); - else - continue; - // check concavity and curvature and limit data._stepSize - const double minCurvature = - 1. / ( eos->_hyp.GetTotalThickness() * ( 1 + theThickToIntersection )); - size_t iStep = Max( 1, eos->_edges.size() / nbTestPnt ); - for ( size_t i = 0; i < eos->_edges.size(); i += iStep ) - { - gp_XY uv = helper.GetNodeUV( F, eos->_edges[ i ]->_nodes[0] ); - surfProp.SetParameters( uv.X(), uv.Y() ); - if ( !surfProp.IsCurvatureDefined() ) - continue; - if ( surfProp.MaxCurvature() * oriFactor > minCurvature ) - { - limitStepSize( data, 0.9 / surfProp.MaxCurvature() * oriFactor ); - isTooCurved = true; - } - if ( surfProp.MinCurvature() * oriFactor > minCurvature ) - { - limitStepSize( data, 0.9 / surfProp.MinCurvature() * oriFactor ); - isTooCurved = true; - } - } - } // loop on sub-shapes of the FACE + cnvFace._face = F; + cnvFace._normalsFixed = false; + cnvFace._isTooCurved = false; - if ( !isTooCurved ) continue; + double maxCurvature = cnvFace.GetMaxCurvature( data, eof, surfProp, helper ); + if ( maxCurvature > 0 ) + { + limitStepSize( data, 0.9 / maxCurvature ); + findEdgesToUpdateNormalNearConvexFace( cnvFace, data, helper ); + } + if ( !cnvFace._isTooCurved ) continue; _ConvexFace & convFace = data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second; - convFace._face = F; - convFace._normalsFixed = false; - // skip a closed surface (data._convexFaces is useful anyway) bool isClosedF = false; helper.SetSubShape( F ); @@ -2850,6 +2869,7 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) if ( isClosedF ) { // limit _LayerEdge::_maxLen on the FACE + const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. ); const double minCurvature = 1. / ( eof._hyp.GetTotalThickness() * ( 1 + theThickToIntersection )); map< TGeomID, _EdgesOnShape* >::iterator id2eos = cnvFace._subIdToEOS.find( faceID ); @@ -2861,14 +2881,13 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) _LayerEdge* ledge = eos._edges[ i ]; gp_XY uv = helper.GetNodeUV( F, ledge->_nodes[0] ); surfProp.SetParameters( uv.X(), uv.Y() ); - if ( !surfProp.IsCurvatureDefined() ) - continue; - - if ( surfProp.MaxCurvature() * oriFactor > minCurvature ) - ledge->_maxLen = Min( ledge->_maxLen, 1. / surfProp.MaxCurvature() * oriFactor ); - - if ( surfProp.MinCurvature() * oriFactor > minCurvature ) - ledge->_maxLen = Min( ledge->_maxLen, 1. / surfProp.MinCurvature() * oriFactor ); + if ( surfProp.IsCurvatureDefined() ) + { + double curvature = Max( surfProp.MaxCurvature() * oriFactor, + surfProp.MinCurvature() * oriFactor ); + if ( curvature > minCurvature ) + ledge->SetMaxLen( Min( ledge->_maxLen, 1. / curvature )); + } } } continue; @@ -2888,7 +2907,12 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) for ( size_t j = 0; j < ledge->_simplices.size(); ++j ) if ( ledge->_simplices[j]._nNext->GetPosition()->GetDim() < 2 ) { - convFace._simplexTestEdges.push_back( ledge ); + // do not select _LayerEdge's neighboring sharp EDGEs + bool sharpNbr = false; + for ( size_t iN = 0; iN < ledge->_neibors.size() && !sharpNbr; ++iN ) + sharpNbr = ( ledge->_neibors[iN]->_cosin > theMinSmoothCosin ); + if ( !sharpNbr ) + convFace._simplexTestEdges.push_back( ledge ); break; } } @@ -2938,22 +2962,11 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) // define allowed thickness computeGeomSize( data ); // compute data._geomSize and _LayerEdge::_maxLen - data._maxThickness = 0; - data._minThickness = 1e100; - list< const StdMeshers_ViscousLayers* >::iterator hyp = data._hyps.begin(); - for ( ; hyp != data._hyps.end(); ++hyp ) - { - data._maxThickness = Max( data._maxThickness, (*hyp)->GetTotalThickness() ); - data._minThickness = Min( data._minThickness, (*hyp)->GetTotalThickness() ); - } - //const double tgtThick = /*Min( 0.5 * data._geomSize, */data._maxThickness; // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's // boundary inclined to the shape at a sharp angle - //list< TGeomID > shapesToSmooth; TopTools_MapOfShape edgesOfSmooFaces; - SMESH_MesherHelper helper( *_mesh ); bool ok = true; @@ -2968,32 +2981,34 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) continue; double tgtThick = eos._hyp.GetTotalThickness(); - TopExp_Explorer eExp( edgesByGeom[iS]._shape, TopAbs_EDGE ); - for ( ; eExp.More() && !eos._toSmooth; eExp.Next() ) + SMESH_subMeshIteratorPtr subIt = eos._subMesh->getDependsOnIterator(/*includeSelf=*/false ); + while ( subIt->more() && !eos._toSmooth ) { - TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() ); - vector<_LayerEdge*>& eE = edgesByGeom[ iE ]._edges; - if ( eE.empty() ) continue; + TGeomID iSub = subIt->next()->GetId(); + const vector<_LayerEdge*>& eSub = edgesByGeom[ iSub ]._edges; + if ( eSub.empty() ) continue; double faceSize; - for ( size_t i = 0; i < eE.size() && !eos._toSmooth; ++i ) - if ( eE[i]->_cosin > theMinSmoothCosin ) + for ( size_t i = 0; i < eSub.size() && !eos._toSmooth; ++i ) + if ( eSub[i]->_cosin > theMinSmoothCosin ) { - SMDS_ElemIteratorPtr fIt = eE[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face); + SMDS_ElemIteratorPtr fIt = eSub[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face); while ( fIt->more() && !eos._toSmooth ) { const SMDS_MeshElement* face = fIt->next(); if ( face->getshapeId() == eos._shapeID && - getDistFromEdge( face, eE[i]->_nodes[0], faceSize )) + getDistFromEdge( face, eSub[i]->_nodes[0], faceSize )) { - eos._toSmooth = needSmoothing( eE[i]->_cosin, tgtThick, faceSize ); + eos._toSmooth = needSmoothing( eSub[i]->_cosin, + tgtThick * eSub[i]->_lenFactor, + faceSize); } } } } if ( eos._toSmooth ) { - for ( eExp.ReInit(); eExp.More(); eExp.Next() ) + for ( TopExp_Explorer eExp( edgesByGeom[iS]._shape, TopAbs_EDGE ); eExp.More(); eExp.Next() ) edgesOfSmooFaces.Add( eExp.Current() ); data.PrepareEdgesToSmoothOnFace( &edgesByGeom[iS], /*substituteSrcNodes=*/false ); @@ -3037,16 +3052,16 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) if ( endSeg->getshapeId() == (int) iS ) { double segLen = - SMESH_TNodeXYZ( endSeg->GetNode(0) ).Distance( endSeg->GetNode(1 )); - eos._toSmooth = needSmoothing( cosinAbs, tgtThick, segLen ); + SMESH_TNodeXYZ( endSeg->GetNode( 0 )).Distance( endSeg->GetNode( 1 )); + eos._toSmooth = needSmoothing( cosinAbs, tgtThick * eV[0]->_lenFactor, segLen ); } } if ( eos._toSmooth ) { eos._edgeSmoother = new _Smoother1D( curve, eos ); - for ( size_t i = 0; i < eos._edges.size(); ++i ) - eos._edges[i]->Set( _LayerEdge::TO_SMOOTH ); + // for ( size_t i = 0; i < eos._edges.size(); ++i ) + // eos._edges[i]->Set( _LayerEdge::TO_SMOOTH ); } } } @@ -3062,7 +3077,8 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) if ( !eos._hyp.ToSmooth() ) for ( size_t i = 0; i < eos._edges.size(); ++i ) - eos._edges[i]->SetCosin( 0 ); + //eos._edges[i]->SetCosin( 0 ); // keep _cosin to use in limitMaxLenByCurvature() + eos._edges[i]->_lenFactor = 1; } @@ -3126,19 +3142,21 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) { _EdgesOnShape* eof = data.GetShapeEdges( *face ); if ( !eof ) continue; // other solid - if ( !eos.HasC1( eoe )) - { - eos._eosC1.push_back( eoe ); - eoe->_toSmooth = false; - data.PrepareEdgesToSmoothOnFace( eoe, /*substituteSrcNodes=*/false ); - } - if ( eos._shapeID != eof->_shapeID && !eos.HasC1( eof )) + if ( eos._shapeID == eof->_shapeID ) continue; + if ( !eos.HasC1( eof )) { + // check the FACEs eos._eosC1.push_back( eof ); eof->_toSmooth = false; data.PrepareEdgesToSmoothOnFace( eof, /*substituteSrcNodes=*/false ); smQueue.push_back( eof->_subMesh ); } + if ( !eos.HasC1( eoe )) + { + eos._eosC1.push_back( eoe ); + eoe->_toSmooth = false; + data.PrepareEdgesToSmoothOnFace( eoe, /*substituteSrcNodes=*/false ); + } } } } @@ -3294,6 +3312,7 @@ void _ViscousBuilder::setShapeData( _EdgesOnShape& eos, if ( eos.ShapeType() == TopAbs_FACE ) // get normals to elements on a FACE { SMESHDS_SubMesh* smDS = sm->GetSubMeshDS(); + if ( !smDS ) return; eos._faceNormals.resize( smDS->NbElements() ); SMDS_ElemIteratorPtr eIt = smDS->GetElements(); @@ -3425,11 +3444,12 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, } // find _normal + bool fromVonF = false; if ( useGeometry ) { - bool fromVonF = ( eos.ShapeType() == TopAbs_VERTEX && - eos.SWOLType() == TopAbs_FACE && - totalNbFaces > 1 ); + fromVonF = ( eos.ShapeType() == TopAbs_VERTEX && + eos.SWOLType() == TopAbs_FACE && + totalNbFaces > 1 ); if ( onShrinkShape && !fromVonF ) // one of faces the node is on has no layers { @@ -3531,14 +3551,19 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, break; } case TopAbs_VERTEX: { - //if ( eos.SWOLType() != TopAbs_FACE ) // else _cosin is set by getFaceDir() + if ( fromVonF ) + { + getFaceDir( TopoDS::Face( eos._sWOL ), TopoDS::Vertex( eos._shape ), + node, helper, normOK, &edge._cosin ); + } + else if ( eos.SWOLType() != TopAbs_FACE ) // else _cosin is set by getFaceDir() { TopoDS_Vertex V = TopoDS::Vertex( eos._shape ); gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK ); double angle = inFaceDir.Angle( edge._normal ); // [0,PI] edge._cosin = Cos( angle ); if ( totalNbFaces > 2 || helper.IsSeamShape( node->getshapeId() )) - for ( int iF = totalNbFaces-2; iF >=0; --iF ) + for ( int iF = 1; iF < totalNbFaces; ++iF ) { F = face2Norm[ iF ].first; inFaceDir = getFaceDir( F, V, node, helper, normOK=true ); @@ -3568,7 +3593,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, getMeshDS()->RemoveFreeNode( edge._nodes.back(), 0, /*fromGroups=*/false ); edge._nodes.resize( 1 ); edge._normal.SetCoord( 0,0,0 ); - edge._maxLen = 0; + edge.SetMaxLen( 0 ); } // Set the rest data @@ -3705,7 +3730,7 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node, isOK = false; return p.XYZ(); } - Quantity_Parameter U,V; + Standard_Real U,V; projector.LowerDistanceParameters(U,V); uv.SetCoord( U,V ); } @@ -4116,6 +4141,8 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, { if ( eos.ShapeType() != TopAbs_EDGE ) return; + if ( _curvature && Is( SMOOTHED_C1 )) + return; gp_XYZ pos = SMESH_TNodeXYZ( _nodes[0] ); gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 ); @@ -4369,7 +4396,7 @@ void _ViscousBuilder::computeGeomSize( _SolidData& data ) _EdgesOnShape& eos = data._edgesOnShape[ iS ]; if ( eos._edges.empty() ) continue; - // get neighbor faces intersection with which should not be considered since + // get neighbor faces, intersection with which should not be considered since // collisions are avoided by means of smoothing set< TGeomID > neighborFaces; if ( eos._hyp.ToSmooth() ) @@ -4389,16 +4416,88 @@ void _ViscousBuilder::computeGeomSize( _SolidData& data ) for ( size_t i = 0; i < eos._edges.size(); ++i ) { if ( eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; - eos._edges[i]->_maxLen = thinkness; + eos._edges[i]->SetMaxLen( thinkness ); eos._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon, eos, &face ); if ( intersecDist > 0 && face ) { data._geomSize = Min( data._geomSize, intersecDist ); if ( !neighborFaces.count( face->getshapeId() )) - eos._edges[i]->_maxLen = Min( thinkness, intersecDist / ( face->GetID() < 0 ? 3. : 2. )); + eos[i]->SetMaxLen( Min( thinkness, intersecDist / ( face->GetID() < 0 ? 3. : 2. ))); } } } + + data._maxThickness = 0; + data._minThickness = 1e100; + list< const StdMeshers_ViscousLayers* >::iterator hyp = data._hyps.begin(); + for ( ; hyp != data._hyps.end(); ++hyp ) + { + data._maxThickness = Max( data._maxThickness, (*hyp)->GetTotalThickness() ); + data._minThickness = Min( data._minThickness, (*hyp)->GetTotalThickness() ); + } + + // Limit inflation step size by geometry size found by intersecting + // normals of _LayerEdge's with mesh faces + if ( data._stepSize > 0.3 * data._geomSize ) + limitStepSize( data, 0.3 * data._geomSize ); + + if ( data._stepSize > data._minThickness ) + limitStepSize( data, data._minThickness ); + + + // ------------------------------------------------------------------------- + // Detect _LayerEdge which can't intersect with opposite or neighbor layer, + // so no need in detecting intersection at each inflation step + // ------------------------------------------------------------------------- + + int nbSteps = data._maxThickness / data._stepSize; + if ( nbSteps < 3 || nbSteps * data._n2eMap.size() < 100000 ) + return; + + vector< const SMDS_MeshElement* > closeFaces; + int nbDetected = 0; + + for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS ) + { + _EdgesOnShape& eos = data._edgesOnShape[ iS ]; + if ( eos._edges.empty() || eos.ShapeType() != TopAbs_FACE ) + continue; + + for ( size_t i = 0; i < eos.size(); ++i ) + { + SMESH_NodeXYZ p( eos[i]->_nodes[0] ); + double radius = data._maxThickness + 2 * eos[i]->_maxLen; + closeFaces.clear(); + searcher->GetElementsInSphere( p, radius, SMDSAbs_Face, closeFaces ); + + bool toIgnore = true; + for ( size_t iF = 0; iF < closeFaces.size() && toIgnore; ++iF ) + if ( !( toIgnore = ( closeFaces[ iF ]->getshapeId() == eos._shapeID || + data._ignoreFaceIds.count( closeFaces[ iF ]->getshapeId() )))) + { + // check if a _LayerEdge will inflate in a direction opposite to a direction + // toward a close face + bool allBehind = true; + for ( int iN = 0; iN < closeFaces[ iF ]->NbCornerNodes() && allBehind; ++iN ) + { + SMESH_NodeXYZ pi( closeFaces[ iF ]->GetNode( iN )); + allBehind = (( pi - p ) * eos[i]->_normal < 0.1 * data._stepSize ); + } + toIgnore = allBehind; + } + + + if ( toIgnore ) // no need to detect intersection + { + eos[i]->Set( _LayerEdge::INTERSECTED ); + ++nbDetected; + } + } + } + + debugMsg( "Nb LE to intersect " << data._n2eMap.size()-nbDetected << ", ignore " << nbDetected ); + + return; } //================================================================================ @@ -4411,24 +4510,20 @@ bool _ViscousBuilder::inflate(_SolidData& data) { SMESH_MesherHelper helper( *_mesh ); - // Limit inflation step size by geometry size found by itersecting - // normals of _LayerEdge's with mesh faces - if ( data._stepSize > 0.3 * data._geomSize ) - limitStepSize( data, 0.3 * data._geomSize ); - const double tgtThick = data._maxThickness; - if ( data._stepSize > data._minThickness ) - limitStepSize( data, data._minThickness ); if ( data._stepSize < 1. ) data._epsilon = data._stepSize * 1e-7; debugMsg( "-- geomSize = " << data._geomSize << ", stepSize = " << data._stepSize ); + _pyDump->Pause(); findCollisionEdges( data, helper ); limitMaxLenByCurvature( data, helper ); + _pyDump->Resume(); + // limit length of _LayerEdge's around MULTI_NORMAL _LayerEdge's for ( size_t i = 0; i < data._edgesOnShape.size(); ++i ) if ( data._edgesOnShape[i].ShapeType() == TopAbs_VERTEX && @@ -4533,6 +4628,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) break; } #endif + // new step size limitStepSize( data, 0.25 * distToIntersection ); if ( data._stepSizeNodes[0] ) @@ -4854,6 +4950,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, _LayerEdge* edge = eos._edges[i]; if ( edge->_nodes.size() < 2 ) continue; SMESH_TNodeXYZ tgtXYZ = edge->_nodes.back(); + //SMESH_TNodeXYZ prevXYZ = edge->_nodes[0]; gp_XYZ prevXYZ = edge->PrevCheckPos( &eos ); //const gp_XYZ& prevXYZ = edge->PrevPos(); for ( size_t j = 0; j < edge->_simplices.size(); ++j ) @@ -4893,6 +4990,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, const SMDS_MeshElement* intFace = 0; const SMDS_MeshElement* closestFace = 0; _LayerEdge* le = 0; + bool is1stBlocked = true; // dbg for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS ) { _EdgesOnShape& eos = data._edgesOnShape[ iS ]; @@ -4955,7 +5053,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, // ignore intersection of a _LayerEdge based on a _ConvexFace with a face // lying on this _ConvexFace if ( _ConvexFace* convFace = data.GetConvexFace( intFace->getshapeId() )) - if ( convFace->_subIdToEOS.count ( eos._shapeID )) + if ( convFace->_isTooCurved && convFace->_subIdToEOS.count ( eos._shapeID )) continue; // ignore intersection of a _LayerEdge based on a FACE with an element on this FACE @@ -4964,18 +5062,19 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, continue; // ignore intersection with intFace of an adjacent FACE - if ( dist > 0 ) + if ( dist > 0.1 * eos._edges[i]->_len ) { bool toIgnore = false; - if ( eos._edges[i]->Is( _LayerEdge::TO_SMOOTH )) + if ( eos._toSmooth ) { const TopoDS_Shape& S = getMeshDS()->IndexToShape( intFace->getshapeId() ); if ( !S.IsNull() && S.ShapeType() == TopAbs_FACE ) { - TopExp_Explorer edge( eos._shape, TopAbs_EDGE ); - for ( ; !toIgnore && edge.More(); edge.Next() ) - // is adjacent - has a common EDGE - toIgnore = ( helper.IsSubShape( edge.Current(), S )); + TopExp_Explorer sub( eos._shape, + eos.ShapeType() == TopAbs_FACE ? TopAbs_EDGE : TopAbs_VERTEX ); + for ( ; !toIgnore && sub.More(); sub.Next() ) + // is adjacent - has a common EDGE or VERTEX + toIgnore = ( helper.IsSubShape( sub.Current(), S )); if ( toIgnore ) // check angle between normals { @@ -5002,6 +5101,9 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, if ( toBlockInfaltion && dist < ( eos._edges[i]->_len * theThickToIntersection )) { + if ( is1stBlocked ) { is1stBlocked = false; // debug + dumpFunction(SMESH_Comment("blockIntersected") <Set( _LayerEdge::INTERSECTED ); // not to intersect eos._edges[i]->Block( data ); // not to inflate @@ -5035,11 +5137,14 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, } // loop on eos._edges } // loop on data._edgesOnShape + if ( !is1stBlocked ) + dumpFunctionEnd(); + if ( closestFace && le ) { #ifdef __myDEBUG SMDS_MeshElement::iterator nIt = closestFace->begin_nodes(); - cout << "Shortest distance: _LayerEdge nodes: tgt " << le->_nodes.back()->GetID() + cout << "#Shortest distance: _LayerEdge nodes: tgt " << le->_nodes.back()->GetID() << " src " << le->_nodes[0]->GetID()<< ", intersection with face (" << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID() << ") distance = " << distToIntersection<< endl; @@ -5233,7 +5338,8 @@ void _ViscousBuilder::makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper& try { - BRepOffsetAPI_MakeOffsetShape offsetMaker( eos._shape, -offset, Precision::Confusion() ); + BRepOffsetAPI_MakeOffsetShape offsetMaker; + offsetMaker.PerformByJoin( eos._shape, -offset, Precision::Confusion() ); if ( !offsetMaker.IsDone() ) return; TopExp_Explorer fExp( offsetMaker.Shape(), TopAbs_FACE ); @@ -5260,7 +5366,7 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, int infStep, vector< _EdgesOnShape* >& eosC1, int smooStep, - bool moveAll ) + int moveAll ) { _EdgesOnShape * eof = & eos; if ( eos.ShapeType() != TopAbs_FACE ) // eos is a boundary of C1 FACE, look for the FACE eos @@ -5291,8 +5397,13 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, edge->Unset( _LayerEdge::MARKED ); if ( edge->Is( _LayerEdge::BLOCKED ) || !edge->_curvature ) continue; - if ( !moveAll && !edge->Is( _LayerEdge::MOVED )) + if ( moveAll == _LayerEdge::UPD_NORMAL_CONV ) + { + if ( !edge->Is( _LayerEdge::UPD_NORMAL_CONV )) continue; + } + else if ( !moveAll && !edge->Is( _LayerEdge::MOVED )) + continue; int nbBlockedAround = 0; for ( size_t iN = 0; iN < edge->_neibors.size(); ++iN ) @@ -5302,7 +5413,7 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, gp_Pnt tgtP = SMESH_TNodeXYZ( edge->_nodes.back() ); gp_Pnt2d uv = eof->_offsetSurf->NextValueOfUV( edge->_curvature->_uv, tgtP, preci ); - if ( eof->_offsetSurf->Gap() > edge->_len ) continue; // NextValueOfUV() bug + if ( eof->_offsetSurf->Gap() > edge->_len ) continue; // NextValueOfUV() bug edge->_curvature->_uv = uv; if ( eof->_offsetSurf->Gap() < 10 * preci ) continue; // same pos @@ -5321,9 +5432,15 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, edge->_pos.back() = newP; edge->Set( _LayerEdge::MARKED ); + if ( moveAll == _LayerEdge::UPD_NORMAL_CONV ) + { + edge->_normal = ( newP - prevP ).Normalized(); + } } } + + #ifdef _DEBUG_ // dumpMove() for debug size_t i = 0; @@ -5332,7 +5449,7 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, break; if ( i < eos._edges.size() ) { - dumpFunction(SMESH_Comment("putOnOffsetSurface_F") << eos._shapeID + dumpFunction(SMESH_Comment("putOnOffsetSurface_S") << eos._shapeID << "_InfStep" << infStep << "_" << smooStep ); for ( ; i < eos._edges.size(); ++i ) { @@ -5342,6 +5459,26 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, dumpFunctionEnd(); } #endif + + _ConvexFace* cnvFace; + if ( moveAll != _LayerEdge::UPD_NORMAL_CONV && + eos.ShapeType() == TopAbs_FACE && + (cnvFace = eos.GetData().GetConvexFace( eos._shapeID )) && + !cnvFace->_normalsFixedOnBorders ) + { + // put on the surface nodes built on FACE boundaries + SMESH_subMeshIteratorPtr smIt = eos._subMesh->getDependsOnIterator(/*includeSelf=*/false); + while ( smIt->more() ) + { + SMESH_subMesh* sm = smIt->next(); + _EdgesOnShape* subEOS = eos.GetData().GetShapeEdges( sm->GetId() ); + if ( !subEOS->_sWOL.IsNull() ) continue; + if ( std::find( eosC1.begin(), eosC1.end(), subEOS ) != eosC1.end() ) continue; + + putOnOffsetSurface( *subEOS, infStep, eosC1, smooStep, _LayerEdge::UPD_NORMAL_CONV ); + } + cnvFace->_normalsFixedOnBorders = true; + } } //================================================================================ @@ -5439,6 +5576,109 @@ Handle(Geom_Curve) _Smoother1D::CurveForSmooth( const TopoDS_Edge& E, return Handle(Geom_Curve)(); } +//================================================================================ +/*! + * \brief Smooth edges on EDGE + */ +//================================================================================ + +bool _Smoother1D::Perform(_SolidData& data, + Handle(ShapeAnalysis_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper ) +{ + if ( _leParams.empty() || ( !isAnalytic() && _offPoints.empty() )) + prepare( data ); + + findEdgesToSmooth(); + if ( isAnalytic() ) + return smoothAnalyticEdge( data, surface, F, helper ); + else + return smoothComplexEdge ( data, surface, F, helper ); +} + +//================================================================================ +/*! + * \brief Find edges to smooth + */ +//================================================================================ + +void _Smoother1D::findEdgesToSmooth() +{ + _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + if ( leOnV[iEnd]->Is( _LayerEdge::NORMAL_UPDATED )) + _leOnV[iEnd]._cosin = Abs( _edgeDir[iEnd].Normalized() * leOnV[iEnd]->_normal ); + + _eToSmooth[0].first = _eToSmooth[0].second = 0; + + for ( size_t i = 0; i < _eos.size(); ++i ) + { + if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) + { + if ( needSmoothing( _leOnV[0]._cosin, + _eos[i]->_len * leOnV[0]->_lenFactor, _curveLen * _leParams[i] ) || + isToSmooth( i ) + ) + _eos[i]->Set( _LayerEdge::TO_SMOOTH ); + else + break; + } + _eToSmooth[0].second = i+1; + } + + _eToSmooth[1].first = _eToSmooth[1].second = _eos.size(); + + for ( int i = _eos.size() - 1; i >= _eToSmooth[0].second; --i ) + { + if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) + { + if ( needSmoothing( _leOnV[1]._cosin, + _eos[i]->_len * leOnV[1]->_lenFactor, _curveLen * ( 1.-_leParams[i] )) || + isToSmooth( i )) + _eos[i]->Set( _LayerEdge::TO_SMOOTH ); + else + break; + } + _eToSmooth[1].first = i; + } +} + +//================================================================================ +/*! + * \brief Check if iE-th _LayerEdge needs smoothing + */ +//================================================================================ + +bool _Smoother1D::isToSmooth( int iE ) +{ + SMESH_NodeXYZ pi( _eos[iE]->_nodes[0] ); + SMESH_NodeXYZ p0( _eos[iE]->_2neibors->srcNode(0) ); + SMESH_NodeXYZ p1( _eos[iE]->_2neibors->srcNode(1) ); + gp_XYZ seg0 = pi - p0; + gp_XYZ seg1 = p1 - pi; + gp_XYZ tangent = seg0 + seg1; + double tangentLen = tangent.Modulus(); + double segMinLen = Min( seg0.Modulus(), seg1.Modulus() ); + if ( tangentLen < std::numeric_limits::min() ) + return false; + tangent /= tangentLen; + + for ( size_t i = 0; i < _eos[iE]->_neibors.size(); ++i ) + { + _LayerEdge* ne = _eos[iE]->_neibors[i]; + if ( !ne->Is( _LayerEdge::TO_SMOOTH ) || + ne->_nodes.size() < 2 || + ne->_nodes[0]->GetPosition()->GetDim() != 2 ) + continue; + gp_XYZ edgeVec = SMESH_NodeXYZ( ne->_nodes.back() ) - SMESH_NodeXYZ( ne->_nodes[0] ); + double proj = edgeVec * tangent; + if ( needSmoothing( 1., proj, segMinLen )) + return true; + } + return false; +} + //================================================================================ /*! * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE @@ -5452,80 +5692,121 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, { if ( !isAnalytic() ) return false; - const size_t iFrom = 0, iTo = _eos._edges.size(); + size_t iFrom = 0, iTo = _eos._edges.size(); if ( _anaCurve->IsKind( STANDARD_TYPE( Geom_Line ))) { if ( F.IsNull() ) // 3D { - SMESH_TNodeXYZ p0 ( _eos._edges[iFrom]->_2neibors->tgtNode(0) ); - SMESH_TNodeXYZ p1 ( _eos._edges[iTo-1]->_2neibors->tgtNode(1) ); SMESH_TNodeXYZ pSrc0( _eos._edges[iFrom]->_2neibors->srcNode(0) ); SMESH_TNodeXYZ pSrc1( _eos._edges[iTo-1]->_2neibors->srcNode(1) ); - gp_XYZ newPos, lineDir = pSrc1 - pSrc0; - _LayerEdge* vLE0 = _eos._edges[iFrom]->_2neibors->_edges[0]; - _LayerEdge* vLE1 = _eos._edges[iTo-1]->_2neibors->_edges[1]; - bool shiftOnly = ( vLE0->Is( _LayerEdge::NORMAL_UPDATED ) || - vLE0->Is( _LayerEdge::BLOCKED ) || - vLE1->Is( _LayerEdge::NORMAL_UPDATED ) || - vLE1->Is( _LayerEdge::BLOCKED )); - for ( size_t i = iFrom; i < iTo; ++i ) - { - _LayerEdge* edge = _eos._edges[i]; - SMDS_MeshNode* tgtNode = const_cast( edge->_nodes.back() ); - newPos = p0 * ( 1. - _leParams[i] ) + p1 * _leParams[i]; - - if ( shiftOnly || edge->Is( _LayerEdge::NORMAL_UPDATED )) + //const gp_XYZ lineDir = pSrc1 - pSrc0; + //_LayerEdge* vLE0 = getLEdgeOnV( 0 ); + //_LayerEdge* vLE1 = getLEdgeOnV( 1 ); + // bool shiftOnly = ( vLE0->Is( _LayerEdge::NORMAL_UPDATED ) || + // vLE0->Is( _LayerEdge::BLOCKED ) || + // vLE1->Is( _LayerEdge::NORMAL_UPDATED ) || + // vLE1->Is( _LayerEdge::BLOCKED )); + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + { + iFrom = _eToSmooth[ iEnd ].first, iTo = _eToSmooth[ iEnd ].second; + if ( iFrom >= iTo ) continue; + SMESH_TNodeXYZ p0( _eos[iFrom]->_2neibors->tgtNode(0) ); + SMESH_TNodeXYZ p1( _eos[iTo-1]->_2neibors->tgtNode(1) ); + double param0 = ( iFrom == 0 ) ? 0. : _leParams[ iFrom-1 ]; + double param1 = _leParams[ iTo ]; + for ( size_t i = iFrom; i < iTo; ++i ) { - gp_XYZ curPos = SMESH_TNodeXYZ ( tgtNode ); - double shift = ( lineDir * ( newPos - pSrc0 ) - - lineDir * ( curPos - pSrc0 )); - newPos = curPos + lineDir * shift / lineDir.SquareModulus(); - } - if ( edge->Is( _LayerEdge::BLOCKED )) - { - SMESH_TNodeXYZ pSrc( edge->_nodes[0] ); - double curThick = pSrc.SquareDistance( tgtNode ); - double newThink = ( pSrc - newPos ).SquareModulus(); - if ( newThink > curThick ) - continue; + _LayerEdge* edge = _eos[i]; + SMDS_MeshNode* tgtNode = const_cast( edge->_nodes.back() ); + double param = ( _leParams[i] - param0 ) / ( param1 - param0 ); + gp_XYZ newPos = p0 * ( 1. - param ) + p1 * param; + + // if ( shiftOnly || edge->Is( _LayerEdge::NORMAL_UPDATED )) + // { + // gp_XYZ curPos = SMESH_TNodeXYZ ( tgtNode ); + // double shift = ( lineDir * ( newPos - pSrc0 ) - + // lineDir * ( curPos - pSrc0 )); + // newPos = curPos + lineDir * shift / lineDir.SquareModulus(); + // } + if ( edge->Is( _LayerEdge::BLOCKED )) + { + SMESH_TNodeXYZ pSrc( edge->_nodes[0] ); + double curThick = pSrc.SquareDistance( tgtNode ); + double newThink = ( pSrc - newPos ).SquareModulus(); + if ( newThink > curThick ) + continue; + } + edge->_pos.back() = newPos; + tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); + dumpMove( tgtNode ); } - edge->_pos.back() = newPos; - tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); - dumpMove( tgtNode ); } } else // 2D { - _LayerEdge* e0 = getLEdgeOnV( 0 ); - _LayerEdge* e1 = getLEdgeOnV( 1 ); - gp_XY uv0 = e0->LastUV( F, *data.GetShapeEdges( e0 )); - gp_XY uv1 = e1->LastUV( F, *data.GetShapeEdges( e1 )); - if ( e0->_nodes.back() == e1->_nodes.back() ) // closed edge + _LayerEdge* eV0 = getLEdgeOnV( 0 ); + _LayerEdge* eV1 = getLEdgeOnV( 1 ); + gp_XY uvV0 = eV0->LastUV( F, *data.GetShapeEdges( eV0 )); + gp_XY uvV1 = eV1->LastUV( F, *data.GetShapeEdges( eV1 )); + if ( eV0->_nodes.back() == eV1->_nodes.back() ) // 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 ); + uvV1.SetCoord( iPeriodic, helper.GetOtherParam( uvV1.Coord( iPeriodic ))); + if ( uvV0.Coord( iPeriodic ) > uvV1.Coord( iPeriodic )) + std::swap( uvV0, uvV1 ); } } - const gp_XY rangeUV = uv1 - uv0; - for ( size_t i = iFrom; i < iTo; ++i ) - { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; - gp_XY newUV = uv0 + _leParams[i] * rangeUV; - _eos._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + { + iFrom = _eToSmooth[ iEnd ].first, iTo = _eToSmooth[ iEnd ].second; + if ( iFrom >= iTo ) continue; + _LayerEdge* e0 = _eos[iFrom]->_2neibors->_edges[0]; + _LayerEdge* e1 = _eos[iTo-1]->_2neibors->_edges[1]; + gp_XY uv0 = ( e0 == eV0 ) ? uvV0 : e0->LastUV( F, _eos ); + gp_XY uv1 = ( e1 == eV1 ) ? uvV1 : e1->LastUV( F, _eos ); + double param0 = ( iFrom == 0 ) ? 0. : _leParams[ iFrom-1 ]; + double param1 = _leParams[ iTo ]; + gp_XY rangeUV = uv1 - uv0; + for ( size_t i = iFrom; i < iTo; ++i ) + { + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + double param = ( _leParams[i] - param0 ) / ( param1 - param0 ); + gp_XY newUV = uv0 + param * rangeUV; - gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() ); - SMDS_MeshNode* tgtNode = const_cast( _eos._edges[i]->_nodes.back() ); - tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); - dumpMove( tgtNode ); + gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() ); + SMDS_MeshNode* tgtNode = const_cast( _eos[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() ); + SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); + pos->SetUParameter( newUV.X() ); + pos->SetVParameter( newUV.Y() ); + + gp_XYZ newUV0( newUV.X(), newUV.Y(), 0 ); + + if ( !_eos[i]->Is( _LayerEdge::SMOOTHED )) + { + _eos[i]->Set( _LayerEdge::SMOOTHED ); // to check in refine() (IPAL54237) + if ( _eos[i]->_pos.size() > 2 ) + { + // modify previous positions to make _LayerEdge less sharply bent + vector& uvVec = _eos[i]->_pos; + const gp_XYZ uvShift = newUV0 - uvVec.back(); + const double len2 = ( uvVec.back() - uvVec[ 0 ] ).SquareModulus(); + int iPrev = uvVec.size() - 2; + while ( iPrev > 0 ) + { + double r = ( uvVec[ iPrev ] - uvVec[0] ).SquareModulus() / len2; + uvVec[ iPrev ] += uvShift * r; + --iPrev; + } + } + } + _eos[i]->_pos.back() = newUV0; + } } } return true; @@ -5564,9 +5845,10 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, if ( uLast < 0 ) uLast += 2 * M_PI; - for ( size_t i = iFrom; i < iTo; ++i ) + for ( size_t i = 0; i < _eos.size(); ++i ) { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue; double u = uLast * _leParams[i]; gp_Pnt p = ElCLib::Value( u, newCirc ); _eos._edges[i]->_pos.back() = p.XYZ(); @@ -5598,9 +5880,10 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, gp_Ax2d axis( center, vec0 ); gp_Circ2d circ( axis, radius ); - for ( size_t i = iFrom; i < iTo; ++i ) + for ( size_t i = 0; i < _eos.size(); ++i ) { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue; double newU = uLast * _leParams[i]; gp_Pnt2d newUV = ElCLib::Value( newU, circ ); _eos._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); @@ -5613,6 +5896,8 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); pos->SetUParameter( newUV.X() ); pos->SetVParameter( newUV.Y() ); + + _eos[i]->Set( _LayerEdge::SMOOTHED ); // to check in refine() (IPAL54237) } } return true; @@ -5635,7 +5920,9 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, if ( _offPoints.empty() ) return false; + // ---------------------------------------------- // move _offPoints along normals of _LayerEdge's + // ---------------------------------------------- _LayerEdge* e[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; if ( e[0]->Is( _LayerEdge::NORMAL_UPDATED )) @@ -5679,18 +5966,23 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, } } + // ----------------------------------------------------------------- // project tgt nodes of extreme _LayerEdge's to the offset segments + // ----------------------------------------------------------------- - if ( e[0]->Is( _LayerEdge::NORMAL_UPDATED )) _iSeg[0] = 0; - if ( e[1]->Is( _LayerEdge::NORMAL_UPDATED )) _iSeg[1] = _offPoints.size()-2; + const int updatedOrBlocked = _LayerEdge::NORMAL_UPDATED | _LayerEdge::BLOCKED; + if ( e[0]->Is( updatedOrBlocked )) _iSeg[0] = 0; + if ( e[1]->Is( updatedOrBlocked )) _iSeg[1] = _offPoints.size()-2; gp_Pnt pExtreme[2], pProj[2]; + bool isProjected[2]; for ( int is2nd = 0; is2nd < 2; ++is2nd ) { pExtreme[ is2nd ] = SMESH_TNodeXYZ( e[is2nd]->_nodes.back() ); int i = _iSeg[ is2nd ]; int di = is2nd ? -1 : +1; - bool projected = false; + bool & projected = isProjected[ is2nd ]; + projected = false; double uOnSeg, distMin = Precision::Infinite(), dist, distPrev = 0; int nbWorse = 0; do { @@ -5737,42 +6029,52 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, gp_Vec vDiv0( pExtreme[0], pProj[0] ); gp_Vec vDiv1( pExtreme[1], pProj[1] ); double d0 = vDiv0.Magnitude(); - double d1 = vDiv1.Magnitude(); - if ( e[0]->_normal * vDiv0.XYZ() < 0 ) e[0]->_len += d0; - else e[0]->_len -= d0; - if ( e[1]->_normal * vDiv1.XYZ() < 0 ) e[1]->_len += d1; - else e[1]->_len -= d1; + double d1 = isProjected[1] ? vDiv1.Magnitude() : 0; + if ( e[0]->Is( _LayerEdge::BLOCKED )) { + if ( e[0]->_normal * vDiv0.XYZ() < 0 ) e[0]->_len += d0; + else e[0]->_len -= d0; + } + if ( e[1]->Is( _LayerEdge::BLOCKED )) { + if ( e[1]->_normal * vDiv1.XYZ() < 0 ) e[1]->_len += d1; + else e[1]->_len -= d1; + } + // --------------------------------------------------------------------------------- // compute normalized length of the offset segments located between the projections + // --------------------------------------------------------------------------------- + + // temporary replace extreme _offPoints by pExtreme + gp_XYZ opXYZ[2] = { _offPoints[ _iSeg[0] ]._xyz, + _offPoints[ _iSeg[1]+1 ]._xyz }; + _offPoints[ _iSeg[0] ]._xyz = pExtreme[0].XYZ(); + _offPoints[ _iSeg[1]+ 1]._xyz = pExtreme[1].XYZ(); size_t iSeg = 0, nbSeg = _iSeg[1] - _iSeg[0] + 1; vector< double > len( nbSeg + 1 ); len[ iSeg++ ] = 0; - len[ iSeg++ ] = pProj[ 0 ].Distance( _offPoints[ _iSeg[0]+1 ]._xyz )/* * e[0]->_lenFactor*/; + len[ iSeg++ ] = pProj[ 0 ].Distance( _offPoints[ _iSeg[0]+1 ]._xyz ); for ( size_t i = _iSeg[0]+1; i <= _iSeg[1]; ++i, ++iSeg ) { len[ iSeg ] = len[ iSeg-1 ] + _offPoints[i].Distance( _offPoints[i+1] ); } - len[ nbSeg ] -= pProj[ 1 ].Distance( _offPoints[ _iSeg[1]+1 ]._xyz )/* * e[1]->_lenFactor*/; + // if ( isProjected[ 1 ]) + // len[ nbSeg ] -= pProj[ 1 ].Distance( _offPoints[ _iSeg[1]+1 ]._xyz ); + // else + // len[ nbSeg ] += pExtreme[ 1 ].Distance( _offPoints[ _iSeg[1]+1 ]._xyz ); - // d0 *= e[0]->_lenFactor; - // d1 *= e[1]->_lenFactor; double fullLen = len.back() - d0 - d1; for ( iSeg = 0; iSeg < len.size(); ++iSeg ) len[iSeg] = ( len[iSeg] - d0 ) / fullLen; - // temporary replace extreme _offPoints by pExtreme - gp_XYZ op[2] = { _offPoints[ _iSeg[0] ]._xyz, - _offPoints[ _iSeg[1]+1 ]._xyz }; - _offPoints[ _iSeg[0] ]._xyz = pExtreme[0].XYZ(); - _offPoints[ _iSeg[1]+ 1]._xyz = pExtreme[1].XYZ(); - + // ------------------------------------------------------------- // distribute tgt nodes of _LayerEdge's between the projections + // ------------------------------------------------------------- iSeg = 0; - for ( size_t i = 0; i < _eos._edges.size(); ++i ) + for ( size_t i = 0; i < _eos.size(); ++i ) { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue; while ( iSeg+2 < len.size() && _leParams[i] > len[ iSeg+1 ] ) iSeg++; double r = ( _leParams[i] - len[ iSeg ]) / ( len[ iSeg+1 ] - len[ iSeg ]); @@ -5781,23 +6083,23 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, if ( surface.IsNull() ) { - _eos._edges[i]->_pos.back() = p; + _eos[i]->_pos.back() = p; } else // project a new node position to a FACE { - gp_Pnt2d uv ( _eos._edges[i]->_pos.back().X(), _eos._edges[i]->_pos.back().Y() ); + gp_Pnt2d uv ( _eos[i]->_pos.back().X(), _eos[i]->_pos.back().Y() ); gp_Pnt2d uv2( surface->NextValueOfUV( uv, p, fTol )); p = surface->Value( uv2 ).XYZ(); - _eos._edges[i]->_pos.back().SetCoord( uv2.X(), uv2.Y(), 0 ); + _eos[i]->_pos.back().SetCoord( uv2.X(), uv2.Y(), 0 ); } - SMDS_MeshNode* tgtNode = const_cast( _eos._edges[i]->_nodes.back() ); + SMDS_MeshNode* tgtNode = const_cast( _eos[i]->_nodes.back() ); tgtNode->setXYZ( p.X(), p.Y(), p.Z() ); dumpMove( tgtNode ); } - _offPoints[ _iSeg[0] ]._xyz = op[0]; - _offPoints[ _iSeg[1]+1 ]._xyz = op[1]; + _offPoints[ _iSeg[0] ]._xyz = opXYZ[0]; + _offPoints[ _iSeg[1]+1 ]._xyz = opXYZ[1]; return true; } @@ -5832,15 +6134,27 @@ void _Smoother1D::prepare(_SolidData& data) double fullLen = _leParams.back() + pPrev.Distance( SMESH_TNodeXYZ( getLEdgeOnV(1)->_nodes[0])); for ( size_t i = 0; i < _leParams.size()-1; ++i ) _leParams[i] = _leParams[i+1] / fullLen; + _leParams.back() = 1.; } + _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; + + // get cosin to use in findEdgesToSmooth() + _edgeDir[0] = getEdgeDir( E, leOnV[0]->_nodes[0], data.GetHelper() ); + _edgeDir[1] = getEdgeDir( E, leOnV[1]->_nodes[0], data.GetHelper() ); + _leOnV[0]._cosin = Abs( leOnV[0]->_cosin ); + _leOnV[1]._cosin = Abs( leOnV[1]->_cosin ); + if ( _eos._sWOL.IsNull() ) // 3D + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + _leOnV[iEnd]._cosin = Abs( _edgeDir[iEnd].Normalized() * leOnV[iEnd]->_normal ); + if ( isAnalytic() ) return; // divide E to have offset segments with low deflection BRepAdaptor_Curve c3dAdaptor( E ); - const double curDeflect = 0.1; //0.3; // 0.01; // Curvature deflection - const double angDeflect = 0.1; //0.2; // 0.09; // Angular deflection + const double curDeflect = 0.1; //0.01; // Curvature deflection == |p1p2]*sin(p1p2,p1pM) + const double angDeflect = 0.1; //0.09; // Angular deflection == sin(p1pM,pMp2) GCPnts_TangentialDeflection discret(c3dAdaptor, angDeflect, curDeflect); if ( discret.NbPoints() <= 2 ) { @@ -5850,17 +6164,40 @@ void _Smoother1D::prepare(_SolidData& data) const double u0 = c3dAdaptor.FirstParameter(); gp_Pnt p; gp_Vec tangent; - _offPoints.resize( discret.NbPoints() ); - for ( size_t i = 0; i < _offPoints.size(); i++ ) + if ( discret.NbPoints() >= (int) _eos.size() + 2 ) { - double u = discret.Parameter( i+1 ); - c3dAdaptor.D1( u, p, tangent ); - _offPoints[i]._xyz = p.XYZ(); - _offPoints[i]._edgeDir = tangent.XYZ(); - _offPoints[i]._param = GCPnts_AbscissaPoint::Length( c3dAdaptor, u0, u ) / _curveLen; + _offPoints.resize( discret.NbPoints() ); + for ( size_t i = 0; i < _offPoints.size(); i++ ) + { + double u = discret.Parameter( i+1 ); + c3dAdaptor.D1( u, p, tangent ); + _offPoints[i]._xyz = p.XYZ(); + _offPoints[i]._edgeDir = tangent.XYZ(); + _offPoints[i]._param = GCPnts_AbscissaPoint::Length( c3dAdaptor, u0, u ) / _curveLen; + } } + else + { + std::vector< double > params( _eos.size() + 2 ); - _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; + params[0] = data.GetHelper().GetNodeU( E, leOnV[0]->_nodes[0] ); + params.back() = data.GetHelper().GetNodeU( E, leOnV[1]->_nodes[0] ); + for ( size_t i = 0; i < _eos.size(); i++ ) + params[i+1] = data.GetHelper().GetNodeU( E, _eos[i]->_nodes[0] ); + + if ( params[1] > params[ _eos.size() ] ) + std::reverse( params.begin() + 1, params.end() - 1 ); + + _offPoints.resize( _eos.size() + 2 ); + for ( size_t i = 0; i < _offPoints.size(); i++ ) + { + const double u = params[i]; + c3dAdaptor.D1( u, p, tangent ); + _offPoints[i]._xyz = p.XYZ(); + _offPoints[i]._edgeDir = tangent.XYZ(); + _offPoints[i]._param = GCPnts_AbscissaPoint::Length( c3dAdaptor, u0, u ) / _curveLen; + } + } // set _2edges _offPoints [0]._2edges.set( &_leOnV[0], &_leOnV[0], 0.5, 0.5 ); @@ -5899,11 +6236,14 @@ void _Smoother1D::prepare(_SolidData& data) int iLBO = _offPoints.size() - 2; // last but one - _edgeDir[0] = getEdgeDir( E, leOnV[0]->_nodes[0], data.GetHelper() ); - _edgeDir[1] = getEdgeDir( E, leOnV[1]->_nodes[0], data.GetHelper() ); - - _leOnV[ 0 ]._normal = getNormalNormal( leOnV[0]->_normal, _edgeDir[0] ); - _leOnV[ 1 ]._normal = getNormalNormal( leOnV[1]->_normal, _edgeDir[1] ); + if ( leOnV[ 0 ]->Is( _LayerEdge::MULTI_NORMAL )) + _leOnV[ 0 ]._normal = getNormalNormal( _eos._edges[1]->_normal, _edgeDir[0] ); + else + _leOnV[ 0 ]._normal = getNormalNormal( leOnV[0]->_normal, _edgeDir[0] ); + if ( leOnV[ 1 ]->Is( _LayerEdge::MULTI_NORMAL )) + _leOnV[ 1 ]._normal = getNormalNormal( _eos._edges.back()->_normal, _edgeDir[1] ); + else + _leOnV[ 1 ]._normal = getNormalNormal( leOnV[1]->_normal, _edgeDir[1] ); _leOnV[ 0 ]._len = 0; _leOnV[ 1 ]._len = 0; _leOnV[ 0 ]._lenFactor = _offPoints[1 ]._2edges._edges[1]->_lenFactor; @@ -5937,7 +6277,7 @@ void _Smoother1D::prepare(_SolidData& data) //================================================================================ /*! - * \brief set _normal of _leOnV[is2nd] to be normal to the EDGE + * \brief return _normal of _leOnV[is2nd] normal to the EDGE */ //================================================================================ @@ -5948,9 +6288,36 @@ gp_XYZ _Smoother1D::getNormalNormal( const gp_XYZ & normal, gp_XYZ norm = edgeDir ^ cross; double size = norm.Modulus(); + // if ( size == 0 ) // MULTI_NORMAL _LayerEdge + // return gp_XYZ( 1e-100, 1e-100, 1e-100 ); + return norm / size; } +//================================================================================ +/*! + * \brief Writes a script creating a mesh composed of _offPoints + */ +//================================================================================ + +void _Smoother1D::offPointsToPython() const +{ + const char* fname = "/tmp/offPoints.py"; + cout << "execfile('"<_edges.size(); ++i ) - { - eos->_edges[i]->SetSmooLen( Precision::Infinite() ); - } - SMESH_subMeshIteratorPtr smIt = eos->_subMesh->getDependsOnIterator(/*includeSelf=*/false); - while ( smIt->more() ) // loop on sub-shapes of the FACE - { - _EdgesOnShape* eoe = GetShapeEdges( smIt->next()->GetId() ); - if ( !eoe ) continue; - - vector<_LayerEdge*>& eE = eoe->_edges; - for ( size_t iE = 0; iE < eE.size(); ++iE ) // loop on _LayerEdge's on EDGE or VERTEX - { - if ( eE[iE]->_cosin <= theMinSmoothCosin ) - continue; + // for ( size_t i = 0; i < eos->_edges.size(); ++i ) + // { + // eos->_edges[i]->SetSmooLen( Precision::Infinite() ); + // } + // SMESH_subMeshIteratorPtr smIt = eos->_subMesh->getDependsOnIterator(/*includeSelf=*/false); + // while ( smIt->more() ) // loop on sub-shapes of the FACE + // { + // _EdgesOnShape* eoe = GetShapeEdges( smIt->next()->GetId() ); + // if ( !eoe ) continue; - SMDS_ElemIteratorPtr segIt = eE[iE]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge); - while ( segIt->more() ) - { - const SMDS_MeshElement* seg = segIt->next(); - if ( !eos->_subMesh->DependsOn( seg->getshapeId() )) - continue; - if ( seg->GetNode(0) != eE[iE]->_nodes[0] ) - continue; // not to check a seg twice - for ( size_t iN = 0; iN < eE[iE]->_neibors.size(); ++iN ) - { - _LayerEdge* eN = eE[iE]->_neibors[iN]; - if ( eN->_nodes[0]->getshapeId() != eos->_shapeID ) - continue; - double dist = SMESH_MeshAlgos::GetDistance( seg, SMESH_TNodeXYZ( eN->_nodes[0] )); - double smooLen = getSmoothingThickness( eE[iE]->_cosin, dist ); - eN->SetSmooLen( Min( smooLen, eN->GetSmooLen() )); - eN->Set( _LayerEdge::NEAR_BOUNDARY ); - } - } - } - } + // vector<_LayerEdge*>& eE = eoe->_edges; + // for ( size_t iE = 0; iE < eE.size(); ++iE ) // loop on _LayerEdge's on EDGE or VERTEX + // { + // if ( eE[iE]->_cosin <= theMinSmoothCosin ) + // continue; + + // SMDS_ElemIteratorPtr segIt = eE[iE]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge); + // while ( segIt->more() ) + // { + // const SMDS_MeshElement* seg = segIt->next(); + // if ( !eos->_subMesh->DependsOn( seg->getshapeId() )) + // continue; + // if ( seg->GetNode(0) != eE[iE]->_nodes[0] ) + // continue; // not to check a seg twice + // for ( size_t iN = 0; iN < eE[iE]->_neibors.size(); ++iN ) + // { + // _LayerEdge* eN = eE[iE]->_neibors[iN]; + // if ( eN->_nodes[0]->getshapeId() != eos->_shapeID ) + // continue; + // double dist = SMESH_MeshAlgos::GetDistance( seg, SMESH_TNodeXYZ( eN->_nodes[0] )); + // double smooLen = getSmoothingThickness( eE[iE]->_cosin, dist ); + // eN->SetSmooLen( Min( smooLen, eN->GetSmooLen() )); + // eN->Set( _LayerEdge::NEAR_BOUNDARY ); + // } + // } + // } + // } } // if ( eos->ShapeType() == TopAbs_FACE ) for ( size_t i = 0; i < eos->_edges.size(); ++i ) @@ -6132,6 +6499,7 @@ void _SolidData::PrepareEdgesToSmoothOnFace( _EdgesOnShape* eos, bool substitute avgLen /= edge->_simplices.size(); if (( edge->_curvature = _Curvature::New( avgNormProj, avgLen ))) { + edge->Set( _LayerEdge::SMOOTHED_C1 ); isCurved = true; SMDS_FacePosition* fPos = dynamic_cast( edge->_nodes[0]->GetPosition() ); if ( !fPos ) @@ -6227,7 +6595,7 @@ void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelp if ( eI->_nodes[0]->GetID() < eN->_nodes[0]->GetID() ) // treat this pair once { _EdgesOnShape* eosN = data.GetShapeEdges( eN ); - limitMaxLenByCurvature( eI, eN, eosI, *eosN, helper ); + limitMaxLenByCurvature( eI, eN, eosI, *eosN, eosI._hyp.ToSmooth() ); } } } @@ -6241,7 +6609,7 @@ void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelp for ( size_t i = 1; i < eosI._edges.size(); ++i ) { _LayerEdge* eI = eosI._edges[i]; - limitMaxLenByCurvature( eI, e0, eosI, eosI, helper ); + limitMaxLenByCurvature( eI, e0, eosI, eosI, eosI._hyp.ToSmooth() ); e0 = eI; } } @@ -6254,12 +6622,17 @@ void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelp */ //================================================================================ -void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge* e1, - _LayerEdge* e2, - _EdgesOnShape& eos1, - _EdgesOnShape& eos2, - SMESH_MesherHelper& helper ) +void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge* e1, + _LayerEdge* e2, + _EdgesOnShape& eos1, + _EdgesOnShape& eos2, + const bool isSmoothable ) { + if (( e1->_nodes[0]->GetPosition()->GetDim() != + e2->_nodes[0]->GetPosition()->GetDim() ) && + ( e1->_cosin < 0.75 )) + return; // angle > 90 deg at e1 + gp_XYZ plnNorm = e1->_normal ^ e2->_normal; double norSize = plnNorm.SquareModulus(); if ( norSize < std::numeric_limits::min() ) @@ -6279,9 +6652,10 @@ void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge* e1, double ovl = ( u1 * e1->_normal * dir12 - u2 * e2->_normal * dir12 ) / dir12.SquareModulus(); if ( ovl > theSmoothThickToElemSizeRatio ) - { - e1->_maxLen = Min( e1->_maxLen, 0.75 * u1 / e1->_lenFactor ); - e2->_maxLen = Min( e2->_maxLen, 0.75 * u2 / e2->_lenFactor ); + { + const double coef = 0.75; + e1->SetMaxLen( Min( e1->_maxLen, coef * u1 / e1->_lenFactor )); + e2->SetMaxLen( Min( e2->_maxLen, coef * u2 / e2->_lenFactor )); } } } @@ -6456,7 +6830,7 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper& // else // { // double shortLen = 0.75 * ( Min( dist1, dist2 ) / edge->_lenFactor ); - // edge->_maxLen = Min( shortLen, edge->_maxLen ); + // edge->SetMaxLen( Min( shortLen, edge->_maxLen )); // } } @@ -6486,6 +6860,80 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper& } } +//================================================================================ +/*! + * \brief Find _LayerEdge's located on boundary of a convex FACE whose normal + * will be updated at each inflation step + */ +//================================================================================ + +void _ViscousBuilder::findEdgesToUpdateNormalNearConvexFace( _ConvexFace & convFace, + _SolidData& data, + SMESH_MesherHelper& helper ) +{ + const TGeomID convFaceID = getMeshDS()->ShapeToIndex( convFace._face ); + const double preci = BRep_Tool::Tolerance( convFace._face ); + Handle(ShapeAnalysis_Surface) surface = helper.GetSurface( convFace._face ); + + bool edgesToUpdateFound = false; + + map< TGeomID, _EdgesOnShape* >::iterator id2eos = convFace._subIdToEOS.begin(); + for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos ) + { + _EdgesOnShape& eos = * id2eos->second; + if ( !eos._sWOL.IsNull() ) continue; + if ( !eos._hyp.ToSmooth() ) continue; + for ( size_t i = 0; i < eos._edges.size(); ++i ) + { + _LayerEdge* ledge = eos._edges[ i ]; + if ( ledge->Is( _LayerEdge::UPD_NORMAL_CONV )) continue; // already checked + if ( ledge->Is( _LayerEdge::MULTI_NORMAL )) continue; // not inflatable + + gp_XYZ tgtPos = ( SMESH_NodeXYZ( ledge->_nodes[0] ) + + ledge->_normal * ledge->_lenFactor * ledge->_maxLen ); + + // the normal must be updated if distance from tgtPos to surface is less than + // target thickness + + // find an initial UV for search of a projection of tgtPos to surface + const SMDS_MeshNode* nodeInFace = 0; + SMDS_ElemIteratorPtr fIt = ledge->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() && !nodeInFace ) + { + const SMDS_MeshElement* f = fIt->next(); + if ( convFaceID != f->getshapeId() ) continue; + + SMDS_ElemIteratorPtr nIt = f->nodesIterator(); + while ( nIt->more() && !nodeInFace ) + { + const SMDS_MeshElement* n = nIt->next(); + if ( n->getshapeId() == convFaceID ) + nodeInFace = static_cast< const SMDS_MeshNode* >( n ); + } + } + if ( !nodeInFace ) + continue; + gp_XY uv = helper.GetNodeUV( convFace._face, nodeInFace ); + + // projection + surface->NextValueOfUV( uv, tgtPos, preci ); + double dist = surface->Gap(); + if ( dist < 0.95 * ledge->_maxLen ) + { + ledge->Set( _LayerEdge::UPD_NORMAL_CONV ); + if ( !ledge->_curvature ) ledge->_curvature = new _Curvature; + ledge->_curvature->_uv.SetCoord( uv.X(), uv.Y() ); + edgesToUpdateFound = true; + } + } + } + + if ( !convFace._isTooCurved && edgesToUpdateFound ) + { + data._convexFaces.insert( make_pair( convFaceID, convFace )).first->second; + } +} + //================================================================================ /*! * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with @@ -6593,8 +7041,8 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, // compute new _normals for ( size_t i = 0; i < intEdgesDist.size(); ++i ) { - _LayerEdge* edge2 = intEdgesDist[i].first; - double distWgt = edge1->_len / intEdgesDist[i].second; + _LayerEdge* edge2 = intEdgesDist[i].first; + double distWgt = edge1->_len / intEdgesDist[i].second; // if ( edge1->Is( _LayerEdge::BLOCKED ) && // edge2->Is( _LayerEdge::BLOCKED )) continue; if ( edge2->Is( _LayerEdge::MARKED )) continue; @@ -6632,12 +7080,17 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, e2neIt = edge2newEdge.insert( make_pair( edge1, zeroEdge )).first; e2neIt->second._normal += distWgt * newNormal; e2neIt->second._cosin = newCos; - e2neIt->second._maxLen = 0.7 * minIntDist / edge1->_lenFactor; + e2neIt->second.SetMaxLen( 0.7 * minIntDist / edge1->_lenFactor ); if ( iter > 0 && sgn1 * sgn2 < 0 && edge1->_cosin < 0 ) e2neIt->second._normal += dir2; + e2neIt = edge2newEdge.insert( make_pair( edge2, zeroEdge )).first; e2neIt->second._normal += distWgt * newNormal; - e2neIt->second._cosin = edge2->_cosin; + if ( Precision::IsInfinite( zeroEdge._maxLen )) + { + e2neIt->second._cosin = edge2->_cosin; + e2neIt->second.SetMaxLen( 1.3 * minIntDist / edge1->_lenFactor ); + } if ( iter > 0 && sgn1 * sgn2 < 0 && edge2->_cosin < 0 ) e2neIt->second._normal += dir1; } @@ -6654,9 +7107,10 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, for ( e2neIt = edge2newEdge.begin(); e2neIt != edge2newEdge.end(); ++e2neIt ) { _LayerEdge* edge = e2neIt->first; - if ( edge->Is( _LayerEdge::BLOCKED )) continue; _LayerEdge& newEdge = e2neIt->second; _EdgesOnShape* eos = data.GetShapeEdges( edge ); + if ( edge->Is( _LayerEdge::BLOCKED && newEdge._maxLen > edge->_len )) + continue; // Check if a new _normal is OK: newEdge._normal.Normalize(); @@ -6665,7 +7119,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, if ( newEdge._maxLen < edge->_len && iter > 0 ) // limit _maxLen { edge->InvalidateStep( stepNb + 1, *eos, /*restoreLength=*/true ); - edge->_maxLen = newEdge._maxLen; + edge->SetMaxLen( newEdge._maxLen ); edge->SetNewLength( newEdge._maxLen, *eos, helper ); } continue; // the new _normal is bad @@ -6938,6 +7392,8 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData& data, for ( ; id2face != data._convexFaces.end(); ++id2face ) { _ConvexFace & convFace = (*id2face).second; + convFace._normalsFixedOnBorders = false; // to update at each inflation step + if ( convFace._normalsFixed ) continue; // already fixed if ( convFace.CheckPrisms() ) @@ -7300,6 +7756,59 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData& data, return true; } +//================================================================================ +/*! + * \brief Return max curvature of a FACE + */ +//================================================================================ + +double _ConvexFace::GetMaxCurvature( _SolidData& data, + _EdgesOnShape& eof, + BRepLProp_SLProps& surfProp, + SMESH_MesherHelper& helper) +{ + double maxCurvature = 0; + + TopoDS_Face F = TopoDS::Face( eof._shape ); + + const int nbTestPnt = 5; + const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. ); + SMESH_subMeshIteratorPtr smIt = eof._subMesh->getDependsOnIterator(/*includeSelf=*/true); + while ( smIt->more() ) + { + SMESH_subMesh* sm = smIt->next(); + const TGeomID subID = sm->GetId(); + + // find _LayerEdge's of a sub-shape + _EdgesOnShape* eos; + if (( eos = data.GetShapeEdges( subID ))) + this->_subIdToEOS.insert( make_pair( subID, eos )); + else + continue; + + // check concavity and curvature and limit data._stepSize + const double minCurvature = + 1. / ( eos->_hyp.GetTotalThickness() * ( 1 + theThickToIntersection )); + size_t iStep = Max( 1, eos->_edges.size() / nbTestPnt ); + for ( size_t i = 0; i < eos->_edges.size(); i += iStep ) + { + gp_XY uv = helper.GetNodeUV( F, eos->_edges[ i ]->_nodes[0] ); + surfProp.SetParameters( uv.X(), uv.Y() ); + if ( surfProp.IsCurvatureDefined() ) + { + double curvature = Max( surfProp.MaxCurvature() * oriFactor, + surfProp.MinCurvature() * oriFactor ); + maxCurvature = Max( maxCurvature, curvature ); + + if ( curvature > minCurvature ) + this->_isTooCurved = true; + } + } + } // loop on sub-shapes of the FACE + + return maxCurvature; +} + //================================================================================ /*! * \brief Finds a center of curvature of a surface at a _LayerEdge @@ -7508,7 +8017,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher& searcher, gp_XYZ _LayerEdge::PrevCheckPos( _EdgesOnShape* eos ) const { - size_t i = Is( NORMAL_UPDATED ) ? _pos.size()-2 : 0; + size_t i = Is( NORMAL_UPDATED ) && IsOnFace() ? _pos.size()-2 : 0; if ( !eos || eos->_sWOL.IsNull() ) return _pos[ i ]; @@ -7582,13 +8091,14 @@ gp_Ax1 _LayerEdge::LastSegment(double& segLen, _EdgesOnShape& eos) const //================================================================================ /*! - * \brief Return the last position of the target node on a FACE. + * \brief Return the last (or \a which) position of the target node on a FACE. * \param [in] F - the FACE this _LayerEdge is inflated along + * \param [in] which - index of position * \return gp_XY - result UV */ //================================================================================ -gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const +gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos, int which ) const { if ( F.IsSame( eos._sWOL )) // F is my FACE return gp_XY( _pos.back().X(), _pos.back().Y() ); @@ -7597,7 +8107,7 @@ gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const return gp_XY( 1e100, 1e100 ); // _sWOL is EDGE of F; _pos.back().X() is the last U on the EDGE - double f, l, u = _pos.back().X(); + double f, l, u = _pos[ which < 0 ? _pos.size()-1 : which ].X(); Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( TopoDS::Edge(eos._sWOL), F, f,l); if ( !C2d.IsNull() && f <= u && u <= l ) return C2d->Value( u ).XY(); @@ -7675,7 +8185,7 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment, * \param [in] eov - EOS of the VERTEX * \param [in] eos - EOS of the FACE * \param [in] step - inflation step - * \param [in,out] badSmooEdges - not untangled _LayerEdge's + * \param [in,out] badSmooEdges - tangled _LayerEdge's */ //================================================================================ @@ -9027,7 +9537,12 @@ void _LayerEdge::Block( _SolidData& data ) //if ( Is( BLOCKED )) return; Set( BLOCKED ); - _maxLen = _len; + SMESH_Comment msg( "#BLOCK shape="); + msg << data.GetShapeEdges( this )->_shapeID + << ", nodes " << _nodes[0]->GetID() << ", " << _nodes.back()->GetID(); + dumpCmd( msg + " -- BEGIN"); + + SetMaxLen( _len ); std::queue<_LayerEdge*> queue; queue.push( this ); @@ -9049,18 +9564,21 @@ void _LayerEdge::Block( _SolidData& data ) minDist = Min( pSrc.SquareDistance( pTgtN ), minDist ); minDist = Min( pTgt.SquareDistance( pSrcN ), minDist ); double newMaxLen = edge->_maxLen + 0.5 * Sqrt( minDist ); - if ( edge->_nodes[0]->getshapeId() == neibor->_nodes[0]->getshapeId() ) + //if ( edge->_nodes[0]->getshapeId() == neibor->_nodes[0]->getshapeId() ) viscous_layers_00/A3 { - newMaxLen *= edge->_lenFactor / neibor->_lenFactor; + //newMaxLen *= edge->_lenFactor / neibor->_lenFactor; + // newMaxLen *= Min( edge->_lenFactor / neibor->_lenFactor, + // neibor->_lenFactor / edge->_lenFactor ); } if ( neibor->_maxLen > newMaxLen ) { - neibor->_maxLen = newMaxLen; + neibor->SetMaxLen( newMaxLen ); if ( neibor->_maxLen < neibor->_len ) { _EdgesOnShape* eos = data.GetShapeEdges( neibor ); + int lastStep = neibor->Is( BLOCKED ) ? 1 : 0; while ( neibor->_len > neibor->_maxLen && - neibor->NbSteps() > 1 ) + neibor->NbSteps() > lastStep ) neibor->InvalidateStep( neibor->NbSteps(), *eos, /*restoreLength=*/true ); neibor->SetNewLength( neibor->_maxLen, *eos, data.GetHelper() ); //neibor->Block( data ); @@ -9069,6 +9587,7 @@ void _LayerEdge::Block( _SolidData& data ) } } } + dumpCmd( msg + " -- END"); } //================================================================================ @@ -9111,9 +9630,15 @@ void _LayerEdge::InvalidateStep( size_t curStep, const _EdgesOnShape& eos, bool if ( restoreLength ) { - _len -= ( nXYZ.XYZ() - curXYZ ).Modulus() / _lenFactor; + if ( NbSteps() == 0 ) + _len = 0.; + else if ( IsOnFace() && Is( MOVED )) + _len = ( nXYZ.XYZ() - SMESH_NodeXYZ( _nodes[0] )) * _normal; + else + _len -= ( nXYZ.XYZ() - curXYZ ).Modulus() / _lenFactor; } } + return; } //================================================================================ @@ -9149,10 +9674,23 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol ) int iSmoothed = GetSmoothedPos( tol ); if ( !iSmoothed ) return; - //if ( 1 || Is( DISTORTED )) + gp_XYZ normal = _normal; + if ( Is( NORMAL_UPDATED )) { - gp_XYZ normal = _normal; - if ( Is( NORMAL_UPDATED )) + double minDot = 1; + for ( size_t i = 0; i < _neibors.size(); ++i ) + { + if ( _neibors[i]->IsOnFace() ) + { + double dot = _normal * _neibors[i]->_normal; + if ( dot < minDot ) + { + normal = _neibors[i]->_normal; + minDot = dot; + } + } + } + if ( minDot == 1. ) for ( size_t i = 1; i < _pos.size(); ++i ) { normal = _pos[i] - _pos[0]; @@ -9163,42 +9701,29 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol ) break; } } - const double r = 0.2; - for ( int iter = 0; iter < 50; ++iter ) + } + const double r = 0.2; + for ( int iter = 0; iter < 50; ++iter ) + { + double minDot = 1; + for ( size_t i = Max( 1, iSmoothed-1-iter ); i < _pos.size()-1; ++i ) { - double minDot = 1; - for ( size_t i = Max( 1, iSmoothed-1-iter ); i < _pos.size()-1; ++i ) - { - gp_XYZ midPos = 0.5 * ( _pos[i-1] + _pos[i+1] ); - gp_XYZ newPos = ( 1-r ) * midPos + r * _pos[i]; - _pos[i] = newPos; - double midLen = 0.5 * ( segLen[i-1] + segLen[i+1] ); - double newLen = ( 1-r ) * midLen + r * segLen[i]; - const_cast< double& >( segLen[i] ) = newLen; - // check angle between normal and (_pos[i+1], _pos[i] ) - gp_XYZ posDir = _pos[i+1] - _pos[i]; - double size = posDir.SquareModulus(); - if ( size > RealSmall() ) - minDot = Min( minDot, ( normal * posDir ) * ( normal * posDir ) / size ); - } - if ( minDot > 0.5 * 0.5 ) - break; + gp_XYZ midPos = 0.5 * ( _pos[i-1] + _pos[i+1] ); + gp_XYZ newPos = ( 1-r ) * midPos + r * _pos[i]; + _pos[i] = newPos; + double midLen = 0.5 * ( segLen[i-1] + segLen[i+1] ); + double newLen = ( 1-r ) * midLen + r * segLen[i]; + const_cast< double& >( segLen[i] ) = newLen; + // check angle between normal and (_pos[i+1], _pos[i] ) + gp_XYZ posDir = _pos[i+1] - _pos[i]; + double size = posDir.SquareModulus(); + if ( size > RealSmall() ) + minDot = Min( minDot, ( normal * posDir ) * ( normal * posDir ) / size ); } + if ( minDot > 0.5 * 0.5 ) + break; } - // else - // { - // for ( size_t i = 1; i < _pos.size()-1; ++i ) - // { - // if ((int) i < iSmoothed && ( segLen[i] / segLen.back() < 0.5 )) - // continue; - - // double wgt = segLen[i] / segLen.back(); - // gp_XYZ normPos = _pos[0] + _normal * wgt * _len; - // gp_XYZ tgtPos = ( 1 - wgt ) * _pos[0] + wgt * _pos.back(); - // gp_XYZ newPos = ( 1 - wgt ) * normPos + wgt * tgtPos; - // _pos[i] = newPos; - // } - // } + return; } //================================================================================ @@ -9223,6 +9748,7 @@ std::string _LayerEdge::DumpFlags() const case BLOCKED: dump << "BLOCKED"; break; case INTERSECTED: dump << "INTERSECTED"; break; case NORMAL_UPDATED: dump << "NORMAL_UPDATED"; break; + case UPD_NORMAL_CONV: dump << "UPD_NORMAL_CONV"; break; case MARKED: dump << "MARKED"; break; case MULTI_NORMAL: dump << "MULTI_NORMAL"; break; case NEAR_BOUNDARY: dump << "NEAR_BOUNDARY"; break; @@ -9238,11 +9764,11 @@ std::string _LayerEdge::DumpFlags() const return dump; } + //================================================================================ /*! - case brief: - default: -*/ + * \brief Create layers of prisms + */ //================================================================================ bool _ViscousBuilder::refine(_SolidData& data) @@ -9258,7 +9784,7 @@ bool _ViscousBuilder::refine(_SolidData& data) double f,l, u = 0; gp_XY uv; vector< gp_XYZ > pos3D; - bool isOnEdge; + bool isOnEdge, isTooConvexFace = false; TGeomID prevBaseId = -1; TNode2Edge* n2eMap = 0; TNode2Edge::iterator n2e; @@ -9297,11 +9823,11 @@ bool _ViscousBuilder::refine(_SolidData& data) surface = helper.GetSurface( geomFace ); // propagate _toSmooth back to _eosC1, which was unset in findShapesToSmooth() for ( size_t i = 0; i < eos._eosC1.size(); ++i ) - { eos._eosC1[ i ]->_toSmooth = true; - for ( size_t j = 0; j < eos._eosC1[i]->_edges.size(); ++j ) - eos._eosC1[i]->_edges[j]->Set( _LayerEdge::SMOOTHED_C1 ); - } + + isTooConvexFace = false; + if ( _ConvexFace* cf = data.GetConvexFace( eos._shapeID )) + isTooConvexFace = cf->_isTooCurved; } vector< double > segLen; @@ -9317,8 +9843,8 @@ bool _ViscousBuilder::refine(_SolidData& data) if ( eos._sWOL.IsNull() ) { bool useNormal = true; - bool usePos = false; - bool smoothed = false; + bool usePos = false; + bool smoothed = false; double preci = 0.1 * edge._len; if ( eos._toSmooth && edge._pos.size() > 2 ) { @@ -9326,8 +9852,7 @@ bool _ViscousBuilder::refine(_SolidData& data) } if ( smoothed ) { - if ( !surface.IsNull() && - !data._convexFaces.count( eos._shapeID )) // edge smoothed on FACE + if ( !surface.IsNull() && !isTooConvexFace ) // edge smoothed on FACE { useNormal = usePos = false; gp_Pnt2d uv = helper.GetNodeUV( geomFace, edge._nodes[0] ); @@ -9402,8 +9927,20 @@ bool _ViscousBuilder::refine(_SolidData& data) } else if ( eos._isRegularSWOL ) // usual SWOL { - for ( size_t j = 1; j < edge._pos.size(); ++j ) - segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus(); + if ( edge.Is( _LayerEdge::SMOOTHED )) + { + SMESH_NodeXYZ p0( edge._nodes[0] ); + for ( size_t j = 1; j < edge._pos.size(); ++j ) + { + gp_XYZ pj = surface->Value( edge._pos[j].X(), edge._pos[j].Y() ).XYZ(); + segLen[j] = ( pj - p0 ) * edge._normal; + } + } + else + { + for ( size_t j = 1; j < edge._pos.size(); ++j ) + segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus(); + } } else if ( !surface.IsNull() ) // SWOL surface with singularities { @@ -11136,11 +11673,11 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) const SMDS_MeshNode* tgtN0 = ledges[0]->_nodes.back(); const SMDS_MeshNode* tgtN1 = ledges[1]->_nodes.back(); int nbSharedPyram = 0; - SMDS_ElemIteratorPtr vIt = tgtN0->GetInverseElementIterator(SMDSAbs_Volume); + SMDS_ElemIteratorPtr vIt = tgtN1->GetInverseElementIterator(SMDSAbs_Volume); while ( vIt->more() ) { const SMDS_MeshElement* v = vIt->next(); - nbSharedPyram += int( v->GetNodeIndex( tgtN1 ) >= 0 ); + nbSharedPyram += int( v->GetNodeIndex( tgtN0 ) >= 0 ); } if ( nbSharedPyram > 1 ) continue; // not free border of the pyramid