-// void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
-// const int theMethodFlags)
-// {
-// // sdt-like iterator on coordinates of nodes of mesh element
-// typedef SMDS_StdIterator< TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator;
-// NXyzIterator xyzEnd;
-
-// SMESH_MesherHelper helper( *GetMesh());
-
-// TIDSortedElemSet::const_iterator elem = theElems.begin();
-// for ( ; elem != theElems.end(); ++elem )
-// {
-// SMDSAbs_EntityType geomType = (*elem)->GetEntityType();
-// if ( geomType <= SMDSEntity_Quad_Tetra )
-// continue; // tetra or face or edge
-
-// if ( (*elem)->IsQuadratic() )
-// {
-// // add quadratic links to the helper
-// SMDS_VolumeTool vol( *elem );
-// for ( int iF = 0; iF < vol.NbFaces(); ++iF )
-// {
-// const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iF );
-// for ( int iN = 0; iN < vol.NbFaceNodes( iF ); iN += 2)
-// helper.AddTLinkNode( fNodes[iF], fNodes[iF+2], fNodes[iF+1] );
-// }
-// helper.SetIsQuadratic( true );
-// }
-// else
-// {
-// helper.SetIsQuadratic( false );
-// }
-
-// vector<const SMDS_MeshElement* > tetras; // splits of a volume
-
-// if ( geomType == SMDSEntity_Polyhedra )
-// {
-// // Each face of a polyhedron is split into triangles and
-// // each of triangles and a cell barycenter form a tetrahedron.
-
-// SMDS_VolumeTool vol( *elem );
-
-// // make a node at barycenter
-// gp_XYZ gc = std::accumulate( NXyzIterator((*elem)->nodesIterator()), xyzEnd,gp_XYZ(0,0,0));
-// gc /= vol.NbNodes();
-// SMDS_MeshNode* gcNode = GetMeshDS()->AddNode( gc.X(), gc.Y(), gc.Z() );
-
-// for ( int iF = 0; iF < vol.NbFaces(); ++iF )
-// {
-// const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iF );
-// int nbFNodes = vol.NbFaceNodes( iF );
-// int nbTria = nbFNodes - 2;
-// bool extFace = vol.IsFaceExternal( iF );
-// SMDS_MeshElement* tet;
-// for ( int i = 0; i < nbTria; ++i )
-// {
-// if ( extFace )
-// tet = helper.AddVolume( fNodes[0], fNodes[i+1], fNodes[i+2], gcNode );
-// else
-// tet = helper.AddVolume( fNodes[0], fNodes[i+2], fNodes[i+1], gcNode );
-// tetras.push_back( tet );
-// }
-// }
-
-// }
-// else
-// {
-
-// TSplitMethod splitMethod = getSplitMethod( *elem, theMethodFlags );
-// if ( splitMethod._nbTetra < 1 ) continue;
-
-// vector<const SMDS_MeshNode*> volNodes( (*elem)->begin_nodes(), (*elem)->end_nodes());
-// }
-// }
-// }
+void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
+ const int theMethodFlags)
+{
+ // std-like iterator on coordinates of nodes of mesh element
+ typedef SMDS_StdIterator< TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator;
+ NXyzIterator xyzEnd;
+
+ SMDS_VolumeTool volTool;
+ SMESH_MesherHelper helper( *GetMesh());
+
+ SMESHDS_SubMesh* subMesh = GetMeshDS()->MeshElements(1);
+ SMESHDS_SubMesh* fSubMesh = subMesh;
+
+ SMESH_SequenceOfElemPtr newNodes, newElems;
+
+ TIDSortedElemSet::const_iterator elem = theElems.begin();
+ for ( ; elem != theElems.end(); ++elem )
+ {
+ SMDSAbs_EntityType geomType = (*elem)->GetEntityType();
+ if ( geomType <= SMDSEntity_Quad_Tetra )
+ continue; // tetra or face or ...
+
+ if ( !volTool.Set( *elem )) continue; // not volume? strange...
+
+ TSplitMethod splitMethod = getSplitMethod( volTool, theMethodFlags );
+ if ( splitMethod._nbTetra < 1 ) continue;
+
+ // find submesh to add new tetras in
+ if ( !subMesh || !subMesh->Contains( *elem ))
+ {
+ int shapeID = FindShape( *elem );
+ helper.SetSubShape( shapeID ); // helper will add tetras to the found submesh
+ subMesh = GetMeshDS()->MeshElements( shapeID );
+ }
+ int iQ;
+ if ( (*elem)->IsQuadratic() )
+ {
+ iQ = 2;
+ // add quadratic links to the helper
+ for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
+ {
+ const SMDS_MeshNode** fNodes = volTool.GetFaceNodes( iF );
+ for ( int iN = 0; iN < volTool.NbFaceNodes( iF ); iN += iQ )
+ helper.AddTLinkNode( fNodes[iF], fNodes[iF+2], fNodes[iF+1] );
+ }
+ helper.SetIsQuadratic( true );
+ }
+ else
+ {
+ iQ = 1;
+ helper.SetIsQuadratic( false );
+ }
+ vector<const SMDS_MeshNode*> nodes( (*elem)->begin_nodes(), (*elem)->end_nodes() );
+ if ( splitMethod._baryNode )
+ {
+ // make a node at barycenter
+ gp_XYZ gc( 0,0,0 );
+ gc = accumulate( NXyzIterator((*elem)->nodesIterator()), xyzEnd, gc ) / nodes.size();
+ SMDS_MeshNode* gcNode = helper.AddNode( gc.X(), gc.Y(), gc.Z() );
+ nodes.push_back( gcNode );
+ newNodes.Append( gcNode );
+ }
+
+ // make tetras
+ helper.SetElementsOnShape( true );
+ vector<const SMDS_MeshElement* > tetras( splitMethod._nbTetra ); // splits of a volume
+ const int* tetConn = splitMethod._connectivity;
+ for ( int i = 0; i < splitMethod._nbTetra; ++i, tetConn += 4 )
+ newElems.Append( tetras[ i ] = helper.AddVolume( nodes[ tetConn[0] ],
+ nodes[ tetConn[1] ],
+ nodes[ tetConn[2] ],
+ nodes[ tetConn[3] ]));
+
+ ReplaceElemInGroups( *elem, tetras, GetMeshDS() );
+
+ // Split faces on sides of the split volume
+
+ const SMDS_MeshNode** volNodes = volTool.GetNodes();
+ for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
+ {
+ const int nbNodes = volTool.NbFaceNodes( iF ) / iQ;
+ if ( nbNodes < 4 ) continue;
+
+ // find an existing face
+ vector<const SMDS_MeshNode*> fNodes( volTool.GetFaceNodes( iF ),
+ volTool.GetFaceNodes( iF ) + nbNodes*iQ );
+ while ( const SMDS_MeshElement* face = GetMeshDS()->FindFace( fNodes ))
+ {
+ // among possible triangles create ones discribed by split method
+ const int* nInd = volTool.GetFaceNodesIndices( iF );
+ int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
+ int iCom = 0; // common node of triangle faces to split into
+ list< TTriangleFacet > facets;
+ for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCom )
+ {
+ TTriangleFacet t012( nInd[ iQ * ( iCom )],
+ nInd[ iQ * ( (iCom+1)%nbNodes )],
+ nInd[ iQ * ( (iCom+2)%nbNodes )]);
+ TTriangleFacet t023( nInd[ iQ * ( iCom )],
+ nInd[ iQ * ( (iCom+2)%nbNodes )],
+ nInd[ iQ * ( (iCom+3)%nbNodes )]);
+ if ( splitMethod.hasFacet( t012 ) && splitMethod.hasFacet( t023 ))
+ {
+ facets.push_back( t012 );
+ facets.push_back( t023 );
+ for ( int iLast = iCom+4; iLast < iCom+nbNodes; ++iLast )
+ facets.push_back( TTriangleFacet( nInd[ iQ * ( iCom )],
+ nInd[ iQ * ((iLast-1)%nbNodes )],
+ nInd[ iQ * ((iLast )%nbNodes )]));
+ break;
+ }
+ }
+ // find submesh to add new faces in
+ if ( !fSubMesh || !fSubMesh->Contains( face ))
+ {
+ int shapeID = FindShape( face );
+ fSubMesh = GetMeshDS()->MeshElements( shapeID );
+ }
+ // make triangles
+ helper.SetElementsOnShape( false );
+ vector< const SMDS_MeshElement* > triangles;
+ list< TTriangleFacet >::iterator facet = facets.begin();
+ for ( ; facet != facets.end(); ++facet )
+ {
+ if ( !volTool.IsFaceExternal( iF ))
+ swap( facet->_n2, facet->_n3 );
+ triangles.push_back( helper.AddFace( volNodes[ facet->_n1 ],
+ volNodes[ facet->_n2 ],
+ volNodes[ facet->_n3 ]));
+ if ( triangles.back() && fSubMesh )
+ fSubMesh->AddElement( triangles.back());
+ newElems.Append( triangles.back() );
+ }
+ ReplaceElemInGroups( face, triangles, GetMeshDS() );
+ GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false );
+ }
+
+ } // loop on volume faces to split them into triangles
+
+ GetMeshDS()->RemoveFreeElement( *elem, subMesh, /*fromGroups=*/false );
+
+ } // loop on volumes to split
+
+ myLastCreatedNodes = newNodes;
+ myLastCreatedElems = newElems;
+}