+
+ // add not split hexadrons to the mesh
+ int nbAdded = 0;
+ vector<int> intHexInd( nbIntHex );
+ nbIntHex = 0;
+ for ( size_t i = 0; i < intersectedHex.size(); ++i )
+ {
+ Hexahedron * & hex = intersectedHex[ i ];
+ if ( hex )
+ {
+ intHexInd[ nbIntHex++ ] = i;
+ if ( hex->_nbIntNodes > 0 ) continue;
+ init( hex->_i, hex->_j, hex->_k );
+ }
+ else
+ {
+ init( i );
+ }
+ if ( _nbCornerNodes == 8 && ( _nbBndNodes < _nbCornerNodes || !isInHole() ))
+ {
+ // order of _hexNodes is defined by enum SMESH_Block::TShapeID
+ SMDS_MeshElement* el =
+ mesh->AddVolume( _hexNodes[0].Node(), _hexNodes[2].Node(),
+ _hexNodes[3].Node(), _hexNodes[1].Node(),
+ _hexNodes[4].Node(), _hexNodes[6].Node(),
+ _hexNodes[7].Node(), _hexNodes[5].Node() );
+ mesh->SetMeshElementOnShape( el, helper.GetSubShapeID() );
+ ++nbAdded;
+ if ( hex )
+ {
+ delete hex;
+ intersectedHex[ i ] = 0;
+ --nbIntHex;
+ }
+ }
+ else if ( _nbCornerNodes > 3 && !hex )
+ {
+ // all intersection of hex with geometry are at grid nodes
+ hex = new Hexahedron( *this );
+ hex->init( i );
+ intHexInd.push_back(0);
+ intHexInd[ nbIntHex++ ] = i;
+ }
+ }
+
+ // add elements resulted from hexadron intersection
+#ifdef WITH_TBB
+ intHexInd.resize( nbIntHex );
+ tbb::parallel_for ( tbb::blocked_range<size_t>( 0, nbIntHex ),
+ ParallelHexahedron( intersectedHex, intHexInd ),
+ tbb::simple_partitioner()); // ComputeElements() is called here
+ for ( size_t i = 0; i < intHexInd.size(); ++i )
+ if ( Hexahedron * hex = intersectedHex[ intHexInd[ i ]] )
+ nbAdded += hex->addElements( helper );
+#else
+ for ( size_t i = 0; i < intHexInd.size(); ++i )
+ if ( Hexahedron * hex = intersectedHex[ intHexInd[ i ]] )
+ {
+ hex->ComputeElements();
+ nbAdded += hex->addElements( helper );
+ }
+#endif
+
+ for ( size_t i = 0; i < intersectedHex.size(); ++i )
+ if ( intersectedHex[ i ] )
+ delete intersectedHex[ i ];
+
+ return nbAdded;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Adds computed elements to the mesh
+ */
+ int Hexahedron::addElements(SMESH_MesherHelper& helper)
+ {
+ int nbAdded = 0;
+ // add elements resulted from hexahedron intersection
+ //for ( size_t i = 0; i < _volumeDefs.size(); ++i )
+ {
+ vector< const SMDS_MeshNode* >& nodes = _volumeDefs._nodes;
+
+ if ( !_volumeDefs._quantities.empty() )
+ {
+ helper.AddPolyhedralVolume( nodes, _volumeDefs._quantities );
+ }
+ else
+ {
+ switch ( nodes.size() )
+ {
+ case 8: helper.AddVolume( nodes[0],nodes[1],nodes[2],nodes[3],
+ nodes[4],nodes[5],nodes[6],nodes[7] );
+ break;
+ case 4: helper.AddVolume( nodes[0],nodes[1],nodes[2],nodes[3] );
+ break;
+ case 6: helper.AddVolume( nodes[0],nodes[1],nodes[2],nodes[3], nodes[4],nodes[5] );
+ break;
+ case 5:
+ helper.AddVolume( nodes[0],nodes[1],nodes[2],nodes[3],nodes[4] );
+ break;
+ }
+ }
+ nbAdded += int ( _volumeDefs._nodes.size() > 0 );
+ }
+
+ return nbAdded;
+ }
+ //================================================================================
+ /*!
+ * \brief Return true if the element is in a hole
+ */
+ bool Hexahedron::isInHole() const
+ {
+ const int ijk[3] = { _i, _j, _k };
+ IntersectionPoint curIntPnt;
+
+ // consider a cell to be in a hole if all links in any direction
+ // comes OUT of geometry
+ for ( int iDir = 0; iDir < 3; ++iDir )
+ {
+ const vector<double>& coords = _grid->_coords[ iDir ];
+ LineIndexer li = _grid->GetLineIndexer( iDir );
+ li.SetIJK( _i,_j,_k );
+ size_t lineIndex[4] = { li.LineIndex (),
+ li.LineIndex10(),
+ li.LineIndex01(),
+ li.LineIndex11() };
+ bool allLinksOut = true, hasLinks = false;
+ for ( int iL = 0; iL < 4 && allLinksOut; ++iL ) // loop on 4 links parallel to iDir
+ {
+ const _Link& link = _hexLinks[ iL + 4*iDir ];
+ // 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];
+ const GridLine& line = _grid->_lines[ iDir ][ lineIndex[ iL ]];
+ multiset< IntersectionPoint >::const_iterator ip =
+ line._intPoints.upper_bound( curIntPnt );
+ --ip;
+ firstIntPnt = &(*ip);
+ }
+ else if ( !link._intNodes.empty() )
+ {
+ firstIntPnt = link._intNodes[0]._intPoint;
+ }
+
+ if ( firstIntPnt )
+ {
+ hasLinks = true;
+ allLinksOut = ( firstIntPnt->_transition == Trans_OUT );
+ }
+ }
+ if ( hasLinks && allLinksOut )
+ return true;
+ }
+ return false;