From ab008b333babd39d821e2dbf34de4f72f1c7e348 Mon Sep 17 00:00:00 2001 From: eap Date: Mon, 19 Jun 2017 18:42:49 +0300 Subject: [PATCH] 0023444: EDF 14244 - Viscous Layer + 54237: Invalid viscous layers (in Grids/smesh/viscous_layers_02/A3) + Improve viscous layers performance - by avoiding unnecessary intersection checks - by using std::vector instead of std::set in SMESH_ElementSearcher + Fix minor errors in ObjectPool which is now used in SMESH_ElementSearcher --- src/SMDS/ObjectPool.hxx | 116 +++++++++-- src/SMESHUtils/SMESH_MeshAlgos.cxx | 195 ++++++++++++------ src/StdMeshers/StdMeshers_ViscousLayers.cxx | 211 ++++++++++++++++---- 3 files changed, 399 insertions(+), 123 deletions(-) diff --git a/src/SMDS/ObjectPool.hxx b/src/SMDS/ObjectPool.hxx index 5e5f0d289..ef514b353 100644 --- a/src/SMDS/ObjectPool.hxx +++ b/src/SMDS/ObjectPool.hxx @@ -21,9 +21,10 @@ #define _OBJECTPOOL_HXX_ #include -//#include #include +#include "SMDS_Iterator.hxx" + namespace { // assure deallocation of memory of a vector @@ -33,18 +34,22 @@ namespace } } +template class ObjectPoolIterator; + template class ObjectPool { private: - std::vector _chunkList; + std::vector _chunkList; std::vector _freeList; - int _nextFree; - int _maxAvail; - int _chunkSize; - int _maxOccupied; - int _nbHoles; - int _lastDelChunk; + int _nextFree; // either the 1st hole or last added + int _maxAvail; // nb allocated elements + int _chunkSize; + int _maxOccupied; // max used ID + int _nbHoles; + int _lastDelChunk; + + friend class ObjectPoolIterator; int getNextFree() { @@ -76,16 +81,16 @@ private: } public: - ObjectPool(int nblk) + ObjectPool(int nblk = 1024) { - _chunkSize = nblk; - _nextFree = 0; - _maxAvail = 0; - _maxOccupied = 0; - _nbHoles = 0; + _chunkSize = nblk; + _nextFree = 0; + _maxAvail = 0; + _maxOccupied = -1; + _nbHoles = 0; + _lastDelChunk = 0; _chunkList.clear(); _freeList.clear(); - _lastDelChunk = 0; } virtual ~ObjectPool() @@ -105,16 +110,16 @@ public: _freeList.insert(_freeList.end(), _chunkSize, true); _maxAvail += _chunkSize; _freeList[_nextFree] = false; - obj = newChunk; // &newChunk[0]; + obj = newChunk; } else { int chunkId = _nextFree / _chunkSize; int rank = _nextFree - chunkId * _chunkSize; _freeList[_nextFree] = false; - obj = _chunkList[chunkId] + rank; // &_chunkList[chunkId][rank]; + obj = _chunkList[chunkId] + rank; } - if (_nextFree < _maxOccupied) + if (_nextFree <= _maxOccupied) { _nbHoles-=1; } @@ -122,7 +127,6 @@ public: { _maxOccupied = _nextFree; } - //obj->init(); return obj; } @@ -148,10 +152,10 @@ public: if (toFree < _nextFree) _nextFree = toFree; if (toFree < _maxOccupied) - _nbHoles += 1; + ++_nbHoles; + else + --_maxOccupied; _lastDelChunk = i; - //obj->clean(); - //checkDelete(i); compactage non fait } void clear() @@ -167,6 +171,37 @@ public: clearVector( _freeList ); } + // nb allocated elements + size_t size() const + { + return _freeList.size(); + } + + // nb used elements + size_t nbElements() const + { + return _maxOccupied + 1 - _nbHoles; + } + + // return an element w/o any check + const X* operator[]( size_t i ) const // i < size() + { + int chunkId = i / _chunkSize; + int rank = i - chunkId * _chunkSize; + return _chunkList[ chunkId ] + rank; + } + + // return only being used element + const X* at( size_t i ) const // i < size() + { + if ( i >= size() || _freeList[ i ] ) + return 0; + + int chunkId = i / _chunkSize; + int rank = i - chunkId * _chunkSize; + return _chunkList[ chunkId ] + rank; + } + // void destroy(int toFree) // { // // no control 0<= toFree < _freeList.size() @@ -177,4 +212,41 @@ public: }; +template class ObjectPoolIterator : public SMDS_Iterator +{ + const ObjectPool& _pool; + int _i, _nbFound; +public: + + ObjectPoolIterator( const ObjectPool& pool ) : _pool( pool ), _i( 0 ), _nbFound( 0 ) + { + if ( more() && _pool._freeList[ _i ] == true ) + { + next(); + --_nbFound; + } + } + + virtual bool more() + { + return ( _i <= _pool._maxOccupied && _nbFound < (int)_pool.nbElements() ); + } + + virtual const X* next() + { + const X* x = 0; + if ( more() ) + { + x = _pool[ _i ]; + + ++_nbFound; + + for ( ++_i; _i <= _pool._maxOccupied; ++_i ) + if ( _pool._freeList[ _i ] == false ) + break; + } + return x; + } +}; + #endif diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index 1eae5f781..03e6f7cf5 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.cxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.cxx @@ -228,15 +228,15 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(), double tolerance = NodeRadius ); - void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems ); - void getElementsNearLine ( const gp_Ax1& line, TIDSortedElemSet& foundElems); - void getElementsInSphere ( const gp_XYZ& center, - const double radius, TIDSortedElemSet& foundElems); - size_t getSize() { return std::max( _size, _elements.size() ); } - virtual ~ElementBndBoxTree(); + void prepare(); // !!!call it before calling the following methods!!! + void getElementsNearPoint( const gp_Pnt& point, vector& foundElems ); + void getElementsNearLine ( const gp_Ax1& line, vector& foundElems); + void getElementsInSphere ( const gp_XYZ& center, + const double radius, + vector& foundElems); protected: - ElementBndBoxTree():_size(0) {} + ElementBndBoxTree() {} SMESH_Octree* newChild() const { return new ElementBndBoxTree; } void buildChildrenData(); Bnd_B3d* buildRootBox(); @@ -245,11 +245,25 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() struct ElementBox : public Bnd_B3d { const SMDS_MeshElement* _element; - int _refCount; // an ElementBox can be included in several tree branches - ElementBox(const SMDS_MeshElement* elem, double tolerance); + bool _isMarked; + void init(const SMDS_MeshElement* elem, double tolerance); }; vector< ElementBox* > _elements; - size_t _size; + + typedef ObjectPool< ElementBox > TElementBoxPool; + + //!< allocator of ElementBox's and SMESH_TreeLimit + struct LimitAndPool : public SMESH_TreeLimit + { + TElementBoxPool _elBoPool; + std::vector< ElementBox* > _markedElems; + LimitAndPool():SMESH_TreeLimit( MaxLevel, /*minSize=*/0. ) {} + }; + LimitAndPool* getLimitAndPool() const + { + SMESH_TreeLimit* limitAndPool = const_cast< SMESH_TreeLimit* >( myLimit ); + return static_cast< LimitAndPool* >( limitAndPool ); + } }; //================================================================================ @@ -258,32 +272,27 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() */ //================================================================================ - ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt, double tolerance) - :SMESH_Octree( new SMESH_TreeLimit( MaxLevel, /*minSize=*/0. )) + ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, + SMDSAbs_ElementType elemType, + SMDS_ElemIteratorPtr theElemIt, + double tolerance) + :SMESH_Octree( new LimitAndPool() ) { int nbElems = mesh.GetMeshInfo().NbElements( elemType ); _elements.reserve( nbElems ); + TElementBoxPool& elBoPool = getLimitAndPool()->_elBoPool; + SMDS_ElemIteratorPtr elemIt = theElemIt ? theElemIt : mesh.elementsIterator( elemType ); while ( elemIt->more() ) - _elements.push_back( new ElementBox( elemIt->next(),tolerance )); - + { + ElementBox* eb = elBoPool.getNew(); + eb->init( elemIt->next(), tolerance ); + _elements.push_back( eb ); + } compute(); } - //================================================================================ - /*! - * \brief Destructor - */ - //================================================================================ - - ElementBndBoxTree::~ElementBndBoxTree() - { - for ( size_t i = 0; i < _elements.size(); ++i ) - if ( --_elements[i]->_refCount <= 0 ) - delete _elements[i]; - } - //================================================================================ /*! * \brief Return the maximal box @@ -311,14 +320,10 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() for (int j = 0; j < 8; j++) { if ( !_elements[i]->IsOut( *myChildren[j]->getBox() )) - { - _elements[i]->_refCount++; ((ElementBndBoxTree*)myChildren[j])->_elements.push_back( _elements[i]); - } } - _elements[i]->_refCount--; } - _size = _elements.size(); + //_size = _elements.size(); SMESHUtils::FreeVector( _elements ); // = _elements.clear() + free memory for (int j = 0; j < 8; j++) @@ -327,33 +332,61 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() if ((int) child->_elements.size() <= MaxNbElemsInLeaf ) child->myIsLeaf = true; - if ( child->_elements.capacity() - child->_elements.size() > 1000 ) + if ( child->isLeaf() && child->_elements.capacity() > child->_elements.size() ) SMESHUtils::CompactVector( child->_elements ); } } + //================================================================================ + /*! + * \brief Un-mark all elements + */ + //================================================================================ + + void ElementBndBoxTree::prepare() + { + // TElementBoxPool& elBoPool = getElementBoxPool(); + // for ( size_t i = 0; i < elBoPool.nbElements(); ++i ) + // const_cast< ElementBox* >( elBoPool[ i ])->_isMarked = false; + } + //================================================================================ /*! * \brief Return elements which can include the point */ //================================================================================ - void ElementBndBoxTree::getElementsNearPoint( const gp_Pnt& point, - TIDSortedElemSet& foundElems) + void ElementBndBoxTree::getElementsNearPoint( const gp_Pnt& point, + vector& foundElems) { if ( getBox()->IsOut( point.XYZ() )) return; if ( isLeaf() ) { + LimitAndPool* pool = getLimitAndPool(); + for ( size_t i = 0; i < _elements.size(); ++i ) - if ( !_elements[i]->IsOut( point.XYZ() )) - foundElems.insert( _elements[i]->_element ); + if ( !_elements[i]->IsOut( point.XYZ() ) && + !_elements[i]->_isMarked ) + { + foundElems.push_back( _elements[i]->_element ); + _elements[i]->_isMarked = true; + pool->_markedElems.push_back( _elements[i] ); + } } else { for (int i = 0; i < 8; i++) ((ElementBndBoxTree*) myChildren[i])->getElementsNearPoint( point, foundElems ); + + if ( level() == 0 ) + { + LimitAndPool* pool = getLimitAndPool(); + for ( size_t i = 0; i < pool->_markedElems.size(); ++i ) + pool->_markedElems[i]->_isMarked = false; + pool->_markedElems.clear(); + } } } @@ -363,22 +396,37 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() */ //================================================================================ - void ElementBndBoxTree::getElementsNearLine( const gp_Ax1& line, - TIDSortedElemSet& foundElems) + void ElementBndBoxTree::getElementsNearLine( const gp_Ax1& line, + vector& foundElems) { if ( getBox()->IsOut( line )) return; if ( isLeaf() ) { + LimitAndPool* pool = getLimitAndPool(); + for ( size_t i = 0; i < _elements.size(); ++i ) - if ( !_elements[i]->IsOut( line )) - foundElems.insert( _elements[i]->_element ); + if ( !_elements[i]->IsOut( line ) && + !_elements[i]->_isMarked ) + { + foundElems.push_back( _elements[i]->_element ); + _elements[i]->_isMarked = true; + pool->_markedElems.push_back( _elements[i] ); + } } else { for (int i = 0; i < 8; i++) ((ElementBndBoxTree*) myChildren[i])->getElementsNearLine( line, foundElems ); + + if ( level() == 0 ) + { + LimitAndPool* pool = getLimitAndPool(); + for ( size_t i = 0; i < pool->_markedElems.size(); ++i ) + pool->_markedElems[i]->_isMarked = false; + pool->_markedElems.clear(); + } } } @@ -388,23 +436,38 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() */ //================================================================================ - void ElementBndBoxTree::getElementsInSphere ( const gp_XYZ& center, - const double radius, - TIDSortedElemSet& foundElems) + void ElementBndBoxTree::getElementsInSphere ( const gp_XYZ& center, + const double radius, + vector& foundElems) { if ( getBox()->IsOut( center, radius )) return; if ( isLeaf() ) { + LimitAndPool* pool = getLimitAndPool(); + for ( size_t i = 0; i < _elements.size(); ++i ) - if ( !_elements[i]->IsOut( center, radius )) - foundElems.insert( _elements[i]->_element ); + if ( !_elements[i]->IsOut( center, radius ) && + !_elements[i]->_isMarked ) + { + foundElems.push_back( _elements[i]->_element ); + _elements[i]->_isMarked = true; + pool->_markedElems.push_back( _elements[i] ); + } } else { for (int i = 0; i < 8; i++) ((ElementBndBoxTree*) myChildren[i])->getElementsInSphere( center, radius, foundElems ); + + if ( level() == 0 ) + { + LimitAndPool* pool = getLimitAndPool(); + for ( size_t i = 0; i < pool->_markedElems.size(); ++i ) + pool->_markedElems[i]->_isMarked = false; + pool->_markedElems.clear(); + } } } @@ -414,13 +477,13 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() */ //================================================================================ - ElementBndBoxTree::ElementBox::ElementBox(const SMDS_MeshElement* elem, double tolerance) + void ElementBndBoxTree::ElementBox::init(const SMDS_MeshElement* elem, double tolerance) { _element = elem; - _refCount = 1; + _isMarked = false; SMDS_ElemIteratorPtr nIt = elem->nodesIterator(); while ( nIt->more() ) - Add( SMESH_TNodeXYZ( nIt->next() )); + Add( SMESH_NodeXYZ( nIt->next() )); Enlarge( tolerance ); } @@ -771,9 +834,13 @@ FindElementsByPoint(const gp_Pnt& point, { _ebbTree[_elementType] = new ElementBndBoxTree( *_mesh, type, _meshPartIt, tolerance ); } - TIDSortedElemSet suspectElems; + else + { + _ebbTree[ type ]->prepare(); + } + vector< const SMDS_MeshElement* > suspectElems; _ebbTree[ type ]->getElementsNearPoint( point, suspectElems ); - TIDSortedElemSet::iterator elem = suspectElems.begin(); + vector< const SMDS_MeshElement* >::iterator elem = suspectElems.begin(); for ( ; elem != suspectElems.end(); ++elem ) if ( !SMESH_MeshAlgos::IsOut( *elem, point, tolerance )) foundElements.push_back( *elem ); @@ -801,8 +868,10 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt& point, ElementBndBoxTree*& ebbTree = _ebbTree[ type ]; if ( !ebbTree ) ebbTree = new ElementBndBoxTree( *_mesh, type, _meshPartIt ); + else + ebbTree->prepare(); - TIDSortedElemSet suspectElems; + vector suspectElems; ebbTree->getElementsNearPoint( point, suspectElems ); if ( suspectElems.empty() && ebbTree->maxSize() > 0 ) @@ -816,13 +885,14 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt& point, radius = ebbTree->maxSize() / pow( 2., getTreeHeight()) / 2; while ( suspectElems.empty() ) { + ebbTree->prepare(); ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems ); radius *= 1.1; } } double minDist = std::numeric_limits::max(); multimap< double, const SMDS_MeshElement* > dist2face; - TIDSortedElemSet::iterator elem = suspectElems.begin(); + vector::iterator elem = suspectElems.begin(); for ( ; elem != suspectElems.end(); ++elem ) { double dist = SMESH_MeshAlgos::GetDistance( *elem, point ); @@ -886,6 +956,8 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point) ElementBndBoxTree*& ebbTree = _ebbTree[ SMDSAbs_Face ]; if ( !ebbTree ) ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt ); + else + ebbTree->prepare(); // Algo: analyse transition of a line starting at the point through mesh boundary; // try three lines parallel to axis of the coordinate system and perform rough @@ -901,13 +973,14 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point) gp_Ax1 lineAxis( point, axisDir[axis]); gp_Lin line ( lineAxis ); - TIDSortedElemSet suspectFaces; // faces possibly intersecting the line + vector suspectFaces; // faces possibly intersecting the line + if ( axis > 0 ) ebbTree->prepare(); ebbTree->getElementsNearLine( lineAxis, suspectFaces ); // Intersect faces with the line map< double, TInters > & u2inters = paramOnLine2TInters[ axis ]; - TIDSortedElemSet::iterator face = suspectFaces.begin(); + vector::iterator face = suspectFaces.begin(); for ( ; face != suspectFaces.end(); ++face ) { // get face plane @@ -1114,10 +1187,10 @@ void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1& ElementBndBoxTree*& ebbTree = _ebbTree[ type ]; if ( !ebbTree ) ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt ); + else + ebbTree->prepare(); - TIDSortedElemSet suspectFaces; // elements possibly intersecting the line - ebbTree->getElementsNearLine( line, suspectFaces ); - foundElems.assign( suspectFaces.begin(), suspectFaces.end()); + ebbTree->getElementsNearLine( line, foundElems ); } //======================================================================= @@ -1135,10 +1208,10 @@ void SMESH_ElementSearcherImpl::GetElementsInSphere( const gp_XYZ& ElementBndBoxTree*& ebbTree = _ebbTree[ type ]; if ( !ebbTree ) ebbTree = new ElementBndBoxTree( *_mesh, _elementType, _meshPartIt ); + else + ebbTree->prepare(); - TIDSortedElemSet suspectFaces; // elements possibly intersecting the line - ebbTree->getElementsInSphere( center, radius, suspectFaces ); - foundElems.assign( suspectFaces.begin(), suspectFaces.end() ); + ebbTree->getElementsInSphere( center, radius, foundElems ); } //======================================================================= diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 8212b0a8b..665cad6a5 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -2926,15 +2926,6 @@ 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 @@ -4363,7 +4354,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() ) @@ -4393,6 +4384,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; } //================================================================================ @@ -4405,14 +4468,7 @@ 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; @@ -4530,6 +4586,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) break; } #endif + // new step size limitStepSize( data, 0.25 * distToIntersection ); if ( data._stepSizeNodes[0] ) @@ -5654,17 +5711,16 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, 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 ]; - const gp_XY rangeUV = uv1 - uv0; + 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; - _eos[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() ); SMDS_MeshNode* tgtNode = const_cast( _eos[i]->_nodes.back() ); @@ -5674,6 +5730,28 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, 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; } } } @@ -5764,6 +5842,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; @@ -6011,8 +6091,8 @@ void _Smoother1D::prepare(_SolidData& data) // 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 ) { @@ -6022,14 +6102,39 @@ 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 ) + { + _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 { - 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; + std::vector< double > params( _eos.size() + 2 ); + + 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 @@ -6069,8 +6174,14 @@ void _Smoother1D::prepare(_SolidData& data) int iLBO = _offPoints.size() - 2; // last but one - _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; @@ -6104,7 +6215,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 */ //================================================================================ @@ -6115,6 +6226,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; } @@ -6834,8 +6948,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; @@ -6876,9 +6990,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; } @@ -9710,8 +9829,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 { -- 2.30.2