-// Copyright (C) 2007-2008 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE
//
// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
+
// SMESH SMESH : idl implementation based on 'SMESH' unit's classes
// File : SMESH_MeshEditor.cxx
// Created : Mon Apr 12 16:10:22 2004
int nbnode = node.size();
SMESHDS_Mesh* mesh = GetMeshDS();
switch ( type ) {
+ case SMDSAbs_0DElement:
+ if ( nbnode == 1 )
+ if ( ID ) e = mesh->Add0DElementWithID(node[0], ID);
+ else e = mesh->Add0DElement (node[0] );
+ break;
case SMDSAbs_Edge:
if ( nbnode == 2 )
if ( ID ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
{
// Methods of splitting volumes into tetra
- const int theHexTo5[5*4] =
+ const int theHexTo5_1[5*4+1] =
+ {
+ 0, 1, 2, 5, 0, 4, 5, 7, 0, 2, 3, 7, 2, 5, 6, 7, 0, 5, 2, 7, -1
+ };
+ const int theHexTo5_2[5*4+1] =
+ {
+ 1, 2, 3, 6, 1, 4, 5, 6, 0, 1, 3, 4, 3, 4, 6, 7, 1, 3, 4, 6, -1
+ };
+ const int* theHexTo5[2] = { theHexTo5_1, theHexTo5_2 };
+
+ const int theHexTo6_1[6*4+1] =
{
- 0, 1, 5, 2,
- 0, 4, 5, 7,
- 0, 3, 7, 2,
- 5, 6, 7, 2,
- 0, 2, 5, 7
+ 1, 5, 6, 0, 0, 1, 2, 6, 0, 4, 5, 6, 0, 4, 6, 7, 0, 2, 3, 6, 0, 3, 7, 6, -1
};
- const int theHexTo6[6*4] =
+ const int theHexTo6_2[6*4+1] =
{
- 0, 1, 5, 2,
- 0, 4, 5, 7,
- 0, 3, 7, 2,
- 5, 6, 7, 2,
- 0, 2, 5, 7
+ 2, 6, 7, 1, 1, 2, 3, 7, 1, 5, 6, 7, 1, 5, 7, 4, 1, 3, 0, 7, 1, 0, 4, 7, -1
};
- const int thePyraTo2[2*4] =
+ const int theHexTo6_3[6*4+1] =
{
- 0, 1, 2, 4,
- 0, 2, 3, 4
+ 3, 7, 4, 2, 2, 3, 0, 4, 2, 6, 7, 4, 2, 6, 4, 5, 2, 0, 1, 4, 2, 1, 5, 4, -1
};
+ const int theHexTo6_4[6*4+1] =
+ {
+ 0, 4, 5, 3, 3, 0, 1, 5, 3, 7, 4, 5, 3, 7, 5, 6, 3, 1, 2, 5, 3, 2, 6, 5, -1
+ };
+ const int* theHexTo6[4] = { theHexTo6_1, theHexTo6_2, theHexTo6_3, theHexTo6_4 };
- const int thePentaTo8[8*4] =
+ const int thePyraTo2_1[2*4+1] =
+ {
+ 0, 1, 2, 4, 0, 2, 3, 4, -1
+ };
+ const int thePyraTo2_2[2*4+1] =
{
- 0, 1, 2, 6,
- 3, 5, 4, 6,
- 0, 3, 4, 6,
- 0, 4, 1, 6,
- 1, 4, 5, 6,
- 1, 5, 2, 6,
- 2, 5, 3, 6,
- 2, 3, 0, 6
+ 1, 2, 3, 4, 1, 3, 0, 4, -1
};
+ const int* thePyraTo2[2] = { thePyraTo2_1, thePyraTo2_2 };
+ const int thePentaTo3_1[3*4+1] =
+ {
+ 0, 1, 2, 3, 1, 3, 4, 2, 2, 3, 4, 5, -1
+ };
+ const int thePentaTo3_2[3*4+1] =
+ {
+ 1, 2, 0, 4, 2, 4, 5, 0, 0, 4, 5, 3, -1
+ };
+ const int thePentaTo3_3[3*4+1] =
+ {
+ 2, 0, 1, 5, 0, 5, 3, 1, 1, 5, 3, 4, -1
+ };
+ const int thePentaTo3_4[3*4+1] =
+ {
+ 0, 1, 2, 3, 1, 3, 4, 5, 2, 3, 1, 5, -1
+ };
+ const int thePentaTo3_5[3*4+1] =
+ {
+ 1, 2, 0, 4, 2, 4, 5, 3, 0, 4, 2, 3, -1
+ };
+ const int thePentaTo3_6[3*4+1] =
+ {
+ 2, 0, 1, 5, 0, 5, 3, 4, 1, 5, 0, 4, -1
+ };
+ const int* thePentaTo3[6] = { thePentaTo3_1, thePentaTo3_2, thePentaTo3_3,
+ thePentaTo3_4, thePentaTo3_5, thePentaTo3_6 };
+
+ struct TTriangleFacet //!< stores indices of three nodes of tetra facet
+ {
+ int _n1, _n2, _n3;
+ TTriangleFacet(int n1, int n2, int n3): _n1(n1), _n2(n2), _n3(n3) {}
+ bool contains(int n) const { return ( n == _n1 || n == _n2 || n == _n3 ); }
+ bool hasAdjacentTetra( const SMDS_MeshElement* elem ) const;
+ };
struct TSplitMethod
{
int _nbTetra;
- const int* _connectivity;
- bool _addNode; // additional node is to be created
+ const int* _connectivity; //!< foursomes of tetra connectivy finished by -1
+ bool _baryNode; //!< additional node is to be created at cell barycenter
+ bool _ownConn; //!< to delete _connectivity in destructor
+
TSplitMethod( int nbTet=0, const int* conn=0, bool addNode=false)
- : _nbTetra(nbTet), _connectivity(conn), _addNode(addNode) {}
+ : _nbTetra(nbTet), _connectivity(conn), _baryNode(addNode), _ownConn(false) {}
+ ~TSplitMethod() { if ( _ownConn ) delete [] _connectivity; _connectivity = 0; }
+ bool hasFacet( const TTriangleFacet& facet ) const
+ {
+ const int* tetConn = _connectivity;
+ for ( ; tetConn[0] >= 0; tetConn += 4 )
+ if (( facet.contains( tetConn[0] ) +
+ facet.contains( tetConn[1] ) +
+ facet.contains( tetConn[2] ) +
+ facet.contains( tetConn[3] )) == 3 )
+ return true;
+ return false;
+ }
};
+ //=======================================================================
/*!
* \brief return TSplitMethod for the given element
*/
- TSplitMethod getSplitMethod( const SMDS_MeshElement* vol, const int theMethodFlags)
+ //=======================================================================
+
+ TSplitMethod getSplitMethod( SMDS_VolumeTool& vol, const int theMethodFlags)
{
+ int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
+
+ // Find out how adjacent volumes are split
+
+ vector < list< TTriangleFacet > > triaSplitsByFace( vol.NbFaces() ); // splits of each side
+ int hasAdjacentSplits = 0, maxTetConnSize = 0;
+ for ( int iF = 0; iF < vol.NbFaces(); ++iF )
+ {
+ int nbNodes = vol.NbFaceNodes( iF ) / iQ;
+ maxTetConnSize += 4 * ( nbNodes - 2 );
+ if ( nbNodes < 4 ) continue;
+
+ list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
+ const int* nInd = vol.GetFaceNodesIndices( iF );
+ if ( nbNodes == 4 )
+ {
+ TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
+ TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
+ if ( t012.hasAdjacentTetra( vol.Element() )) triaSplits.push_back( t012 );
+ else if ( t123.hasAdjacentTetra( vol.Element() )) triaSplits.push_back( t123 );
+ }
+ else
+ {
+ int iCom = 0; // common node of triangle faces to split into
+ for ( int iVar = 0; iVar < nbNodes; ++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 ( t012.hasAdjacentTetra( vol.Element() ) && t023.hasAdjacentTetra( vol.Element() ))
+ {
+ triaSplits.push_back( t012 );
+ triaSplits.push_back( t023 );
+ break;
+ }
+ }
+ }
+ if ( !triaSplits.empty() )
+ hasAdjacentSplits = true;
+ }
+
+ // Among variants of split method select one compliant with adjacent volumes
+
TSplitMethod method;
- if ( vol->GetType() == SMDSAbs_Volume && !vol->IsPoly())
- switch ( vol->NbNodes() )
+ if ( !vol.Element()->IsPoly() )
+ {
+ int nbVariants = 2, nbTet = 0;
+ const int** connVariants = 0;
+ switch ( vol.Element()->GetEntityType() )
{
- case 8:
- case 20:
+ case SMDSEntity_Hexa:
+ case SMDSEntity_Quad_Hexa:
if ( theMethodFlags & SMESH_MeshEditor::HEXA_TO_5 )
- method = TSplitMethod( 5, theHexTo5 );
+ connVariants = theHexTo5, nbTet = 5;
else
- method = TSplitMethod( 6, theHexTo6 );
+ connVariants = theHexTo6, nbTet = 6, nbVariants = 4;
break;
- case 5:
- case 13:
- method = TSplitMethod( 2, thePyraTo2 );
+ case SMDSEntity_Pyramid:
+ case SMDSEntity_Quad_Pyramid:
+ connVariants = thePyraTo2; nbTet = 2;
break;
- case 6:
- case 15:
- method = TSplitMethod( 8, thePentaTo8, /*addNode=*/true );
+ case SMDSEntity_Penta:
+ case SMDSEntity_Quad_Penta:
+ connVariants = thePentaTo3; nbTet = 3; nbVariants = 6;
break;
- default:;
+ default:
+ nbVariants = 0;
+ }
+ for ( int variant = 0; variant < nbVariants && method._nbTetra == 0; ++variant )
+ {
+ // check method compliancy with adjacent tetras,
+ // all found splits must be among facets of tetras described by this method
+ method = TSplitMethod( nbTet, connVariants[variant] );
+ if ( hasAdjacentSplits && method._nbTetra > 0 )
+ {
+ bool facetCreated = true;
+ for ( int iF = 0; facetCreated && iF < triaSplitsByFace.size(); ++iF )
+ {
+ list< TTriangleFacet >::const_iterator facet = triaSplitsByFace[iF].begin();
+ for ( ; facetCreated && facet != triaSplitsByFace[iF].end(); ++facet )
+ facetCreated = method.hasFacet( *facet );
+ }
+ if ( !facetCreated )
+ method = TSplitMethod(0); // incompatible method
+ }
}
+ }
+ if ( method._nbTetra < 1 )
+ {
+ // No standard method is applicable, use a generic solution:
+ // each facet of a volume is split into triangles and
+ // each of triangles and a volume barycenter form a tetrahedron.
+
+ int* connectivity = new int[ maxTetConnSize + 1 ];
+ method._connectivity = connectivity;
+ method._ownConn = true;
+ method._baryNode = true;
+
+ int connSize = 0;
+ int baryCenInd = vol.NbNodes();
+ for ( int iF = 0; iF < vol.NbFaces(); ++iF )
+ {
+ const int nbNodes = vol.NbFaceNodes( iF ) / iQ;
+ const int* nInd = vol.GetFaceNodesIndices( iF );
+ // find common node of triangle facets of tetra to create
+ int iCommon = 0; // index in linear numeration
+ const list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
+ if ( !triaSplits.empty() )
+ {
+ // by found facets
+ const TTriangleFacet* facet = &triaSplits.front();
+ for ( ; iCommon < nbNodes-1 ; ++iCommon )
+ if ( facet->contains( nInd[ iQ * iCommon ]) &&
+ facet->contains( nInd[ iQ * ((iCommon+2)%nbNodes) ]))
+ break;
+ }
+ else if ( nbNodes > 3 )
+ {
+ // find the best method of splitting into triangles by aspect ratio
+ SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio);
+ map< double, int > badness2iCommon;
+ const SMDS_MeshNode** nodes = vol.GetFaceNodes( iF );
+ int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
+ for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCommon )
+ for ( int iLast = iCommon+2; iLast < iCommon+nbNodes; ++iLast )
+ {
+ SMDS_FaceOfNodes tria ( nodes[ iQ*( iCommon )],
+ nodes[ iQ*((iLast-1)%nbNodes)],
+ nodes[ iQ*((iLast )%nbNodes)]);
+ double badness = getBadRate( &tria, aspectRatio );
+ badness2iCommon.insert( make_pair( badness, iCommon ));
+ }
+ // use iCommon with lowest badness
+ iCommon = badness2iCommon.begin()->second;
+ }
+ if ( iCommon >= nbNodes )
+ iCommon = 0; // something wrong
+ // fill connectivity of tetra
+ int nbTet = nbNodes - 2;
+ for ( int i = 0; i < nbTet; ++i )
+ {
+ int i1 = (iCommon+1+i) % nbNodes, i2 = (iCommon+2+i) % nbNodes;
+ if ( !vol.IsFaceExternal( iF )) swap( i1, i2 );
+ connectivity[ connSize++ ] = nInd[ iQ * iCommon ];
+ connectivity[ connSize++ ] = nInd[ iQ * i1 ];
+ connectivity[ connSize++ ] = nInd[ iQ * i2 ];
+ connectivity[ connSize++ ] = baryCenInd;
+ ++method._nbTetra;
+ }
+ }
+ connectivity[ connSize++ ] = -1;
+ }
return method;
}
-}
+ //================================================================================
+ /*!
+ * \brief Check if there is a tetraherdon adjacent to the given element via this facet
+ */
+ //================================================================================
+
+ bool TTriangleFacet::hasAdjacentTetra( const SMDS_MeshElement* elem ) const
+ {
+ // find the tetrahedron including the three nodes of facet
+ const SMDS_MeshNode* n1 = elem->GetNode(_n1);
+ const SMDS_MeshNode* n2 = elem->GetNode(_n2);
+ const SMDS_MeshNode* n3 = elem->GetNode(_n3);
+ SMDS_ElemIteratorPtr volIt1 = n1->GetInverseElementIterator(SMDSAbs_Volume);
+ while ( volIt1->more() )
+ {
+ const SMDS_MeshElement* v = volIt1->next();
+ if ( v->GetEntityType() != ( v->IsQuadratic() ? SMDSEntity_Quad_Tetra : SMDSEntity_Tetra ))
+ continue;
+ SMDS_ElemIteratorPtr volIt2 = n2->GetInverseElementIterator(SMDSAbs_Volume);
+ while ( volIt2->more() )
+ if ( v != volIt2->next() )
+ continue;
+ SMDS_ElemIteratorPtr volIt3 = n3->GetInverseElementIterator(SMDSAbs_Volume);
+ while ( volIt3->more() )
+ if ( v == volIt3->next() )
+ return true;
+ }
+ return false;
+ }
+} // namespace
//=======================================================================
//function : SplitVolumesIntoTetra
void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
const int theMethodFlags)
{
- // sdt-like iterator on coordinates of nodes of mesh element
+ // 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 edge
+ 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
- SMDS_VolumeTool vol( *elem );
- for ( int iF = 0; iF < vol.NbFaces(); ++iF )
+ for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
{
- const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iF );
- for ( int iN = 0; iN < vol.NbFaceNodes( iF ); iN += 2)
+ 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 );
+ }
- vector<const SMDS_MeshElement* > tetras; // splits of a volume
+ // 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] ]));
- if ( geomType == SMDSEntity_Polyhedra )
- {
- // Each face of a polyhedron is split into triangles and
- // each of triangles and a cell barycenter form a tetrahedron.
+ ReplaceElemInGroups( *elem, tetras, GetMeshDS() );
- SMDS_VolumeTool vol( *elem );
+ // Split faces on sides of the split volume
- // 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() );
+ 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;
- for ( int iF = 0; iF < vol.NbFaces(); ++iF )
+ // 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 ))
{
- 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 )
+ // 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 )
{
- 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 );
+ 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 );
}
- }
- else
- {
+ } // loop on volume faces to split them into triangles
- TSplitMethod splitMethod = getSplitMethod( *elem, theMethodFlags );
- if ( splitMethod._nbTetra < 1 ) continue;
+ GetMeshDS()->RemoveFreeElement( *elem, subMesh, /*fromGroups=*/false );
- vector<const SMDS_MeshNode*> volNodes( (*elem)->begin_nodes(), (*elem)->end_nodes());
- }
- }
+ } // loop on volumes to split
+
+ myLastCreatedNodes = newNodes;
+ myLastCreatedElems = newElems;
}
//=======================================================================
}
}
-//=======================================================================
-//function : ReplaceElemInGroups
-//purpose : replace elemToRm by elemToAdd in the all groups
-//=======================================================================
+//================================================================================
+/*!
+ * \brief Replace elemToRm by elemToAdd in the all groups
+ */
+//================================================================================
void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
const SMDS_MeshElement* elemToAdd,
}
}
+//================================================================================
+/*!
+ * \brief Replace elemToRm by elemToAdd in the all groups
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
+ const vector<const SMDS_MeshElement*>& elemToAdd,
+ SMESHDS_Mesh * aMesh)
+{
+ const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
+ if (!groups.empty())
+ {
+ set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
+ for ( ; grIt != groups.end(); grIt++ ) {
+ SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
+ if ( group && group->SMDSGroup().Remove( elemToRm ) )
+ for ( int i = 0; i < elemToAdd.size(); ++i )
+ group->SMDSGroup().Add( elemToAdd[ i ] );
+ }
+ }
+}
+
//=======================================================================
//function : QuadToTri
//purpose : Cut quadrangles into triangles.
SMESH_NodeSearcherImpl* _nodeSearcher;
SMDSAbs_ElementType _elementType;
double _tolerance;
- set<const SMDS_MeshElement*> _internalFaces;
+ bool _outerFacesFound;
+ set<const SMDS_MeshElement*> _outerFaces; // empty means "no internal faces at all"
SMESH_ElementSearcherImpl( SMESHDS_Mesh& mesh )
- : _mesh(&mesh),_ebbTree(0),_nodeSearcher(0), _tolerance(-1) {}
+ : _mesh(&mesh),_ebbTree(0),_nodeSearcher(0), _tolerance(-1), _outerFacesFound(false) {}
~SMESH_ElementSearcherImpl()
{
if ( _ebbTree ) delete _ebbTree; _ebbTree = 0;
vector< const SMDS_MeshElement* >& foundElements);
virtual TopAbs_State GetPointState(const gp_Pnt& point);
- struct TInters //!< data of intersection of the line and the mesh face
+ double getTolerance();
+ bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face,
+ const double tolerance, double & param);
+ void findOuterBoundary(const SMDS_MeshElement* anyOuterFace);
+ bool isOuterBoundary(const SMDS_MeshElement* face) const
+ {
+ return _outerFaces.empty() || _outerFaces.count(face);
+ }
+ struct TInters //!< data of intersection of the line and the mesh face used in GetPointState()
{
const SMDS_MeshElement* _face;
gp_Vec _faceNorm;
TInters(const SMDS_MeshElement* face, const gp_Vec& faceNorm, bool coinc=false)
: _face(face), _faceNorm( faceNorm ), _coincides( coinc ) {}
};
- double getTolerance();
- bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face,
- const double tolerance, double & param);
- void findOuterBoundary();
- bool isOuterBoundary(const SMDS_MeshElement* face) const { return !_internalFaces.count(face);}
+ struct TFaceLink //!< link and faces sharing it (used in findOuterBoundary())
+ {
+ SMESH_TLink _link;
+ TIDSortedElemSet _faces;
+ TFaceLink( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshElement* face)
+ : _link( n1, n2 ), _faces( &face, &face + 1) {}
+ };
};
+ostream& operator<< (ostream& out, const SMESH_ElementSearcherImpl::TInters& i)
+{
+ return out << "TInters(face=" << ( i._face ? i._face->GetID() : 0)
+ << ", _coincides="<<i._coincides << ")";
+}
+
//=======================================================================
/*!
* \brief define tolerance for search
*/
//================================================================================
-void SMESH_ElementSearcherImpl::findOuterBoundary()
+void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerFace)
{
-
+ if ( _outerFacesFound ) return;
+
+ // Collect all outer faces by passing from one outer face to another via their links
+ // and BTW find out if there are internal faces at all.
+
+ // checked links and links where outer boundary meets internal one
+ set< SMESH_TLink > visitedLinks, seamLinks;
+
+ // links to treat with already visited faces sharing them
+ list < TFaceLink > startLinks;
+
+ // load startLinks with the first outerFace
+ startLinks.push_back( TFaceLink( outerFace->GetNode(0), outerFace->GetNode(1), outerFace));
+ _outerFaces.insert( outerFace );
+
+ TIDSortedElemSet emptySet;
+ while ( !startLinks.empty() )
+ {
+ const SMESH_TLink& link = startLinks.front()._link;
+ TIDSortedElemSet& faces = startLinks.front()._faces;
+
+ outerFace = *faces.begin();
+ // find other faces sharing the link
+ const SMDS_MeshElement* f;
+ while (( f = SMESH_MeshEditor::FindFaceInSet(link.node1(), link.node2(), emptySet, faces )))
+ faces.insert( f );
+
+ // select another outer face among the found
+ const SMDS_MeshElement* outerFace2 = 0;
+ if ( faces.size() == 2 )
+ {
+ outerFace2 = (outerFace == *faces.begin() ? *faces.rbegin() : *faces.begin());
+ }
+ else if ( faces.size() > 2 )
+ {
+ seamLinks.insert( link );
+
+ // link direction within the outerFace
+ gp_Vec n1n2( SMESH_MeshEditor::TNodeXYZ( link.node1()),
+ SMESH_MeshEditor::TNodeXYZ( link.node2()));
+ int i1 = outerFace->GetNodeIndex( link.node1() );
+ int i2 = outerFace->GetNodeIndex( link.node2() );
+ bool rev = ( abs(i2-i1) == 1 ? i1 > i2 : i2 > i1 );
+ if ( rev ) n1n2.Reverse();
+ // outerFace normal
+ gp_XYZ ofNorm, fNorm;
+ if ( SMESH_Algo::FaceNormal( outerFace, ofNorm, /*normalized=*/false ))
+ {
+ // direction from the link inside outerFace
+ gp_Vec dirInOF = gp_Vec( ofNorm ) ^ n1n2;
+ // sort all other faces by angle with the dirInOF
+ map< double, const SMDS_MeshElement* > angle2Face;
+ set< const SMDS_MeshElement*, TIDCompare >::const_iterator face = faces.begin();
+ for ( ; face != faces.end(); ++face )
+ {
+ if ( !SMESH_Algo::FaceNormal( *face, fNorm, /*normalized=*/false ))
+ continue;
+ gp_Vec dirInF = gp_Vec( fNorm ) ^ n1n2;
+ double angle = dirInOF.AngleWithRef( dirInF, n1n2 );
+ if ( angle < 0 ) angle += 2*PI;
+ angle2Face.insert( make_pair( angle, *face ));
+ }
+ if ( !angle2Face.empty() )
+ outerFace2 = angle2Face.begin()->second;
+ }
+ }
+ // store the found outer face and add its links to continue seaching from
+ if ( outerFace2 )
+ {
+ _outerFaces.insert( outerFace );
+ int nbNodes = outerFace2->NbNodes()/( outerFace2->IsQuadratic() ? 2 : 1 );
+ for ( int i = 0; i < nbNodes; ++i )
+ {
+ SMESH_TLink link2( outerFace2->GetNode(i), outerFace2->GetNode((i+1)%nbNodes));
+ if ( visitedLinks.insert( link2 ).second )
+ startLinks.push_back( TFaceLink( link2.node1(), link2.node2(), outerFace2 ));
+ }
+ }
+ startLinks.pop_front();
+ }
+ _outerFacesFound = true;
+
+ if ( !seamLinks.empty() )
+ {
+ // There are internal boundaries touching the outher one,
+ // find all faces of internal boundaries in order to find
+ // faces of boundaries of holes, if any.
+
+ }
+ else
+ {
+ _outerFaces.clear();
+ }
}
//=======================================================================
if ( _ebbTree ) delete _ebbTree;
_ebbTree = new ElementBndBoxTree( *_mesh, _elementType = SMDSAbs_Face );
}
- // algo: analyse transition of a line starting at the point through mesh boundary;
- // try several lines, if none of attemps gives a clear answer, we give up as the
- // task can be too complex including internal boundaries, concave surfaces etc.
+ // Algo: analyse transition of a line starting at the point through mesh boundary;
+ // try three lines parallel to axis of the coordinate system and perform rough
+ // analysis. If solution is not clear perform thorough analysis.
const int nbAxes = 3;
gp_Dir axisDir[ nbAxes ] = { gp::DX(), gp::DY(), gp::DZ() };
if ( !SMESH_Algo::FaceNormal( *face, fNorm, /*normalized=*/false)) continue;
gp_Pln facePlane( SMESH_MeshEditor::TNodeXYZ( (*face)->GetNode(0)), fNorm );
- // intersection
+ // perform intersection
IntAna_IntConicQuad intersection( line, IntAna_Quadric( facePlane ));
if ( !intersection.IsDone() )
continue;
return TopAbs_OUT;
double f = u2inters.begin()->first, l = u2inters.rbegin()->first;
- if ( nbInter == 1 )
+ if ( nbInter == 1 ) // not closed mesh
return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN;
if ( fabs( f ) < tolerance || fabs( l ) < tolerance )
nbInt2Axis.insert( make_pair( min( nbIntBeforePoint, nbIntAfterPoint ), axis ));
+ if ( _outerFacesFound ) break; // pass to thorough analysis
+
} // three attempts - loop on CS axes
- // Analyse intersections thoroughly
- // We make two loops, on the first one we correctly exclude touching intersections,
- // on the second, we additionally just throw away intersections with small angles
+ // Analyse intersections thoroughly.
+ // We make two loops maximum, on the first one we only exclude touching intersections,
+ // on the second, if situation is still unclear, we gather and use information on
+ // position of faces (internal or outer). If faces position is already gathered,
+ // we make the second loop right away.
- for ( int angleCheck = 0; angleCheck < 2; ++angleCheck )
+ for ( int hasPositionInfo = _outerFacesFound; hasPositionInfo < 2; ++hasPositionInfo )
{
multimap< int, int >::const_iterator nb_axis = nbInt2Axis.begin();
for ( ; nb_axis != nbInt2Axis.end(); ++nb_axis )
tangentInters[ axis ].clear();
// Count intersections before and after the point excluding touching ones.
+ // If hasPositionInfo we count intersections of outer boundary only
int nbIntBeforePoint = 0, nbIntAfterPoint = 0;
double f = numeric_limits<double>::max(), l = -numeric_limits<double>::max();
- map< double, TInters >::iterator u_int2 = u2inters.begin(), u_int1 = u_int2++;
+ map< double, TInters >::iterator u_int1 = u2inters.begin(), u_int2 = u_int1;
bool ok = ! u_int1->second._coincides;
while ( ok && u_int1 != u2inters.end() )
{
- // skip intersections at the same point (if line pass through edge or node)
- int nbSamePnt = 0;
double u = u_int1->first;
- while ( u_int2 != u2inters.end() && fabs( u_int2->first - u ) < tolerance )
+ bool touchingInt = false;
+ if ( ++u_int2 != u2inters.end() )
{
- ++nbSamePnt;
- ++u_int2;
- }
+ // skip intersections at the same point (if the line passes through edge or node)
+ int nbSamePnt = 0;
+ while ( u_int2 != u2inters.end() && fabs( u_int2->first - u ) < tolerance )
+ {
+ ++nbSamePnt;
+ ++u_int2;
+ }
- // skip tangent intersections
- int nbTgt = 0;
- const SMDS_MeshElement* prevFace = u_int1->second._face;
- while ( ok && u_int2->second._coincides )
- {
- if ( SMESH_Algo::GetCommonNodes(prevFace , u_int2->second._face).empty() )
- ok = false;
- else
+ // skip tangent intersections
+ int nbTgt = 0;
+ const SMDS_MeshElement* prevFace = u_int1->second._face;
+ while ( ok && u_int2->second._coincides )
{
- nbTgt++;
- u_int2++;
- ok = ( u_int2 != u2inters.end() );
+ if ( SMESH_Algo::GetCommonNodes(prevFace , u_int2->second._face).empty() )
+ ok = false;
+ else
+ {
+ nbTgt++;
+ u_int2++;
+ ok = ( u_int2 != u2inters.end() );
+ }
}
- }
- if ( !ok ) break;
+ if ( !ok ) break;
- // skip intersections at the same point after tangent intersections
- if ( nbTgt > 0 )
- {
- double u = u_int2->first;
- ++u_int2;
- while ( u_int2 != u2inters.end() && fabs( u_int2->first - u ) < tolerance )
+ // skip intersections at the same point after tangent intersections
+ if ( nbTgt > 0 )
{
- ++nbSamePnt;
+ double u2 = u_int2->first;
++u_int2;
+ while ( u_int2 != u2inters.end() && fabs( u_int2->first - u2 ) < tolerance )
+ {
+ ++nbSamePnt;
+ ++u_int2;
+ }
}
- }
-
- bool touchingInt = false;
- if ( nbSamePnt + nbTgt > 0 )
- {
- double minDot = numeric_limits<double>::max(), maxDot = -numeric_limits<double>::max();
- map< double, TInters >::iterator u_int = u_int1;
- for ( ; u_int != u_int2; ++u_int )
+ // decide if we skipped a touching intersection
+ if ( nbSamePnt + nbTgt > 0 )
{
- if ( u_int->second._coincides ) continue;
- double dot = u_int->second._faceNorm * line.Direction();
- if ( dot > maxDot ) maxDot = dot;
- if ( dot < minDot ) minDot = dot;
+ double minDot = numeric_limits<double>::max(), maxDot = -numeric_limits<double>::max();
+ map< double, TInters >::iterator u_int = u_int1;
+ for ( ; u_int != u_int2; ++u_int )
+ {
+ if ( u_int->second._coincides ) continue;
+ double dot = u_int->second._faceNorm * line.Direction();
+ if ( dot > maxDot ) maxDot = dot;
+ if ( dot < minDot ) minDot = dot;
+ }
+ touchingInt = ( minDot*maxDot < 0 );
}
- touchingInt = ( minDot*maxDot < 0 );
- }
- // throw away intersection with lower angles
- if ( !touchingInt && angleCheck )
- {
- const double angTol = 2 * Standard_PI180, normAng = Standard_PI / 2;
- double angle = u_int1->second._faceNorm.Angle( line.Direction() );
- touchingInt = ( fabs( angle - normAng ) < angTol );
}
if ( !touchingInt )
{
- if ( u < 0 )
- ++nbIntBeforePoint;
- else
- ++nbIntAfterPoint;
-
+ if ( !hasPositionInfo || isOuterBoundary( u_int1->second._face ))
+ {
+ if ( u < 0 )
+ ++nbIntBeforePoint;
+ else
+ ++nbIntAfterPoint;
+ }
if ( u < f ) f = u;
if ( u > l ) l = u;
}
- u_int1 = u_int2++; // to next intersection
+ u_int1 = u_int2; // to next intersection
} // loop on intersections with one line
if ( ok )
{
+ if ( fabs( f ) < tolerance || fabs( l ) < tolerance )
+ return TopAbs_ON;
+
if ( nbIntBeforePoint == 0 || nbIntAfterPoint == 0)
return TopAbs_OUT;
- if ( nbIntBeforePoint + nbIntAfterPoint == 1 )
+ if ( nbIntBeforePoint + nbIntAfterPoint == 1 ) // not closed mesh
return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN;
if ( nbIntBeforePoint == 1 || nbIntAfterPoint == 1 )
- return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN;
-
- if ( fabs( f ) < tolerance || fabs( l ) < tolerance )
- return TopAbs_ON;
+ return TopAbs_IN;
if ( (f<0) == (l<0) )
return TopAbs_OUT;
+
+ if ( hasPositionInfo )
+ return nbIntBeforePoint % 2 ? TopAbs_IN : TopAbs_OUT;
}
} // loop on intersections of the tree lines - thorough analysis
- } // two attempts - with and w/o angleCheck
+
+ if ( !hasPositionInfo )
+ {
+ // gather info on faces position - is face in the outer boundary or not
+ map< double, TInters > & u2inters = paramOnLine2TInters[ 0 ];
+ findOuterBoundary( u2inters.begin()->second._face );
+ }
+
+ } // two attempts - with and w/o faces position info in the mesh
return TopAbs_UNKNOWN;
}
SMESHDS_Mesh* aMesh = GetMeshDS();
if (!aMesh)
return false;
- bool res = false;
+ //bool res = false;
+ int nbFree = 0, nbExisted = 0, nbCreated = 0;
SMDS_VolumeIteratorPtr vIt = aMesh->volumesIterator();
while(vIt->more())
{
{
if (!vTool.IsFreeFace(iface))
continue;
+ nbFree++;
vector<const SMDS_MeshNode *> nodes;
int nbFaceNodes = vTool.NbFaceNodes(iface);
const SMDS_MeshNode** faceNodes = vTool.GetFaceNodes(iface);
nodes.push_back(faceNodes[inode]);
// add new face based on volume nodes
- if (aMesh->FindFace( nodes ) )
+ if (aMesh->FindFace( nodes ) ) {
+ nbExisted++;
continue; // face already exsist
+ }
myLastCreatedElems.Append( AddElement(nodes, SMDSAbs_Face, isPoly && iface == 1) );
- res = true;
+ nbCreated++;
}
}
- return res;
+ return ( nbFree==(nbExisted+nbCreated) );
}
