X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=76f96eae201a68043f9eea3c03060b53273102af;hp=f84fc5b5e791102a8d181afa598a336bcce75639;hb=6466454bba7e473e87ce90458fb2b458dbbc2518;hpb=7211ec8d0ebb65251f41eeb32b74b1e6e322a875 diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index f84fc5b5e..76f96eae2 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -409,7 +409,7 @@ namespace VISCOUS_3D gp_XYZ _normal; // to boundary of solid vector _pos; // points computed during inflation - double _len; // length achived with the last inflation step + double _len; // length achieved with the last inflation step double _maxLen; // maximal possible length double _cosin; // of angle (_normal ^ surface) double _minAngle; // of _simplices @@ -427,24 +427,26 @@ 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 + 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; } void Unset( int flag ) { _flags &= ~flag; } + std::string DumpFlags() const; // debug void SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelper& helper ); bool SetNewLength2d( Handle(Geom_Surface)& surface, @@ -496,8 +498,9 @@ 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 ); } gp_XYZ Copy( _LayerEdge& other, _EdgesOnShape& eos, SMESH_MesherHelper& helper ); void SetCosin( double cosin ); void SetNormal( const gp_XYZ& n ) { _normal = n; } @@ -663,6 +666,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 @@ -677,7 +682,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 */ @@ -691,7 +696,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, @@ -756,8 +768,6 @@ namespace VISCOUS_3D int _nbShapesToSmooth; - //map< TGeomID,Handle(Geom_Curve)> _edge2curve; - vector< _CollisionEdges > _collisionEdges; set< TGeomID > _concaveFaces; @@ -869,6 +879,7 @@ namespace VISCOUS_3D const gp_XY& uvToFix, const double refSign ); }; + struct PyDump; //-------------------------------------------------------------------------------- /*! * \brief Builder of viscous layers @@ -941,8 +952,11 @@ 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, @@ -986,6 +1000,7 @@ namespace VISCOUS_3D TopTools_MapOfShape _shrinkedFaces; int _tmpFaceID; + PyDump* _pyDump; }; //-------------------------------------------------------------------------------- /*! @@ -1033,6 +1048,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, @@ -1046,31 +1062,24 @@ 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 ]; @@ -1662,14 +1671,15 @@ 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__);} @@ -1709,7 +1727,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) @@ -1851,6 +1869,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()) @@ -1864,15 +1883,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; @@ -2772,8 +2795,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(); @@ -2785,58 +2806,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 ); @@ -2849,6 +2839,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 ); @@ -2860,14 +2851,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->_maxLen = Min( ledge->_maxLen, 1. / curvature ); + } } } continue; @@ -2937,18 +2927,9 @@ 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 - // boundry inclined to the shape at a sharp angle + // boundary inclined to the shape at a sharp angle //list< TGeomID > shapesToSmooth; TopTools_MapOfShape edgesOfSmooFaces; @@ -3044,8 +3025,8 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) { 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 ); } } } @@ -3424,11 +3405,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 { @@ -3530,14 +3512,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 ); @@ -4157,7 +4144,7 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, //================================================================================ /*! * \brief Copy data from a _LayerEdge of other SOLID and based on the same node; - * this and other _LayerEdge's are inflated along a FACE or an EDGE + * this and the other _LayerEdge are inflated along a FACE or an EDGE */ //================================================================================ @@ -4368,7 +4355,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() ) @@ -4398,6 +4385,78 @@ void _ViscousBuilder::computeGeomSize( _SolidData& data ) } } } + + 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; } //================================================================================ @@ -4410,24 +4469,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 && @@ -4507,7 +4562,10 @@ bool _ViscousBuilder::inflate(_SolidData& data) const double shapeTgtThick = eos._hyp.GetTotalThickness(); for ( size_t i = 0; i < eos._edges.size(); ++i ) { - avgThick += Min( 1., eos._edges[i]->_len / shapeTgtThick ); + if ( eos._edges[i]->_nodes.size() > 1 ) + avgThick += Min( 1., eos._edges[i]->_len / shapeTgtThick ); + else + avgThick += shapeTgtThick; nbActiveEdges += ( ! eos._edges[i]->Is( _LayerEdge::BLOCKED )); } } @@ -4529,6 +4587,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) break; } #endif + // new step size limitStepSize( data, 0.25 * distToIntersection ); if ( data._stepSizeNodes[0] ) @@ -4850,6 +4909,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 ) @@ -4951,7 +5011,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 @@ -4963,15 +5023,16 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, if ( dist > 0 ) { 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 { @@ -5256,7 +5317,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 @@ -5287,8 +5348,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 ) @@ -5298,7 +5364,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 @@ -5317,9 +5383,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; @@ -5328,7 +5400,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 ) { @@ -5338,6 +5410,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; + } } //================================================================================ @@ -5435,6 +5527,106 @@ 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, _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, _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 @@ -5448,80 +5640,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 )) - { - 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 )) + //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 ) { - 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; @@ -5560,9 +5793,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(); @@ -5594,9 +5828,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 ); @@ -5609,6 +5844,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; @@ -5631,7 +5868,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 )) @@ -5675,10 +5914,13 @@ 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]; for ( int is2nd = 0; is2nd < 2; ++is2nd ) @@ -5734,12 +5976,18 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, 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; + 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 + // --------------------------------------------------------------------------------- size_t iSeg = 0, nbSeg = _iSeg[1] - _iSeg[0] + 1; vector< double > len( nbSeg + 1 ); @@ -5763,12 +6011,15 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, _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 ]); @@ -5777,17 +6028,17 @@ 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 ); } @@ -5828,15 +6079,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 ) { @@ -5846,17 +6109,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 ); @@ -5895,11 +6181,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; @@ -5933,7 +6222,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 */ //================================================================================ @@ -5944,6 +6233,9 @@ 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; } @@ -6482,6 +6774,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 @@ -6589,8 +6955,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; @@ -6631,9 +6997,14 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, e2neIt->second._maxLen = 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._maxLen = 1.3 * minIntDist / edge1->_lenFactor; + } if ( iter > 0 && sgn1 * sgn2 < 0 && edge2->_cosin < 0 ) e2neIt->second._normal += dir1; } @@ -6650,9 +7021,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(); @@ -6934,6 +7306,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() ) @@ -7296,6 +7670,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 @@ -7504,7 +7931,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 ]; @@ -7578,13 +8005,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() ); @@ -7593,7 +8021,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(); @@ -7671,7 +8099,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 */ //================================================================================ @@ -8507,7 +8935,7 @@ gp_XYZ _LayerEdge::smoothNefPolygon() { gp_XYZ newPos(0,0,0); - // get a plane to seach a solution on + // get a plane to search a solution on vector< gp_XYZ > vecs( _simplices.size() + 1 ); size_t i; @@ -8698,7 +9126,7 @@ gp_XYZ _LayerEdge::smoothNefPolygon() { ////////////////////////////////// NEW gp_XYZ newPos(0,0,0); - // get a plane to seach a solution on + // get a plane to search a solution on size_t i; gp_XYZ center(0,0,0); @@ -9023,6 +9451,11 @@ void _LayerEdge::Block( _SolidData& data ) //if ( Is( BLOCKED )) return; Set( BLOCKED ); + SMESH_Comment msg( "#BLOCK shape="); + msg << data.GetShapeEdges( this )->_shapeID + << ", nodes " << _nodes[0]->GetID() << ", " << _nodes.back()->GetID(); + dumpCmd( msg + " -- BEGIN") + _maxLen = _len; std::queue<_LayerEdge*> queue; queue.push( this ); @@ -9045,9 +9478,11 @@ 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 *= Min( edge->_lenFactor / neibor->_lenFactor, + // neibor->_lenFactor / edge->_lenFactor ); } if ( neibor->_maxLen > newMaxLen ) { @@ -9056,7 +9491,7 @@ void _LayerEdge::Block( _SolidData& data ) { _EdgesOnShape* eos = data.GetShapeEdges( neibor ); while ( neibor->_len > neibor->_maxLen && - neibor->NbSteps() > 1 ) + neibor->NbSteps() > 0 ) neibor->InvalidateStep( neibor->NbSteps(), *eos, /*restoreLength=*/true ); neibor->SetNewLength( neibor->_maxLen, *eos, data.GetHelper() ); //neibor->Block( data ); @@ -9065,6 +9500,7 @@ void _LayerEdge::Block( _SolidData& data ) } } } + dumpCmd( msg + " -- END") } //================================================================================ @@ -9107,7 +9543,10 @@ void _LayerEdge::InvalidateStep( size_t curStep, const _EdgesOnShape& eos, bool if ( restoreLength ) { - _len -= ( nXYZ.XYZ() - curXYZ ).Modulus() / _lenFactor; + if ( NbSteps() == 0 ) + _len = 0.; + else + _len -= ( nXYZ.XYZ() - curXYZ ).Modulus() / _lenFactor; } } } @@ -9199,10 +9638,49 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol ) //================================================================================ /*! - * \brief Create layers of prisms + * \brief Print flags */ //================================================================================ +std::string _LayerEdge::DumpFlags() const +{ + SMESH_Comment dump; + for ( int flag = 1; flag < 0x1000000; flag *= 2 ) + if ( _flags & flag ) + { + EFlags f = (EFlags) flag; + switch ( f ) { + case TO_SMOOTH: dump << "TO_SMOOTH"; break; + case MOVED: dump << "MOVED"; break; + case SMOOTHED: dump << "SMOOTHED"; break; + case DIFFICULT: dump << "DIFFICULT"; break; + case ON_CONCAVE_FACE: dump << "ON_CONCAVE_FACE"; break; + 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; + case SMOOTHED_C1: dump << "SMOOTHED_C1"; break; + case DISTORTED: dump << "DISTORTED"; break; + case RISKY_SWOL: dump << "RISKY_SWOL"; break; + case SHRUNK: dump << "SHRUNK"; break; + case UNUSED_FLAG: dump << "UNUSED_FLAG"; break; + } + dump << " "; + } + cout << dump << endl; + return dump; +} + +//================================================================================ +/*! + case brief: + default: +*/ +//================================================================================ + bool _ViscousBuilder::refine(_SolidData& data) { SMESH_MesherHelper& helper = data.GetHelper(); @@ -9216,7 +9694,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; @@ -9260,6 +9738,9 @@ bool _ViscousBuilder::refine(_SolidData& data) 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; @@ -9275,8 +9756,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 ) { @@ -9284,8 +9765,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] ); @@ -9360,8 +9840,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 { @@ -9765,6 +10257,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData) continue; _SolidData& data = *dataList.front(); + _SolidData* data2 = dataList.size() > 1 ? dataList.back() : 0; const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first )); SMESH_subMesh* sm = _mesh->GetSubMesh( F ); SMESHDS_SubMesh* smDS = sm->GetSubMeshDS(); @@ -9778,7 +10271,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData) // Prepare data for shrinking // =========================== - // Collect nodes to smooth, they are marked at the beginning of this method + // Collect nodes to smooth (they are marked at the beginning of this method) vector < const SMDS_MeshNode* > smoothNodes; { SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); @@ -9816,8 +10309,13 @@ bool _ViscousBuilder::shrink(_SolidData& theData) if ( data._noShrinkShapes.count( subID )) continue; _EdgesOnShape* eos = data.GetShapeEdges( subID ); - if ( !eos || eos->_sWOL.IsNull() ) continue; - + if ( !eos || eos->_sWOL.IsNull() ) + if ( data2 ) // check in adjacent SOLID + { + eos = data2->GetShapeEdges( subID ); + if ( !eos || eos->_sWOL.IsNull() ) + continue; + } subEOS.push_back( eos ); for ( size_t i = 0; i < eos->_edges.size(); ++i ) @@ -9877,7 +10375,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData) { const SMDS_MeshNode* n = smoothNodes[i]; nodesToSmooth[ i ]._node = n; - // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices + // src nodes must be already replaced by tgt nodes to have tgt nodes in _simplices _Simplex::GetSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, 0, sortSimplices); // fix up incorrect uv of nodes on the FACE helper.GetNodeUV( F, n, 0, &isOkUV); @@ -10246,6 +10744,8 @@ bool _ViscousBuilder::shrink(_SolidData& theData) // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh VISCOUS_3D::ToClearSubWithMain( sm, data._solid ); + if ( data2 ) + VISCOUS_3D::ToClearSubWithMain( sm, data2->_solid ); } // loop on FACES to srink mesh on @@ -10324,7 +10824,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, double uSrc = helper.GetNodeU( E, srcNode, n2 ); double uTgt = helper.GetNodeU( E, tgtNode, srcNode ); - double u2 = helper.GetNodeU( E, n2, srcNode ); + double u2 = helper.GetNodeU( E, n2, srcNode ); //edge._pos.clear(); @@ -10859,7 +11359,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, _done = false; } // check _LayerEdge - if ( e == _edges[0] || e == _edges[1] ) + if ( e == _edges[0] || e == _edges[1] || e->_nodes.size() < 2 ) return; if ( eos.SWOLType() != TopAbs_EDGE ) throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added")); @@ -11060,7 +11560,8 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) for ( int iE = 1; iE <= geomEdges.Extent(); ++iE ) { const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE)); - if ( data._noShrinkShapes.count( getMeshDS()->ShapeToIndex( E ))) + const TGeomID edgeID = getMeshDS()->ShapeToIndex( E ); + if ( data._noShrinkShapes.count( edgeID )) continue; // Get _LayerEdge's based on E @@ -11109,10 +11610,9 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) // Find out orientation and type of face to create bool reverse = false, isOnFace; - - map< TGeomID, TopoDS_Shape >::iterator e2f = - data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E )); TopoDS_Shape F; + + map< TGeomID, TopoDS_Shape >::iterator e2f = data._shrinkShape2Shape.find( edgeID ); if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() ))) { F = e2f->second.Oriented( TopAbs_FORWARD ); @@ -11138,18 +11638,44 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) if ( !sm ) return error("error in addBoundaryElements()", data._index); + // Find a proxy sub-mesh of the FACE of an adjacent SOLID, which will use the new boundary + // faces for 3D meshing (PAL23414) + SMESHDS_SubMesh* adjSM = 0; + if ( isOnFace ) + { + const TGeomID faceID = sm->GetID(); + PShapeIteratorPtr soIt = helper.GetAncestors( F, *_mesh, TopAbs_SOLID ); + while ( const TopoDS_Shape* solid = soIt->next() ) + if ( !solid->IsSame( data._solid )) + { + size_t iData = _solids.FindIndex( *solid ) - 1; + if ( iData < _sdVec.size() && + _sdVec[ iData ]._ignoreFaceIds.count( faceID ) && + _sdVec[ iData ]._shrinkShape2Shape.count( edgeID ) == 0 ) + { + SMESH_ProxyMesh::SubMesh* proxySub = + _sdVec[ iData ]._proxyMesh->getFaceSubM( TopoDS::Face( F ), /*create=*/false); + if ( proxySub && proxySub->NbElements() > 0 ) + adjSM = proxySub; + } + } + } + // Make faces const int dj1 = reverse ? 0 : 1; const int dj2 = reverse ? 1 : 0; + vector< const SMDS_MeshElement*> ff; // new faces row + SMESHDS_Mesh* m = getMeshDS(); for ( size_t j = 1; j < ledges.size(); ++j ) { vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes; vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes; + ff.resize( std::max( nn1.size(), nn2.size() ), NULL ); if ( nn1.size() == nn2.size() ) { if ( isOnFace ) for ( size_t z = 1; z < nn1.size(); ++z ) - sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] )); + sm->AddElement( ff[z-1] = m->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] )); else for ( size_t z = 1; z < nn1.size(); ++z ) sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z] )); @@ -11158,7 +11684,7 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) { if ( isOnFace ) for ( size_t z = 1; z < nn2.size(); ++z ) - sm->AddElement( getMeshDS()->AddFace( nn1[0], nn2[z-1], nn2[z] )); + sm->AddElement( ff[z-1] = m->AddFace( nn1[0], nn2[z-1], nn2[z] )); else for ( size_t z = 1; z < nn2.size(); ++z ) sm->AddElement( new SMDS_FaceOfNodes( nn1[0], nn2[z-1], nn2[z] )); @@ -11167,11 +11693,19 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) { if ( isOnFace ) for ( size_t z = 1; z < nn1.size(); ++z ) - sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[0], nn1[z] )); + sm->AddElement( ff[z-1] = m->AddFace( nn1[z-1], nn2[0], nn1[z] )); else for ( size_t z = 1; z < nn1.size(); ++z ) sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[0], nn2[z] )); } + + if ( adjSM ) // add faces to a proxy SM of the adjacent SOLID + { + for ( size_t z = 0; z < ff.size(); ++z ) + if ( ff[ z ]) + adjSM->AddElement( ff[ z ]); + ff.clear(); + } } // Make edges @@ -11182,7 +11716,7 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) if ( eos && eos->SWOLType() == TopAbs_EDGE ) { vector< const SMDS_MeshNode*>& nn = edge->_nodes; - if ( nn.size() < 2 || nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() ) + if ( nn.size() < 2 || nn[1]->NbInverseElements( SMDSAbs_Edge ) >= 2 ) continue; helper.SetSubShape( eos->_sWOL ); helper.SetElementsOnShape( true );