+//=======================================================================
+//function : GetHexaFacetsToSplit
+//purpose : For hexahedra that will be split into prisms, finds facets to
+// split into triangles. Only hexahedra adjacent to the one closest
+// to theFacetNormal.Location() are returned.
+//param [in,out] theHexas - the hexahedra
+//param [in] theFacetNormal - facet normal
+//param [out] theFacets - the hexahedra and found facet IDs
+//=======================================================================
+
+void SMESH_MeshEditor::GetHexaFacetsToSplit( TIDSortedElemSet& theHexas,
+ const gp_Ax1& theFacetNormal,
+ TFacetOfElem & theFacets)
+{
+ #define THIS_METHOD "SMESH_MeshEditor::GetHexaFacetsToSplit(): "
+
+ // Find a hexa closest to the location of theFacetNormal
+
+ const SMDS_MeshElement* startHex;
+ {
+ // get SMDS_ElemIteratorPtr on theHexas
+ typedef const SMDS_MeshElement* TValue;
+ typedef TIDSortedElemSet::iterator TSetIterator;
+ typedef SMDS::SimpleAccessor<TValue,TSetIterator> TAccesor;
+ typedef SMDS_MeshElement::GeomFilter TFilter;
+ typedef SMDS_SetIterator < TValue, TSetIterator, TAccesor, TFilter > TElemSetIter;
+ SMDS_ElemIteratorPtr elemIt = SMDS_ElemIteratorPtr
+ ( new TElemSetIter( theHexas.begin(),
+ theHexas.end(),
+ SMDS_MeshElement::GeomFilter( SMDSGeom_HEXA )));
+
+ SMESH_ElementSearcher* searcher =
+ SMESH_MeshAlgos::GetElementSearcher( *myMesh->GetMeshDS(), elemIt );
+
+ startHex = searcher->FindClosestTo( theFacetNormal.Location(), SMDSAbs_Volume );
+
+ delete searcher;
+
+ if ( !startHex )
+ throw SALOME_Exception( THIS_METHOD "startHex not found");
+ }
+
+ // Select a facet of startHex by theFacetNormal
+
+ SMDS_VolumeTool vTool( startHex );
+ double norm[3], dot, maxDot = 0;
+ int facetID = -1;
+ for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+ if ( vTool.GetFaceNormal( iF, norm[0], norm[1], norm[2] ))
+ {
+ dot = Abs( theFacetNormal.Direction().Dot( gp_Dir( norm[0], norm[1], norm[2] )));
+ if ( dot > maxDot )
+ {
+ facetID = iF;
+ maxDot = dot;
+ }
+ }
+ if ( facetID < 0 )
+ throw SALOME_Exception( THIS_METHOD "facet of startHex not found");
+
+ // Fill theFacets starting from facetID of startHex
+
+ // facets used for seach of volumes adjacent to already treated ones
+ typedef pair< TFacetOfElem::iterator, int > TElemFacets;
+ typedef map< TVolumeFaceKey, TElemFacets > TFacetMap;
+ TFacetMap facetsToCheck;
+
+ set<const SMDS_MeshNode*> facetNodes;
+ const SMDS_MeshElement* curHex;
+
+ const bool allHex = ( theHexas.size() == myMesh->NbHexas() );
+
+ while ( startHex )
+ {
+ // move in two directions from startHex via facetID
+ for ( int is2nd = 0; is2nd < 2; ++is2nd )
+ {
+ curHex = startHex;
+ int curFacet = facetID;
+ if ( is2nd ) // do not treat startHex twice
+ {
+ vTool.Set( curHex );
+ if ( vTool.IsFreeFace( curFacet, &curHex ))
+ {
+ curHex = 0;
+ }
+ else
+ {
+ vTool.GetFaceNodes( curFacet, facetNodes );
+ vTool.Set( curHex );
+ curFacet = vTool.GetFaceIndex( facetNodes );
+ }
+ }
+ while ( curHex )
+ {
+ // store a facet to split
+ if ( curHex->GetGeomType() != SMDSGeom_HEXA )
+ {
+ theFacets.insert( make_pair( curHex, -1 ));
+ break;
+ }
+ if ( !allHex && !theHexas.count( curHex ))
+ break;
+
+ pair< TFacetOfElem::iterator, bool > facetIt2isNew =
+ theFacets.insert( make_pair( curHex, curFacet ));
+ if ( !facetIt2isNew.second )
+ break;
+
+ // remember not-to-split facets in facetsToCheck
+ int oppFacet = vTool.GetOppFaceIndexOfHex( curFacet );
+ for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+ {
+ if ( iF == curFacet && iF == oppFacet )
+ continue;
+ TVolumeFaceKey facetKey ( vTool, iF );
+ TElemFacets elemFacet( facetIt2isNew.first, iF );
+ pair< TFacetMap::iterator, bool > it2isnew =
+ facetsToCheck.insert( make_pair( facetKey, elemFacet ));
+ if ( !it2isnew.second )
+ facetsToCheck.erase( it2isnew.first ); // adjacent hex already checked
+ }
+ // pass to a volume adjacent via oppFacet
+ if ( vTool.IsFreeFace( oppFacet, &curHex ))
+ {
+ curHex = 0;
+ }
+ else
+ {
+ // get a new curFacet
+ vTool.GetFaceNodes( oppFacet, facetNodes );
+ vTool.Set( curHex );
+ curFacet = vTool.GetFaceIndex( facetNodes, /*hint=*/curFacet );
+ }
+ }
+ } // move in two directions from startHex via facetID
+
+ // Find a new startHex by facetsToCheck
+
+ startHex = 0;
+ facetID = -1;
+ TFacetMap::iterator fIt = facetsToCheck.begin();
+ while ( !startHex && fIt != facetsToCheck.end() )
+ {
+ const TElemFacets& elemFacets = fIt->second;
+ const SMDS_MeshElement* hex = elemFacets.first->first;
+ int splitFacet = elemFacets.first->second;
+ int lateralFacet = elemFacets.second;
+ facetsToCheck.erase( fIt );
+ fIt = facetsToCheck.begin();
+
+ vTool.Set( hex );
+ if ( vTool.IsFreeFace( lateralFacet, &curHex ) ||
+ curHex->GetGeomType() != SMDSGeom_HEXA )
+ continue;
+ if ( !allHex && !theHexas.count( curHex ))
+ continue;
+
+ startHex = curHex;
+
+ // find a facet of startHex to split
+
+ set<const SMDS_MeshNode*> lateralNodes;
+ vTool.GetFaceNodes( lateralFacet, lateralNodes );
+ vTool.GetFaceNodes( splitFacet, facetNodes );
+ int oppLateralFacet = vTool.GetOppFaceIndexOfHex( lateralFacet );
+ vTool.Set( startHex );
+ lateralFacet = vTool.GetFaceIndex( lateralNodes, oppLateralFacet );
+
+ // look for a facet of startHex having common nodes with facetNodes
+ // but not lateralFacet
+ for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+ {
+ if ( iF == lateralFacet )
+ continue;
+ int nbCommonNodes = 0;
+ const SMDS_MeshNode** nn = vTool.GetFaceNodes( iF );
+ for ( int iN = 0, nbN = vTool.NbFaceNodes( iF ); iN < nbN; ++iN )
+ nbCommonNodes += facetNodes.count( nn[ iN ]);
+
+ if ( nbCommonNodes >= 2 )
+ {
+ facetID = iF;
+ break;
+ }
+ }
+ if ( facetID < 0 )
+ throw SALOME_Exception( THIS_METHOD "facet of a new startHex not found");
+ }
+ } // while ( startHex )
+}
+