From 5c991bb2b770a109e11bd9aa889ff79301446c64 Mon Sep 17 00:00:00 2001 From: eap Date: Tue, 17 Jan 2012 13:17:19 +0000 Subject: [PATCH] 0021468: EDF 2073 SMESH: Body-fitting algo creates elements in hole --- src/StdMeshers/StdMeshers_Cartesian_3D.cxx | 84 +++++++++++++++++----- 1 file changed, 67 insertions(+), 17 deletions(-) diff --git a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx index 800a9735d..dcce058db 100644 --- a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx +++ b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx @@ -173,9 +173,9 @@ namespace _name1 = nv1; _name2 = nv2; _nameConst = nConst; } - size_t I() { return _curInd[0]; } - size_t J() { return _curInd[1]; } - size_t K() { return _curInd[2]; } + size_t I() const { return _curInd[0]; } + size_t J() const { return _curInd[1]; } + size_t K() const { return _curInd[2]; } void SetIJK( size_t i, size_t j, size_t k ) { _curInd[0] = i; _curInd[1] = j; _curInd[2] = k; @@ -204,6 +204,7 @@ namespace double _tol, _minCellSize; vector< const SMDS_MeshNode* > _nodes; // mesh nodes at grid nodes + vector< bool > _isBndNode; // is mesh node at intersection with geometry size_t NodeIndex( size_t i, size_t j, size_t k ) const { @@ -321,6 +322,7 @@ namespace vector< _Node> _intNodes; // _Node's at GridLine intersections vector< _Link > _splits; vector< _Face*> _faces; + const GridLine* _line; }; // -------------------------------------------------------------------------------- struct _OrientedLink @@ -336,11 +338,6 @@ namespace } _Node* FirstNode() const { return _link->_nodes[ _reverse ]; } _Node* LastNode() const { return _link->_nodes[ !_reverse ]; } - // int NbNodes() const { return 2 + _link->_intNodes.size(); } - // _Node* GetNode(const int i) - // { - // return ( 0 < i && i < NbNodes()-1 ) ? _link->_intNodes[i-1] : ( _link->_nodes[bool(i)]); - // } }; // -------------------------------------------------------------------------------- struct _Face @@ -361,13 +358,14 @@ namespace double _sizeThreshold, _sideLength[3]; - int _nbCornerNodes, _nbIntNodes; + int _nbCornerNodes, _nbIntNodes, _nbBndNodes; public: Hexahedron(const double sizeThreshold, Grid* grid); void Init( size_t i, size_t j, size_t k ); int MakeElements(SMESH_MesherHelper& helper); private: + bool isInHole() const; bool checkPolyhedronSize() const; bool addHexa (SMESH_MesherHelper& helper); bool addTetra(SMESH_MesherHelper& helper); @@ -528,8 +526,10 @@ namespace void Grid::ComputeNodes(SMESH_MesherHelper& helper) { // state of each node of the grid relative to the geomerty - vector< bool > isNodeOut( _coords[0].size() * _coords[1].size() * _coords[2].size(), false ); - _nodes.resize( isNodeOut.size(), 0 ); + const size_t nbGridNodes = _coords[0].size() * _coords[1].size() * _coords[2].size(); + vector< bool > isNodeOut( nbGridNodes, false ); + _nodes.resize( nbGridNodes, 0 ); + _isBndNode.resize( nbGridNodes, false ); for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions { @@ -590,7 +590,9 @@ namespace li.SetIndexOnLine( nodeCoord-coord0 ); double xyz[3] = { _coords[0][ li.I() ], _coords[1][ li.J() ], _coords[2][ li.K() ]}; _nodes[ nodeIndex ] = helper.AddNode( xyz[0], xyz[1], xyz[2] ); + _isBndNode[ nodeIndex ] = true; } + ip->_node = _nodes[ nodeIndex ]; if ( ++nodeCoord < coordEnd ) nodeParam = *nodeCoord - *coord0; } @@ -1059,12 +1061,14 @@ namespace void Hexahedron::Init( size_t i, size_t j, size_t k ) { // set nodes of grid to nodes of the hexahedron and - // count nodes at hexahedron corners located IN geometry - _nbCornerNodes = _nbIntNodes = 0; + // count nodes at hexahedron corners located IN and ON geometry + _nbCornerNodes = _nbIntNodes = _nbBndNodes = 0; size_t i000 = _grid->NodeIndex( i,j,k ); for ( int iN = 0; iN < 8; ++iN ) + { _nbCornerNodes += bool(( _hexNodes[iN]._node = _grid->_nodes[ i000 + _nodeShift[iN] ])); - + _nbBndNodes += _grid->_isBndNode[ i000 + _nodeShift[iN] ]; + } // set intersection nodes from GridLine's to hexahedron links int linkID = 0; _Link split; @@ -1087,6 +1091,7 @@ namespace { GridLine& line = _grid->_lines[ iDir ][ lineIndex[ iL ]]; _Link& link = _hexLinks[ linkID++ ]; + link._line = &line; link._intNodes.clear(); link._splits.clear(); split._nodes[ 0 ] = link._nodes[0]; @@ -1123,7 +1128,7 @@ namespace int Hexahedron::MakeElements(SMESH_MesherHelper& helper) { int nbAddedVols = 0; - if ( _nbCornerNodes == 8 && _nbIntNodes == 0 ) + if ( _nbCornerNodes == 8 && _nbIntNodes == 0 && _nbBndNodes < _nbCornerNodes ) { // order of _hexNodes is defined by enum SMESH_Block::TShapeID helper.AddVolume( _hexNodes[0].Node(), _hexNodes[2].Node(), @@ -1135,6 +1140,9 @@ namespace if ( _nbCornerNodes + _nbIntNodes < 4 ) return nbAddedVols; + if ( _nbBndNodes == _nbCornerNodes && isInHole() ) + return nbAddedVols; + _polygons.clear(); vector polyhedraNodes; @@ -1148,7 +1156,7 @@ namespace _Link polyLink; polyLink._faces.reserve( 1 ); - for ( int iF = 0; iF < 6; ++iF ) + for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides on a hexahedron { const _Face& quad = _hexQuads[ iF ] ; @@ -1160,7 +1168,7 @@ namespace // add splits of a link to a polygon and collect info on nodes //int nbIn = 0, nbOut = 0, nbCorners = 0; nodes.clear(); - for ( int iE = 0; iE < 4; ++iE ) + for ( int iE = 0; iE < 4; ++iE ) // loop on 4 sides of a quadrangle { int nbSpits = quad._links[ iE ].NbResultLinks(); for ( int iS = 0; iS < nbSpits; ++iS ) @@ -1306,6 +1314,48 @@ namespace return nbAddedVols; } + //================================================================================ + /*! + * \brief Return true if the element is in a hole + */ + bool Hexahedron::isInHole() const + { + const int ijk[3] = { _lineInd[0].I(), _lineInd[0].J(), _lineInd[0].K() }; + IntersectionPoint curIntPnt; + + for ( int iDir = 0; iDir < 3; ++iDir ) + { + const vector& coords = _grid->_coords[ iDir ]; + bool allLinksOut = true; + int linkID = iDir * 4; + for ( int i = 0; i < 4 && allLinksOut; ++i ) + { + const _Link& link = _hexLinks[ linkID++ ]; + if ( link._splits.empty() ) continue; + // check transition of the first node of a link + const IntersectionPoint* firstIntPnt = 0; + if ( link._nodes[0]->Node() ) // 1st node is a hexa corner + { + curIntPnt._paramOnLine = coords[ ijk[ iDir ]] - coords[0]; + multiset< IntersectionPoint >::const_iterator ip = + link._line->_intPoints.upper_bound( curIntPnt ); + --ip; + firstIntPnt = &(*ip); + } + else if ( !link._intNodes.empty() ) + { + firstIntPnt = link._intNodes[0]._intPoint; + } + + if ( firstIntPnt && firstIntPnt->_transition == Trans_IN ) + allLinksOut = false; + } + if ( allLinksOut ) + return true; + } + return false; + } + //================================================================================ /*! * \brief Return true if a polyhedron passes _sizeThreshold criterion -- 2.39.2