+ bool isClassicElem = false;
+ if ( nbNodes == 8 && _polygons.size() == 6 ) isClassicElem = addHexa();
+ else if ( nbNodes == 4 && _polygons.size() == 4 ) isClassicElem = addTetra();
+ else if ( nbNodes == 6 && _polygons.size() == 5 ) isClassicElem = addPenta();
+ else if ( nbNodes == 5 && _polygons.size() == 5 ) isClassicElem = addPyra ();
+ if ( !isClassicElem )
+ _volumeDefs.set( polyhedraNodes, quantities );
+ }
+ //================================================================================
+ /*!
+ * \brief Create elements in the mesh
+ */
+ int Hexahedron::MakeElements(SMESH_MesherHelper& helper)
+ {
+ SMESHDS_Mesh* mesh = helper.GetMeshDS();
+
+ size_t nbCells[3] = { _grid->_coords[0].size() - 1,
+ _grid->_coords[1].size() - 1,
+ _grid->_coords[2].size() - 1 };
+ const size_t nbGridCells = nbCells[0] *nbCells [1] * nbCells[2];
+ vector< Hexahedron* > intersectedHex( nbGridCells, 0 );
+ int nbIntHex = 0;
+
+ // set intersection nodes from GridLine's to links of intersectedHex
+ int i,j,k, iDirOther[3][2] = {{ 1,2 },{ 0,2 },{ 0,1 }};
+ for ( int iDir = 0; iDir < 3; ++iDir )
+ {
+ int dInd[4][3] = { {0,0,0}, {0,0,0}, {0,0,0}, {0,0,0} };
+ dInd[1][ iDirOther[iDir][0] ] = -1;
+ dInd[2][ iDirOther[iDir][1] ] = -1;
+ dInd[3][ iDirOther[iDir][0] ] = -1; dInd[3][ iDirOther[iDir][1] ] = -1;
+ // loop on GridLine's parallel to iDir
+ LineIndexer lineInd = _grid->GetLineIndexer( iDir );
+ for ( ; lineInd.More(); ++lineInd )
+ {
+ GridLine& line = _grid->_lines[ iDir ][ lineInd.LineIndex() ];
+ multiset< IntersectionPoint >::const_iterator ip = line._intPoints.begin();
+ for ( ; ip != line._intPoints.end(); ++ip )
+ {
+ if ( !ip->_node ) continue;
+ lineInd.SetIndexOnLine( ip->_indexOnLine );
+ for ( int iL = 0; iL < 4; ++iL ) // loop on 4 cells sharing a link
+ {
+ i = int(lineInd.I()) + dInd[iL][0];
+ j = int(lineInd.J()) + dInd[iL][1];
+ k = int(lineInd.K()) + dInd[iL][2];
+ if ( i < 0 || i >= nbCells[0] ||
+ j < 0 || j >= nbCells[1] ||
+ k < 0 || k >= nbCells[2] ) continue;
+
+ const size_t hexIndex = _grid->CellIndex( i,j,k );
+ Hexahedron *& hex = intersectedHex[ hexIndex ];
+ if ( !hex)
+ {
+ hex = new Hexahedron( *this );
+ hex->_i = i;
+ hex->_j = j;
+ hex->_k = k;
+ ++nbIntHex;
+ }
+ const int iLink = iL + iDir * 4;
+ hex->_hexLinks[iLink]._intNodes.push_back( _Node( 0, &(*ip) ));
+ hex->_nbIntNodes++;
+ }
+ }
+ }
+ }
+
+ // 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 )