-// Copyright (C) 2007-2011 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012 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
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// 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
// Author : Edward AGAPOV (eap)
#include <TopTools_SequenceOfShape.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Face.hxx>
+#include <TopoDS_Solid.hxx>
#include <gp.hxx>
#include <gp_Ax1.hxx>
#include <gp_Dir.hxx>
#include <algorithm>
#include <sstream>
+#include <boost/tuple/tuple.hpp>
+
+#include <Standard_Failure.hxx>
+#include <Standard_ErrorHandler.hxx>
+
#define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
using namespace std;
{
}
+//================================================================================
+/*!
+ * \brief Clears myLastCreatedNodes and myLastCreatedElems
+ */
+//================================================================================
+
+void SMESH_MeshEditor::CrearLastCreated()
+{
+ myLastCreatedNodes.Clear();
+ myLastCreatedElems.Clear();
+}
+
+
//=======================================================================
/*!
* \brief Add element
SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
const SMDSAbs_ElementType type,
const bool isPoly,
- const int ID)
+ const int ID,
+ const double ballDiameter)
{
//MESSAGE("AddElement " <<node.size() << " " << type << " " << isPoly << " " << ID);
SMDS_MeshElement* e = 0;
else e = mesh->AddNode (node[0]->X(), node[0]->Y(), node[0]->Z());
break;
+ case SMDSAbs_Ball:
+ if ( ID >= 1 ) e = mesh->AddBallWithID(node[0], ballDiameter, ID);
+ else e = mesh->AddBall (node[0], ballDiameter);
+ break;
+
default:;
}
if ( e ) myLastCreatedElems.Append( e );
return removed;
}
+//================================================================================
+/*!
+ * \brief Create 0D elements on all nodes of the given object except those
+ * nodes on which a 0D element already exists.
+ * \param elements - Elements on whose nodes to create 0D elements; if empty,
+ * the all mesh is treated
+ * \param all0DElems - returns all 0D elements found or created on nodes of \a elements
+ */
+//================================================================================
+
+void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& elements,
+ TIDSortedElemSet& all0DElems )
+{
+ typedef SMDS_SetIterator<const SMDS_MeshElement*, TIDSortedElemSet::const_iterator> TSetIterator;
+ SMDS_ElemIteratorPtr elemIt;
+ if ( elements.empty() )
+ elemIt = GetMeshDS()->elementsIterator( SMDSAbs_Node );
+ else
+ elemIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
+
+ while ( elemIt->more() )
+ {
+ const SMDS_MeshElement* e = elemIt->next();
+ SMDS_ElemIteratorPtr nodeIt = e->nodesIterator();
+ while ( nodeIt->more() )
+ {
+ const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
+ SMDS_ElemIteratorPtr it0D = n->GetInverseElementIterator( SMDSAbs_0DElement );
+ if ( it0D->more() )
+ all0DElems.insert( it0D->next() );
+ else {
+ myLastCreatedElems.Append( GetMeshDS()->Add0DElement( n ));
+ all0DElems.insert( myLastCreatedElems.Last() );
+ }
+ }
+ }
+}
+
//=======================================================================
//function : FindShape
//purpose : Return an index of the shape theElem is on
//=======================================================================
//function : edgeConnectivity
-//purpose : auxilary
+//purpose : auxilary
// return number of the edges connected with the theNode.
// if theEdges has connections with the other type of the
-// elements, return -1
+// elements, return -1
//=======================================================================
static int nbEdgeConnectivity(const SMDS_MeshNode* theNode)
{
return false;
}
+//================================================================================
+/*!
+ * \brief Reorient faces.
+ * \param theFaces - the faces to reorient. If empty the whole mesh is meant
+ * \param theDirection - desired direction of normal of \a theFace
+ * \param theFace - one of \a theFaces that sould be oriented according to
+ * \a theDirection and whose orientation defines orientation of other faces
+ * \return number of reoriented faces.
+ */
+//================================================================================
+
+int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet & theFaces,
+ const gp_Dir& theDirection,
+ const SMDS_MeshElement * theFace)
+{
+ int nbReori = 0;
+ if ( !theFace || theFace->GetType() != SMDSAbs_Face ) return nbReori;
+
+ if ( theFaces.empty() )
+ {
+ SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=*/true);
+ while ( fIt->more() )
+ theFaces.insert( theFaces.end(), fIt->next() );
+ }
+
+ // orient theFace according to theDirection
+ gp_XYZ normal;
+ SMESH_Algo::FaceNormal( theFace, normal, /*normalized=*/false );
+ if ( normal * theDirection.XYZ() < 0 )
+ nbReori += Reorient( theFace );
+
+ // Orient other faces
+
+ set< const SMDS_MeshElement* > startFaces, visitedFaces;
+ TIDSortedElemSet avoidSet;
+ set< SMESH_TLink > checkedLinks;
+ pair< set< SMESH_TLink >::iterator, bool > linkIt_isNew;
+
+ if ( theFaces.size() > 1 )// leave 1 face to prevent finding not selected faces
+ theFaces.erase( theFace );
+ startFaces.insert( theFace );
+
+ int nodeInd1, nodeInd2;
+ const SMDS_MeshElement* otherFace;
+ vector< const SMDS_MeshElement* > facesNearLink;
+ vector< std::pair< int, int > > nodeIndsOfFace;
+
+ set< const SMDS_MeshElement* >::iterator startFace = startFaces.begin();
+ while ( !startFaces.empty() )
+ {
+ startFace = startFaces.begin();
+ theFace = *startFace;
+ startFaces.erase( startFace );
+ if ( !visitedFaces.insert( theFace ).second )
+ continue;
+
+ avoidSet.clear();
+ avoidSet.insert(theFace);
+
+ NLink link( theFace->GetNode( 0 ), 0 );
+
+ const int nbNodes = theFace->NbCornerNodes();
+ for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace
+ {
+ link.second = theFace->GetNode(( i+1 ) % nbNodes );
+ linkIt_isNew = checkedLinks.insert( link );
+ if ( !linkIt_isNew.second )
+ {
+ // link has already been checked and won't be encountered more
+ // if the group (theFaces) is manifold
+ //checkedLinks.erase( linkIt_isNew.first );
+ }
+ else
+ {
+ facesNearLink.clear();
+ nodeIndsOfFace.clear();
+ while (( otherFace = FindFaceInSet( link.first, link.second,
+ theFaces, avoidSet, &nodeInd1, &nodeInd2 )))
+ if ( otherFace != theFace)
+ {
+ facesNearLink.push_back( otherFace );
+ nodeIndsOfFace.push_back( make_pair( nodeInd1, nodeInd2 ));
+ avoidSet.insert( otherFace );
+ }
+ if ( facesNearLink.size() > 1 )
+ {
+ // NON-MANIFOLD mesh shell !
+ // select a face most co-directed with theFace,
+ // other faces won't be visited this time
+ gp_XYZ NF, NOF;
+ SMESH_Algo::FaceNormal( theFace, NF, /*normalized=*/false );
+ double proj, maxProj = -1;
+ for ( size_t i = 0; i < facesNearLink.size(); ++i ) {
+ SMESH_Algo::FaceNormal( facesNearLink[i], NOF, /*normalized=*/false );
+ if (( proj = Abs( NF * NOF )) > maxProj ) {
+ maxProj = proj;
+ otherFace = facesNearLink[i];
+ nodeInd1 = nodeIndsOfFace[i].first;
+ nodeInd2 = nodeIndsOfFace[i].second;
+ }
+ }
+ // not to visit rejected faces
+ for ( size_t i = 0; i < facesNearLink.size(); ++i )
+ if ( facesNearLink[i] != otherFace && theFaces.size() > 1 )
+ visitedFaces.insert( facesNearLink[i] );
+ }
+ else if ( facesNearLink.size() == 1 )
+ {
+ otherFace = facesNearLink[0];
+ nodeInd1 = nodeIndsOfFace.back().first;
+ nodeInd2 = nodeIndsOfFace.back().second;
+ }
+ if ( otherFace && otherFace != theFace)
+ {
+ // link must be reverse in otherFace if orientation ot otherFace
+ // is same as that of theFace
+ if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
+ {
+ nbReori += Reorient( otherFace );
+ }
+ startFaces.insert( otherFace );
+ }
+ }
+ std::swap( link.first, link.second ); // reverse the link
+ }
+ }
+ return nbReori;
+}
+
//=======================================================================
//function : getBadRate
//purpose :
if( !elem->IsQuadratic() ) {
// split liner quadrangle
-
+ // for MaxElementLength2D functor we return minimum diagonal for splitting,
+ // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
if ( aBadRate1 <= aBadRate2 ) {
// tr1 + tr2 is better
newElem1 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
-
+ // for MaxElementLength2D functor we return minimum diagonal for splitting,
+ // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
if (aBadRate1 <= aBadRate2) // tr1 + tr2 is better
return 1; // diagonal 1-3
SMESHDS_SubMesh* subMesh = 0;//GetMeshDS()->MeshElements(1);
SMESHDS_SubMesh* fSubMesh = 0;//subMesh;
-
+
SMESH_SequenceOfElemPtr newNodes, newElems;
// map face of volume to it's baricenrtic node
const SMDS_MeshElement* elem = elemIt->next();
if(elem->GetType() == SMDSAbs_0DElement)
continue;
-
+
SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
if ( elem->GetType() == SMDSAbs_Volume )
{
SMDS_FaceIteratorPtr fIt = aMesh->facesIterator();
while ( fIt->more() ) {
const SMDS_MeshElement* face = fIt->next();
- theElems.insert( face );
+ theElems.insert( theElems.end(), face );
}
}
// get all face ids theElems are on
//MESSAGE("sweepElement " << nbSteps);
SMESHDS_Mesh* aMesh = GetMeshDS();
- const int nbNodes = elem->NbNodes();
+ const int nbNodes = elem->NbNodes();
const int nbCorners = elem->NbCornerNodes();
- SMDSAbs_EntityType baseType = elem->GetEntityType(); /* it can change in case of
+ SMDSAbs_EntityType baseType = elem->GetEntityType(); /* it can change in case of
polyhedron creation !!! */
// Loop on elem nodes:
// find new nodes and detect same nodes indices
break;
}
case SMDSEntity_Quad_Quadrangle: { // sweep Quadratic QUADRANGLE --->
- if ( nbDouble != 4 ) break;
if( nbSame == 0 ) {
+ if ( nbDouble != 4 ) break;
// ---> hexahedron with 20 nodes
aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2], prevNod[3],
nextNod[0], nextNod[1], nextNod[2], nextNod[3],
midlNod[0], midlNod[1], midlNod[2], midlNod[3]);
}
else if(nbSame==1) {
- // ---> pyramid + pentahedron - can not be created since it is needed
+ // ---> pyramid + pentahedron - can not be created since it is needed
// additional middle node at the center of face
INFOS( " Sweep for face " << elem->GetID() << " can not be created" );
return;
}
else if( nbSame == 2 ) {
+ if ( nbDouble != 2 ) break;
// ---> 2d order Pentahedron with 15 nodes
int n1,n2,n4,n5;
if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] ) {
}
break;
}
+ case SMDSEntity_Ball:
+ return;
+
default:
break;
}
{
const SMDS_MeshNode* node =
static_cast<const SMDS_MeshNode*>( nList->first );
+ if ( newElemsMap.count( node ))
+ continue; // node was extruded into edge
SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
int nbInitElems = 0;
const SMDS_MeshElement* el = 0;
vecNewNodes[ 1 ]->second.back())) {
myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
vecNewNodes[ 1 ]->second.back()));
- srcElements.Append( myLastCreatedElems.Last() );
+ srcElements.Append( elem );
}
}
else {
myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
vecNewNodes[ 1 ]->second.back(),
vecNewNodes[ 2 ]->second.back()));
- srcElements.Append( myLastCreatedElems.Last() );
+ srcElements.Append( elem );
}
}
}
int iNext = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1;
const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
- // check if a link is free
+ // check if a link n1-n2 is free
if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
hasFreeLinks = true;
- // make an edge and a ceiling for a new edge
- if ( !aMesh->FindEdge( n1, n2 )) {
- myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // free link edge
+ // make a new edge and a ceiling for a new edge
+ const SMDS_MeshElement* edge;
+ if ( ! ( edge = aMesh->FindEdge( n1, n2 ))) {
+ myLastCreatedElems.Append( edge = aMesh->AddEdge( n1, n2 )); // free link edge
srcElements.Append( myLastCreatedElems.Last() );
}
n1 = vecNewNodes[ iNode ]->second.back();
n2 = vecNewNodes[ iNext ]->second.back();
if ( !aMesh->FindEdge( n1, n2 )) {
- myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // ceiling edge
- srcElements.Append( myLastCreatedElems.Last() );
+ myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // new edge ceiling
+ srcElements.Append( edge );
}
}
}
// find medium node
if ( !aMesh->FindEdge( n1, n2, n3 )) {
myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // free link edge
- srcElements.Append( myLastCreatedElems.Last() );
+ srcElements.Append( elem );
}
n1 = vecNewNodes[ iNode ]->second.back();
n2 = vecNewNodes[ iNext ]->second.back();
n3 = vecNewNodes[ iNode+nbn ]->second.back();
if ( !aMesh->FindEdge( n1, n2, n3 )) {
myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // ceiling edge
- srcElements.Append( myLastCreatedElems.Last() );
+ srcElements.Append( elem );
}
}
}
}
while ( srcElements.Length() < myLastCreatedElems.Length() )
- srcElements.Append( myLastCreatedElems.Last() );
+ srcElements.Append( elem );
}
} // loop on swept elements
}
list<SMESH_MeshEditor_PathPoint> fullList;
const TopoDS_Shape& aS = theTrack->GetSubShape();
- // Sub shape for the Pattern must be an Edge or Wire
+ // Sub-shape for the Pattern must be an Edge or Wire
if( aS.ShapeType() == TopAbs_EDGE ) {
aTrackEdge = TopoDS::Edge( aS );
// the Edge must not be degenerated
list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
int startNid = theN1->GetID();
TColStd_MapOfInteger UsedNums;
-
+
int NbEdges = Edges.Length();
int i = 1;
for(; i<=NbEdges; i++) {
const SMDS_MeshElement* currentElem = NULL;
int totalNbEdges = theTrack->NbEdges();
SMDS_ElemIteratorPtr nIt;
- bool isClosed = false;
//check start node
if( !theTrack->GetMeshDS()->Contains(theN1) ) {
return EXTR_BAD_STARTING_NODE;
}
-
+
conn = nbEdgeConnectivity(theN1);
if(conn > 2)
return EXTR_PATH_NOT_EDGE;
if(currentNode == prevNode)
currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
aNodesList.push_back(currentNode);
- } else {
+ } else {
nIt = currentElem->nodesIterator();
while( nIt->more() ) {
currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
if(currentNode == prevNode)
currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
aNodesList.push_back(currentNode);
-
+
//case of the closed mesh
if(currentNode == theN1) {
nbEdges++;
- isClosed = true;
break;
}
conn = nbEdgeConnectivity(currentNode);
if(conn > 2) {
- return EXTR_PATH_NOT_EDGE;
+ return EXTR_PATH_NOT_EDGE;
}else if( conn == 1 && nbEdges > 0 ) {
//End of the path
nbEdges++;
nbEdges++;
}
}
- }
-
+ }
+
if(nbEdges != totalNbEdges)
return EXTR_PATH_NOT_EDGE;
x1 = aNodesList[i-1]->X();x2 = aNodesList[i]->X();
y1 = aNodesList[i-1]->Y();y2 = aNodesList[i]->Y();
z1 = aNodesList[i-1]->Z();z2 = aNodesList[i]->Z();
- TopoDS_Edge e = BRepBuilderAPI_MakeEdge(gp_Pnt(x1,y1,z1),gp_Pnt(x2,y2,z2));
+ TopoDS_Edge e = BRepBuilderAPI_MakeEdge(gp_Pnt(x1,y1,z1),gp_Pnt(x2,y2,z2));
list<SMESH_MeshEditor_PathPoint> LPP;
aPrms.clear();
MakeEdgePathPoints(aPrms, e, (aNodesList[i-1]->GetID()==startNid), LPP);
fullList.pop_back();
}
fullList.push_back(PP1);
-
- } // Sub shape for the Pattern must be an Edge or Wire
+
+ } // Sub-shape for the Pattern must be an Edge or Wire
else if( aS.ShapeType() == TopAbs_EDGE ) {
aTrackEdge = TopoDS::Edge( aS );
// the Edge must not be degenerated
if ( BRep_Tool::Degenerated( aTrackEdge ) )
return EXTR_BAD_PATH_SHAPE;
TopExp::Vertices( aTrackEdge, aV1, aV2 );
- aItN = theTrack->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
- const SMDS_MeshNode* aN1 = aItN->next();
- aItN = theTrack->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
- const SMDS_MeshNode* aN2 = aItN->next();
+ const SMDS_MeshNode* aN1 = 0;
+ const SMDS_MeshNode* aN2 = 0;
+ if ( theTrack->GetSubMesh( aV1 ) && theTrack->GetSubMesh( aV1 )->GetSubMeshDS() ) {
+ aItN = theTrack->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
+ aN1 = aItN->next();
+ }
+ if ( theTrack->GetSubMesh( aV2 ) && theTrack->GetSubMesh( aV2 )->GetSubMeshDS() ) {
+ aItN = theTrack->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
+ aN2 = aItN->next();
+ }
// starting node must be aN1 or aN2
if ( !( aN1 == theN1 || aN2 == theN1 ) )
return EXTR_BAD_STARTING_NODE;
}
}
list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
- int startNid = theN1->GetID();
+ TopoDS_Vertex aVprev;
TColStd_MapOfInteger UsedNums;
int NbEdges = Edges.Length();
int i = 1;
SMESH_subMesh* locTrack = *itLSM;
SMESHDS_SubMesh* locMeshDS = locTrack->GetSubMeshDS();
TopExp::Vertices( aTrackEdge, aV1, aV2 );
- aItN = locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS()->GetNodes();
- const SMDS_MeshNode* aN1 = aItN->next();
- aItN = locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS()->GetNodes();
- const SMDS_MeshNode* aN2 = aItN->next();
- // starting node must be aN1 or aN2
- if ( !( aN1->GetID() == startNid || aN2->GetID() == startNid ) ) continue;
+ bool aN1isOK = false, aN2isOK = false;
+ if ( aVprev.IsNull() ) {
+ // if previous vertex is not yet defined, it means that we in the beginning of wire
+ // and we have to find initial vertex corresponding to starting node theN1
+ const SMDS_MeshNode* aN1 = 0;
+ const SMDS_MeshNode* aN2 = 0;
+
+ if ( locTrack->GetFather()->GetSubMesh(aV1) && locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS() ) {
+ aItN = locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS()->GetNodes();
+ aN1 = aItN->next();
+ }
+ if ( locTrack->GetFather()->GetSubMesh(aV2) && locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS() ) {
+ aItN = locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS()->GetNodes();
+ aN2 = aItN->next();
+ }
+ // starting node must be aN1 or aN2
+ aN1isOK = ( aN1 && aN1 == theN1 );
+ aN2isOK = ( aN2 && aN2 == theN1 );
+ }
+ else {
+ // we have specified ending vertex of the previous edge on the previous iteration
+ // and we have just to check that it corresponds to any vertex in current segment
+ aN1isOK = aVprev.IsSame( aV1 );
+ aN2isOK = aVprev.IsSame( aV2 );
+ }
+ if ( !aN1isOK && !aN2isOK ) continue;
// 2. Collect parameters on the track edge
aPrms.clear();
aItN = locMeshDS->GetNodes();
while ( aItN->more() ) {
- const SMDS_MeshNode* pNode = aItN->next();
+ const SMDS_MeshNode* pNode = aItN->next();
const SMDS_EdgePosition* pEPos =
static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
double aT = pEPos->GetUParameter();
}
list<SMESH_MeshEditor_PathPoint> LPP;
//Extrusion_Error err =
- MakeEdgePathPoints(aPrms, aTrackEdge,(aN1->GetID()==startNid), LPP);
+ MakeEdgePathPoints(aPrms, aTrackEdge, aN1isOK, LPP);
LLPPs.push_back(LPP);
UsedNums.Add(k);
// update startN for search following egde
- if( aN1->GetID() == startNid ) startNid = aN2->GetID();
- else startNid = aN1->GetID();
+ if ( aN1isOK ) aVprev = aV2;
+ else aVprev = aV1;
break;
}
}
SMESH_MeshEditor_PathPoint PP2 = currList.front();
gp_Dir D1 = PP1.Tangent();
gp_Dir D2 = PP2.Tangent();
- gp_Dir Dnew( gp_Vec( (D1.X()+D2.X())/2, (D1.Y()+D2.Y())/2,
- (D1.Z()+D2.Z())/2 ) );
+ gp_Dir Dnew( ( D1.XYZ() + D2.XYZ() ) / 2 );
PP1.SetTangent(Dnew);
fullList.push_back(PP1);
itPP++;
for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
{
const SMDS_MeshElement* elem = *itElem;
- if ( !elem || elem->GetType() == SMDSAbs_Node )
- continue;
+ if ( !elem ) continue;
- int nbNodes = elem->NbNodes();
- int elemType = elem->GetType();
+ SMDSAbs_GeometryType geomType = elem->GetGeomType();
+ int nbNodes = elem->NbNodes();
+ if ( geomType == SMDSGeom_NONE ) continue; // node
- if (elem->IsPoly()) {
+ switch ( geomType ) {
- // polygon or polyhedral volume
- switch ( elemType ) {
- case SMDSAbs_Face:
- {
- vector<const SMDS_MeshNode*> poly_nodes (nbNodes);
- int iNode = 0;
- SMDS_ElemIteratorPtr itN = elem->nodesIterator();
- while (itN->more()) {
- const SMDS_MeshNode* node =
- static_cast<const SMDS_MeshNode*>(itN->next());
- TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
- if (nodeMapIt == nodeMap.end())
- break; // not all nodes transformed
- if (needReverse) {
- // reverse mirrored faces and volumes
- poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second;
- } else {
- poly_nodes[iNode] = (*nodeMapIt).second;
- }
- iNode++;
+ case SMDSGeom_POLYGON: // ---------------------- polygon
+ {
+ vector<const SMDS_MeshNode*> poly_nodes (nbNodes);
+ int iNode = 0;
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while (itN->more()) {
+ const SMDS_MeshNode* node =
+ static_cast<const SMDS_MeshNode*>(itN->next());
+ TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
+ if (nodeMapIt == nodeMap.end())
+ break; // not all nodes transformed
+ if (needReverse) {
+ // reverse mirrored faces and volumes
+ poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second;
+ } else {
+ poly_nodes[iNode] = (*nodeMapIt).second;
}
- if ( iNode != nbNodes )
- continue; // not all nodes transformed
+ iNode++;
+ }
+ if ( iNode != nbNodes )
+ continue; // not all nodes transformed
- if ( theTargetMesh ) {
- myLastCreatedElems.Append(aTgtMesh->AddPolygonalFace(poly_nodes));
- srcElems.Append( elem );
- }
- else if ( theCopy ) {
- myLastCreatedElems.Append(aMesh->AddPolygonalFace(poly_nodes));
- srcElems.Append( elem );
- }
- else {
- aMesh->ChangePolygonNodes(elem, poly_nodes);
- }
+ if ( theTargetMesh ) {
+ myLastCreatedElems.Append(aTgtMesh->AddPolygonalFace(poly_nodes));
+ srcElems.Append( elem );
}
- break;
- case SMDSAbs_Volume:
- {
- const SMDS_VtkVolume* aPolyedre =
- dynamic_cast<const SMDS_VtkVolume*>( elem );
- if (!aPolyedre) {
- MESSAGE("Warning: bad volumic element");
- continue;
- }
+ else if ( theCopy ) {
+ myLastCreatedElems.Append(aMesh->AddPolygonalFace(poly_nodes));
+ srcElems.Append( elem );
+ }
+ else {
+ aMesh->ChangePolygonNodes(elem, poly_nodes);
+ }
+ }
+ break;
- vector<const SMDS_MeshNode*> poly_nodes; poly_nodes.reserve( nbNodes );
- vector<int> quantities;
+ case SMDSGeom_POLYHEDRA: // ------------------ polyhedral volume
+ {
+ const SMDS_VtkVolume* aPolyedre =
+ dynamic_cast<const SMDS_VtkVolume*>( elem );
+ if (!aPolyedre) {
+ MESSAGE("Warning: bad volumic element");
+ continue;
+ }
- bool allTransformed = true;
- int nbFaces = aPolyedre->NbFaces();
- for (int iface = 1; iface <= nbFaces && allTransformed; iface++) {
- int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
- for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) {
- const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode);
- TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
- if (nodeMapIt == nodeMap.end()) {
- allTransformed = false; // not all nodes transformed
- } else {
- poly_nodes.push_back((*nodeMapIt).second);
- }
- if ( needReverse && allTransformed )
- std::reverse( poly_nodes.end() - nbFaceNodes, poly_nodes.end() );
+ vector<const SMDS_MeshNode*> poly_nodes; poly_nodes.reserve( nbNodes );
+ vector<int> quantities; quantities.reserve( nbNodes );
+
+ bool allTransformed = true;
+ int nbFaces = aPolyedre->NbFaces();
+ for (int iface = 1; iface <= nbFaces && allTransformed; iface++) {
+ int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
+ for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) {
+ const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode);
+ TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
+ if (nodeMapIt == nodeMap.end()) {
+ allTransformed = false; // not all nodes transformed
+ } else {
+ poly_nodes.push_back((*nodeMapIt).second);
}
- quantities.push_back(nbFaceNodes);
+ if ( needReverse && allTransformed )
+ std::reverse( poly_nodes.end() - nbFaceNodes, poly_nodes.end() );
}
- if ( !allTransformed )
- continue; // not all nodes transformed
+ quantities.push_back(nbFaceNodes);
+ }
+ if ( !allTransformed )
+ continue; // not all nodes transformed
- if ( theTargetMesh ) {
- myLastCreatedElems.Append(aTgtMesh->AddPolyhedralVolume(poly_nodes, quantities));
- srcElems.Append( elem );
- }
- else if ( theCopy ) {
- myLastCreatedElems.Append(aMesh->AddPolyhedralVolume(poly_nodes, quantities));
- srcElems.Append( elem );
- }
- else {
- aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
- }
+ if ( theTargetMesh ) {
+ myLastCreatedElems.Append(aTgtMesh->AddPolyhedralVolume(poly_nodes, quantities));
+ srcElems.Append( elem );
+ }
+ else if ( theCopy ) {
+ myLastCreatedElems.Append(aMesh->AddPolyhedralVolume(poly_nodes, quantities));
+ srcElems.Append( elem );
+ }
+ else {
+ aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
}
- break;
- default:;
}
- continue;
+ break;
+
+ case SMDSGeom_BALL: // -------------------- Ball
+ {
+ if ( !theCopy && !theTargetMesh ) continue;
- } // elem->isPoly()
+ TNodeNodeMap::iterator nodeMapIt = nodeMap.find( elem->GetNode(0) );
+ if (nodeMapIt == nodeMap.end())
+ continue; // not all nodes transformed
+
+ double diameter = static_cast<const SMDS_BallElement*>(elem)->GetDiameter();
+ if ( theTargetMesh ) {
+ myLastCreatedElems.Append(aTgtMesh->AddBall( nodeMapIt->second, diameter ));
+ srcElems.Append( elem );
+ }
+ else {
+ myLastCreatedElems.Append(aMesh->AddBall( nodeMapIt->second, diameter ));
+ srcElems.Append( elem );
+ }
+ }
+ break;
- // Regular elements
+ default: // ----------------------- Regular elements
- while ( iForw.size() < nbNodes ) iForw.push_back( iForw.size() );
- const std::vector<int>& iRev = SMDS_MeshCell::reverseSmdsOrder( elem->GetEntityType() );
- const std::vector<int>& i = needReverse ? iRev : iForw;
+ while ( iForw.size() < nbNodes ) iForw.push_back( iForw.size() );
+ const std::vector<int>& iRev = SMDS_MeshCell::reverseSmdsOrder( elem->GetEntityType() );
+ const std::vector<int>& i = needReverse ? iRev : iForw;
- // find transformed nodes
- vector<const SMDS_MeshNode*> nodes(nbNodes);
- int iNode = 0;
- SMDS_ElemIteratorPtr itN = elem->nodesIterator();
- while ( itN->more() ) {
- const SMDS_MeshNode* node =
- static_cast<const SMDS_MeshNode*>( itN->next() );
- TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
- if ( nodeMapIt == nodeMap.end() )
- break; // not all nodes transformed
- nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
- }
- if ( iNode != nbNodes )
- continue; // not all nodes transformed
+ // find transformed nodes
+ vector<const SMDS_MeshNode*> nodes(nbNodes);
+ int iNode = 0;
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() ) {
+ const SMDS_MeshNode* node =
+ static_cast<const SMDS_MeshNode*>( itN->next() );
+ TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
+ if ( nodeMapIt == nodeMap.end() )
+ break; // not all nodes transformed
+ nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
+ }
+ if ( iNode != nbNodes )
+ continue; // not all nodes transformed
- if ( theTargetMesh ) {
- if ( SMDS_MeshElement* copy =
- targetMeshEditor.AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
- myLastCreatedElems.Append( copy );
- srcElems.Append( elem );
+ if ( theTargetMesh ) {
+ if ( SMDS_MeshElement* copy =
+ targetMeshEditor.AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
+ myLastCreatedElems.Append( copy );
+ srcElems.Append( elem );
+ }
}
- }
- else if ( theCopy ) {
- if ( AddElement( nodes, elem->GetType(), elem->IsPoly() ))
- srcElems.Append( elem );
- }
- else {
- // reverse element as it was reversed by transformation
- if ( nbNodes > 2 )
- aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes );
- }
+ else if ( theCopy ) {
+ if ( AddElement( nodes, elem->GetType(), elem->IsPoly() ))
+ srcElems.Append( elem );
+ }
+ else {
+ // reverse element as it was reversed by transformation
+ if ( nbNodes > 2 )
+ aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes );
+ }
+ } // switch ( geomType )
} // loop on elements
// Sort existing groups by types and collect their names
- // to store an old group and a generated new one
- typedef pair< SMESHDS_GroupBase*, SMDS_MeshGroup* > TOldNewGroup;
+ // to store an old group and a generated new ones
+ using boost::tuple;
+ using boost::make_tuple;
+ typedef tuple< SMESHDS_GroupBase*, SMESHDS_Group*, SMESHDS_Group* > TOldNewGroup;
vector< list< TOldNewGroup > > groupsByType( SMDSAbs_NbElementTypes );
+ vector< TOldNewGroup* > orderedOldNewGroups; // in order of old groups
// group names
set< string > groupNames;
- //
- SMDS_MeshGroup* nullNewGroup = (SMDS_MeshGroup*) 0;
+
SMESH_Mesh::GroupIteratorPtr groupIt = GetMesh()->GetGroups();
- while ( groupIt->more() ) {
+ if ( !groupIt->more() ) return newGroupIDs;
+
+ int newGroupID = mesh->GetGroupIds().back()+1;
+ while ( groupIt->more() )
+ {
SMESH_Group * group = groupIt->next();
if ( !group ) continue;
SMESHDS_GroupBase* groupDS = group->GetGroupDS();
if ( !groupDS || groupDS->IsEmpty() ) continue;
- groupNames.insert( group->GetName() );
+ groupNames.insert ( group->GetName() );
groupDS->SetStoreName( group->GetName() );
- groupsByType[ groupDS->GetType() ].push_back( make_pair( groupDS, nullNewGroup ));
+ const SMDSAbs_ElementType type = groupDS->GetType();
+ SMESHDS_Group* newGroup = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(), type );
+ SMESHDS_Group* newTopGroup = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(), type );
+ groupsByType[ groupDS->GetType() ].push_back( make_tuple( groupDS, newGroup, newTopGroup ));
+ orderedOldNewGroups.push_back( & groupsByType[ groupDS->GetType() ].back() );
}
- // Groups creation
+ // Loop on nodes and elements to add them in new groups
- // loop on nodes and elements
for ( int isNodes = 0; isNodes < 2; ++isNodes )
{
const SMESH_SequenceOfElemPtr& gens = isNodes ? nodeGens : elemGens;
const SMESH_SequenceOfElemPtr& elems = isNodes ? myLastCreatedNodes : myLastCreatedElems;
if ( gens.Length() != elems.Length() )
- throw SALOME_Exception(LOCALIZED("invalid args"));
+ throw SALOME_Exception("SMESH_MeshEditor::generateGroups(): invalid args");
// loop on created elements
for (int iElem = 1; iElem <= elems.Length(); ++iElem )
continue;
}
list< TOldNewGroup > & groupsOldNew = groupsByType[ sourceElem->GetType() ];
- if ( groupsOldNew.empty() ) {
+ if ( groupsOldNew.empty() ) { // no groups of this type at all
while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
++iElem; // skip all elements made by sourceElem
continue;
}
- // collect all elements made by sourceElem
+ // collect all elements made by the iElem-th sourceElem
list< const SMDS_MeshElement* > resultElems;
if ( const SMDS_MeshElement* resElem = elems( iElem ))
if ( resElem != sourceElem )
if ( const SMDS_MeshElement* resElem = elems( ++iElem ))
if ( resElem != sourceElem )
resultElems.push_back( resElem );
- // do not generate element groups from node ones
- if ( sourceElem->GetType() == SMDSAbs_Node &&
- elems( iElem )->GetType() != SMDSAbs_Node )
- continue;
// add resultElems to groups made by ones the sourceElem belongs to
list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end();
for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew )
{
- SMESHDS_GroupBase* oldGroup = gOldNew->first;
- if ( oldGroup->Contains( sourceElem )) // sourceElem in oldGroup
+ SMESHDS_GroupBase* oldGroup = gOldNew->get<0>();
+ if ( oldGroup->Contains( sourceElem )) // sourceElem is in oldGroup
{
- SMDS_MeshGroup* & newGroup = gOldNew->second;
- if ( !newGroup )// create a new group
- {
- // make a name
- string name = oldGroup->GetStoreName();
- if ( !targetMesh ) {
- name += "_";
- name += postfix;
- int nb = 0;
- while ( !groupNames.insert( name ).second ) // name exists
- {
- if ( nb == 0 ) {
- name += "_1";
- }
- else {
- TCollection_AsciiString nbStr(nb+1);
- name.resize( name.rfind('_')+1 );
- name += nbStr.ToCString();
- }
- ++nb;
- }
- }
- // make a group
- int id;
- SMESH_Group* group = mesh->AddGroup( resultElems.back()->GetType(),
- name.c_str(), id );
- SMESHDS_Group* groupDS = static_cast<SMESHDS_Group*>(group->GetGroupDS());
- newGroup = & groupDS->SMDSGroup();
- newGroupIDs->push_back( id );
- }
-
// fill in a new group
+ SMDS_MeshGroup & newGroup = gOldNew->get<1>()->SMDSGroup();
+ list< const SMDS_MeshElement* > rejectedElems; // elements of other type
list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt )
- newGroup->Add( *resElemIt );
+ if ( !newGroup.Add( *resElemIt ))
+ rejectedElems.push_back( *resElemIt );
+
+ // fill "top" group
+ if ( !rejectedElems.empty() )
+ {
+ SMDS_MeshGroup & newTopGroup = gOldNew->get<2>()->SMDSGroup();
+ resLast = rejectedElems.end();
+ for ( resElemIt = rejectedElems.begin(); resElemIt != resLast; ++resElemIt )
+ !newTopGroup.Add( *resElemIt );
+ }
}
}
} // loop on created elements
}// loop on nodes and elements
+ // Create new SMESH_Groups from SMESHDS_Groups and remove empty SMESHDS_Groups
+
+ list<int> topGrouIds;
+ for ( size_t i = 0; i < orderedOldNewGroups.size(); ++i )
+ {
+ SMESHDS_GroupBase* oldGroupDS = orderedOldNewGroups[i]->get<0>();
+ SMESHDS_Group* newGroups[2] = { orderedOldNewGroups[i]->get<1>(),
+ orderedOldNewGroups[i]->get<2>() };
+ const int nbNewGroups = !newGroups[0]->IsEmpty() + !newGroups[1]->IsEmpty();
+ for ( int is2nd = 0; is2nd < 2; ++is2nd )
+ {
+ SMESHDS_Group* newGroupDS = newGroups[ is2nd ];
+ if ( newGroupDS->IsEmpty() )
+ {
+ mesh->GetMeshDS()->RemoveGroup( newGroupDS );
+ }
+ else
+ {
+ // set group type
+ newGroupDS->SetType( newGroupDS->GetElements()->next()->GetType() );
+
+ // make a name
+ const bool isTop = ( nbNewGroups == 2 &&
+ newGroupDS->GetType() == oldGroupDS->GetType() );
+ string name = oldGroupDS->GetStoreName();
+ if ( !targetMesh ) {
+ string suffix = ( isTop ? "top": postfix.c_str() );
+ name += "_";
+ name += suffix;
+ int nb = 1;
+ while ( !groupNames.insert( name ).second ) // name exists
+ name = SMESH_Comment( oldGroupDS->GetStoreName() ) << "_" << suffix << "_" << nb++;
+ }
+ else if ( isTop ) {
+ name += "_top";
+ }
+ newGroupDS->SetStoreName( name.c_str() );
+
+ // make a SMESH_Groups
+ mesh->AddGroup( newGroupDS );
+ if ( isTop )
+ topGrouIds.push_back( newGroupDS->GetID() );
+ else
+ newGroupIDs->push_back( newGroupDS->GetID() );
+ }
+ }
+ }
+ newGroupIDs->splice( newGroupIDs->end(), topGrouIds );
+
return newGroupIDs;
}
}
else if ( tree->NbNodes() ) // put a tree to the treeMap
{
- const Bnd_B3d& box = tree->getBox();
+ const Bnd_B3d& box = *tree->getBox();
double sqDist = thePnt.SquareDistance( 0.5 * ( box.CornerMin() + box.CornerMax() ));
pair<TDistTreeMap::iterator,bool> it_in = treeMap.insert( make_pair( sqDist, tree ));
if ( !it_in.second ) // not unique distance to box center
TDistTreeMap::iterator sqDist_tree = treeMap.begin();
if ( treeMap.size() > 5 ) {
SMESH_OctreeNode* closestTree = sqDist_tree->second;
- const Bnd_B3d& box = closestTree->getBox();
+ const Bnd_B3d& box = *closestTree->getBox();
double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() );
sqLimit = limit * limit;
}
*/
//=======================================================================
-SMESH_NodeSearcher* SMESH_MeshEditor::GetNodeSearcher()
+SMESH_NodeSearcher* SMESH_MeshEditor::GetNodeSearcher()
{
return new SMESH_NodeSearcherImpl( GetMeshDS() );
}
{
public:
- ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(), double tolerance = NodeRadius );
- void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems);
+ ElementBndBoxTree(const SMDS_Mesh& mesh,
+ SMDSAbs_ElementType elemType,
+ SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(),
+ double tolerance = NodeRadius );
+ void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems );
void getElementsNearLine ( const gp_Ax1& line, TIDSortedElemSet& foundElems);
+ void getElementsInSphere ( const gp_XYZ& center,
+ const double radius, TIDSortedElemSet& foundElems);
+ size_t getSize() { return std::max( _size, _elements.size() ); }
~ElementBndBoxTree();
protected:
- ElementBndBoxTree() {}
- SMESH_Octree* allocateOctreeChild() const { return new ElementBndBoxTree; }
- void buildChildrenData();
- Bnd_B3d* buildRootBox();
+ ElementBndBoxTree():_size(0) {}
+ SMESH_Octree* newChild() const { return new ElementBndBoxTree; }
+ void buildChildrenData();
+ Bnd_B3d* buildRootBox();
private:
//!< Bounding box of element
struct ElementBox : public Bnd_B3d
ElementBox(const SMDS_MeshElement* elem, double tolerance);
};
vector< ElementBox* > _elements;
+ size_t _size;
};
//================================================================================
//================================================================================
ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt, double tolerance)
- :SMESH_Octree( new SMESH_Octree::Limit( MaxLevel, /*minSize=*/0. ))
+ :SMESH_Octree( new SMESH_TreeLimit( MaxLevel, /*minSize=*/0. ))
{
int nbElems = mesh.GetMeshInfo().NbElements( elemType );
_elements.reserve( nbElems );
while ( elemIt->more() )
_elements.push_back( new ElementBox( elemIt->next(),tolerance ));
- if ( _elements.size() > MaxNbElemsInLeaf )
- compute();
- else
- myIsLeaf = true;
+ compute();
}
//================================================================================
{
for (int j = 0; j < 8; j++)
{
- if ( !_elements[i]->IsOut( myChildren[j]->getBox() ))
+ if ( !_elements[i]->IsOut( *myChildren[j]->getBox() ))
{
_elements[i]->_refCount++;
((ElementBndBoxTree*)myChildren[j])->_elements.push_back( _elements[i]);
}
_elements[i]->_refCount--;
}
- _elements.clear();
+ _size = _elements.size();
+ SMESHUtils::FreeVector( _elements ); // = _elements.clear() + free memory
for (int j = 0; j < 8; j++)
{
child->myIsLeaf = true;
if ( child->_elements.capacity() - child->_elements.size() > 1000 )
- child->_elements.resize( child->_elements.size() ); // compact
+ SMESHUtils::CompactVector( child->_elements );
}
}
void ElementBndBoxTree::getElementsNearPoint( const gp_Pnt& point,
TIDSortedElemSet& foundElems)
{
- if ( level() && getBox().IsOut( point.XYZ() ))
+ if ( getBox()->IsOut( point.XYZ() ))
return;
if ( isLeaf() )
void ElementBndBoxTree::getElementsNearLine( const gp_Ax1& line,
TIDSortedElemSet& foundElems)
{
- if ( level() && getBox().IsOut( line ))
+ if ( getBox()->IsOut( line ))
return;
if ( isLeaf() )
}
}
+ //================================================================================
+ /*!
+ * \brief Return elements from leaves intersecting the sphere
+ */
+ //================================================================================
+
+ void ElementBndBoxTree::getElementsInSphere ( const gp_XYZ& center,
+ const double radius,
+ TIDSortedElemSet& foundElems)
+ {
+ if ( getBox()->IsOut( center, radius ))
+ return;
+
+ if ( isLeaf() )
+ {
+ for ( int i = 0; i < _elements.size(); ++i )
+ if ( !_elements[i]->IsOut( center, radius ))
+ foundElems.insert( _elements[i]->_element );
+ }
+ else
+ {
+ for (int i = 0; i < 8; i++)
+ ((ElementBndBoxTree*) myChildren[i])->getElementsInSphere( center, radius, foundElems );
+ }
+ }
+
//================================================================================
/*!
* \brief Construct the element box
_refCount = 1;
SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
while ( nIt->more() )
- Add( SMESH_TNodeXYZ( cast2Node( nIt->next() )));
+ Add( SMESH_TNodeXYZ( nIt->next() ));
Enlarge( tolerance );
}
SMDSAbs_ElementType type,
vector< const SMDS_MeshElement* >& foundElements);
virtual TopAbs_State GetPointState(const gp_Pnt& point);
+ virtual const SMDS_MeshElement* FindClosestTo( const gp_Pnt& point,
+ SMDSAbs_ElementType type );
void GetElementsNearLine( const gp_Ax1& line,
SMDSAbs_ElementType type,
GeomAPI_ExtremaCurveCurve anExtCC;
Handle(Geom_Curve) lineCurve = new Geom_Line( line );
-
+
int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes();
for ( int i = 0; i < nbNodes && nbInts < 2; ++i )
{
GC_MakeSegment edge( SMESH_TNodeXYZ( face->GetNode( i )),
- SMESH_TNodeXYZ( face->GetNode( (i+1)%nbNodes) ));
+ SMESH_TNodeXYZ( face->GetNode( (i+1)%nbNodes) ));
anExtCC.Init( lineCurve, edge);
if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol)
{
while (( f = SMESH_MeshEditor::FindFaceInSet(link.node1(), link.node2(), emptySet, faces )))
faces.insert( f );
- // select another outer face among the found
+ // select another outer face among the found
const SMDS_MeshElement* outerFace2 = 0;
if ( faces.size() == 2 )
{
// 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
{
* \brief Find elements of given type where the given point is IN or ON.
* Returns nb of found elements and elements them-selves.
*
- * 'ALL' type means elements of any type excluding nodes and 0D elements
+ * 'ALL' type means elements of any type excluding nodes, balls and 0D elements
*/
//=======================================================================
double tolerance = getTolerance();
// =================================================================================
- if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement )
+ if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement || type == SMDSAbs_Ball)
{
if ( !_nodeSearcher )
_nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
}
else
{
- SMDS_ElemIteratorPtr elemIt = closeNode->GetInverseElementIterator( SMDSAbs_0DElement );
+ SMDS_ElemIteratorPtr elemIt = closeNode->GetInverseElementIterator( type );
while ( elemIt->more() )
foundElements.push_back( elemIt->next() );
}
_ebbTree->getElementsNearPoint( point, suspectElems );
TIDSortedElemSet::iterator elem = suspectElems.begin();
for ( ; elem != suspectElems.end(); ++elem )
- if ( !SMESH_MeshEditor::isOut( *elem, point, tolerance ))
+ if ( !SMESH_MeshEditor::IsOut( *elem, point, tolerance ))
foundElements.push_back( *elem );
}
return foundElements.size();
}
-//================================================================================
+//=======================================================================
/*!
- * \brief Classify the given point in the closed 2D mesh
+ * \brief Find an element of given type most close to the given point
+ *
+ * WARNING: Only face search is implemeneted so far
*/
-//================================================================================
+//=======================================================================
-TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
+const SMDS_MeshElement*
+SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt& point,
+ SMDSAbs_ElementType type )
{
- double tolerance = getTolerance();
- if ( !_ebbTree || _elementType != SMDSAbs_Face )
- {
- if ( _ebbTree ) delete _ebbTree;
- _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = SMDSAbs_Face, _meshPartIt );
- }
- // 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 SMDS_MeshElement* closestElem = 0;
- const int nbAxes = 3;
+ if ( type == SMDSAbs_Face )
+ {
+ if ( !_ebbTree || _elementType != type )
+ {
+ if ( _ebbTree ) delete _ebbTree;
+ _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
+ }
+ TIDSortedElemSet suspectElems;
+ _ebbTree->getElementsNearPoint( point, suspectElems );
+
+ if ( suspectElems.empty() && _ebbTree->maxSize() > 0 )
+ {
+ gp_Pnt boxCenter = 0.5 * ( _ebbTree->getBox()->CornerMin() +
+ _ebbTree->getBox()->CornerMax() );
+ double radius;
+ if ( _ebbTree->getBox()->IsOut( point.XYZ() ))
+ radius = point.Distance( boxCenter ) - 0.5 * _ebbTree->maxSize();
+ else
+ radius = _ebbTree->maxSize() / pow( 2., _ebbTree->getHeight()) / 2;
+ while ( suspectElems.empty() )
+ {
+ _ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems );
+ radius *= 1.1;
+ }
+ }
+ double minDist = std::numeric_limits<double>::max();
+ multimap< double, const SMDS_MeshElement* > dist2face;
+ TIDSortedElemSet::iterator elem = suspectElems.begin();
+ for ( ; elem != suspectElems.end(); ++elem )
+ {
+ double dist = SMESH_MeshEditor::GetDistance( dynamic_cast<const SMDS_MeshFace*>(*elem),
+ point );
+ if ( dist < minDist + 1e-10)
+ {
+ minDist = dist;
+ dist2face.insert( dist2face.begin(), make_pair( dist, *elem ));
+ }
+ }
+ if ( !dist2face.empty() )
+ {
+ multimap< double, const SMDS_MeshElement* >::iterator d2f = dist2face.begin();
+ closestElem = d2f->second;
+ // if there are several elements at the same distance, select one
+ // with GC closest to the point
+ typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
+ double minDistToGC = 0;
+ for ( ++d2f; d2f != dist2face.end() && fabs( d2f->first - minDist ) < 1e-10; ++d2f )
+ {
+ if ( minDistToGC == 0 )
+ {
+ gp_XYZ gc(0,0,0);
+ gc = accumulate( TXyzIterator(closestElem->nodesIterator()),
+ TXyzIterator(), gc ) / closestElem->NbNodes();
+ minDistToGC = point.SquareDistance( gc );
+ }
+ gp_XYZ gc(0,0,0);
+ gc = accumulate( TXyzIterator( d2f->second->nodesIterator()),
+ TXyzIterator(), gc ) / d2f->second->NbNodes();
+ double d = point.SquareDistance( gc );
+ if ( d < minDistToGC )
+ {
+ minDistToGC = d;
+ closestElem = d2f->second;
+ }
+ }
+ // cout << "FindClosestTo( " <<point.X()<<", "<<point.Y()<<", "<<point.Z()<<" ) FACE "
+ // <<closestElem->GetID() << " DIST " << minDist << endl;
+ }
+ }
+ else
+ {
+ // NOT IMPLEMENTED SO FAR
+ }
+ return closestElem;
+}
+
+
+//================================================================================
+/*!
+ * \brief Classify the given point in the closed 2D mesh
+ */
+//================================================================================
+
+TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
+{
+ double tolerance = getTolerance();
+ if ( !_ebbTree || _elementType != SMDSAbs_Face )
+ {
+ if ( _ebbTree ) delete _ebbTree;
+ _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = SMDSAbs_Face, _meshPartIt );
+ }
+ // 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() };
map< double, TInters > paramOnLine2TInters[ nbAxes ];
list< TInters > tangentInters[ nbAxes ]; // of faces whose plane includes the line
else if ( ! intersection.IsParallel() && intersection.NbPoints() > 0 )
{
gp_Pnt intersectionPoint = intersection.Point(1);
- if ( !SMESH_MeshEditor::isOut( *face, intersectionPoint, tolerance ))
+ if ( !SMESH_MeshEditor::IsOut( *face, intersectionPoint, tolerance ))
u2inters.insert(make_pair( intersection.ParamOnConic(1), TInters( *face, fNorm )));
}
}
int nbInter = u2inters.size();
if ( nbInter == 0 )
- return TopAbs_OUT;
+ return TopAbs_OUT;
double f = u2inters.begin()->first, l = u2inters.rbegin()->first;
if ( nbInter == 1 ) // not closed mesh
return TopAbs_ON;
if ( nbIntBeforePoint == 0 || nbIntAfterPoint == 0)
- return TopAbs_OUT;
+ return TopAbs_OUT;
if ( nbIntBeforePoint + nbIntAfterPoint == 1 ) // not closed mesh
return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN;
*/
//=======================================================================
-bool SMESH_MeshEditor::isOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol )
+bool SMESH_MeshEditor::IsOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol )
{
if ( element->GetType() == SMDSAbs_Volume)
{
for ( i = 0; i < nbNodes; ++i )
{
SMDS_LinearEdge edge( nodeList[i], nodeList[i+1] );
- if ( !isOut( &edge, point, tol ))
+ if ( !IsOut( &edge, point, tol ))
return false;
}
return true;
return true;
}
+//=======================================================================
+
+namespace
+{
+ // Position of a point relative to a segment
+ // . .
+ // . LEFT .
+ // . .
+ // VERTEX 1 o----ON-----> VERTEX 2
+ // . .
+ // . RIGHT .
+ // . .
+ enum PositionName { POS_LEFT = 1, POS_VERTEX = 2, POS_RIGHT = 4, //POS_ON = 8,
+ POS_ALL = POS_LEFT | POS_RIGHT | POS_VERTEX };
+ struct PointPos
+ {
+ PositionName _name;
+ int _index; // index of vertex or segment
+
+ PointPos( PositionName n, int i=-1 ): _name(n), _index(i) {}
+ bool operator < (const PointPos& other ) const
+ {
+ if ( _name == other._name )
+ return ( _index < 0 || other._index < 0 ) ? false : _index < other._index;
+ return _name < other._name;
+ }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Return of a point relative to a segment
+ * \param point2D - the point to analyze position of
+ * \param xyVec - end points of segments
+ * \param index0 - 0-based index of the first point of segment
+ * \param posToFindOut - flags of positions to detect
+ * \retval PointPos - point position
+ */
+ //================================================================================
+
+ PointPos getPointPosition( const gp_XY& point2D,
+ const gp_XY* segEnds,
+ const int index0 = 0,
+ const int posToFindOut = POS_ALL)
+ {
+ const gp_XY& p1 = segEnds[ index0 ];
+ const gp_XY& p2 = segEnds[ index0+1 ];
+ const gp_XY grad = p2 - p1;
+
+ if ( posToFindOut & POS_VERTEX )
+ {
+ // check if the point2D is at "vertex 1" zone
+ gp_XY pp1[2] = { p1, gp_XY( p1.X() - grad.Y(),
+ p1.Y() + grad.X() ) };
+ if ( getPointPosition( point2D, pp1, 0, POS_LEFT|POS_RIGHT )._name == POS_LEFT )
+ return PointPos( POS_VERTEX, index0 );
+
+ // check if the point2D is at "vertex 2" zone
+ gp_XY pp2[2] = { p2, gp_XY( p2.X() - grad.Y(),
+ p2.Y() + grad.X() ) };
+ if ( getPointPosition( point2D, pp2, 0, POS_LEFT|POS_RIGHT )._name == POS_RIGHT )
+ return PointPos( POS_VERTEX, index0 + 1);
+ }
+ double edgeEquation =
+ ( point2D.X() - p1.X() ) * grad.Y() - ( point2D.Y() - p1.Y() ) * grad.X();
+ return PointPos( edgeEquation < 0 ? POS_LEFT : POS_RIGHT, index0 );
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Return minimal distance from a point to a face
+ *
+ * Currently we ignore non-planarity and 2nd order of face
+ */
+//=======================================================================
+
+double SMESH_MeshEditor::GetDistance( const SMDS_MeshFace* face,
+ const gp_Pnt& point )
+{
+ double badDistance = -1;
+ if ( !face ) return badDistance;
+
+ // coordinates of nodes (medium nodes, if any, ignored)
+ typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
+ vector<gp_XYZ> xyz( TXyzIterator( face->nodesIterator()), TXyzIterator() );
+ xyz.resize( face->NbCornerNodes()+1 );
+
+ // transformation to get xyz[0] lies on the origin, xyz[1] lies on the Z axis,
+ // and xyz[2] lies in the XZ plane. This is to pass to 2D space on XZ plane.
+ gp_Trsf trsf;
+ gp_Vec OZ ( xyz[0], xyz[1] );
+ gp_Vec OX ( xyz[0], xyz[2] );
+ if ( OZ.Magnitude() < std::numeric_limits<double>::min() )
+ {
+ if ( xyz.size() < 4 ) return badDistance;
+ OZ = gp_Vec ( xyz[0], xyz[2] );
+ OX = gp_Vec ( xyz[0], xyz[3] );
+ }
+ gp_Ax3 tgtCS;
+ try {
+ tgtCS = gp_Ax3( xyz[0], OZ, OX );
+ }
+ catch ( Standard_Failure ) {
+ return badDistance;
+ }
+ trsf.SetTransformation( tgtCS );
+
+ // move all the nodes to 2D
+ vector<gp_XY> xy( xyz.size() );
+ for ( size_t i = 0;i < xyz.size()-1; ++i )
+ {
+ gp_XYZ p3d = xyz[i];
+ trsf.Transforms( p3d );
+ xy[i].SetCoord( p3d.X(), p3d.Z() );
+ }
+ xyz.back() = xyz.front();
+ xy.back() = xy.front();
+
+ // // move the point in 2D
+ gp_XYZ tmpPnt = point.XYZ();
+ trsf.Transforms( tmpPnt );
+ gp_XY point2D( tmpPnt.X(), tmpPnt.Z() );
+
+ // loop on segments of the face to analyze point position ralative to the face
+ set< PointPos > pntPosSet;
+ for ( size_t i = 1; i < xy.size(); ++i )
+ {
+ PointPos pos = getPointPosition( point2D, &xy[0], i-1 );
+ pntPosSet.insert( pos );
+ }
+
+ // compute distance
+ PointPos pos = *pntPosSet.begin();
+ // cout << "Face " << face->GetID() << " DIST: ";
+ switch ( pos._name )
+ {
+ case POS_LEFT: {
+ // point is most close to a segment
+ gp_Vec p0p1( point, xyz[ pos._index ] );
+ gp_Vec p1p2( xyz[ pos._index ], xyz[ pos._index+1 ]); // segment vector
+ p1p2.Normalize();
+ double projDist = p0p1 * p1p2; // distance projected to the segment
+ gp_Vec projVec = p1p2 * projDist;
+ gp_Vec distVec = p0p1 - projVec;
+ // cout << distVec.Magnitude() << ", SEG " << face->GetNode(pos._index)->GetID()
+ // << " - " << face->GetNodeWrap(pos._index+1)->GetID() << endl;
+ return distVec.Magnitude();
+ }
+ case POS_RIGHT: {
+ // point is inside the face
+ double distToFacePlane = tmpPnt.Y();
+ // cout << distToFacePlane << ", INSIDE " << endl;
+ return Abs( distToFacePlane );
+ }
+ case POS_VERTEX: {
+ // point is most close to a node
+ gp_Vec distVec( point, xyz[ pos._index ]);
+ // cout << distVec.Magnitude() << " VERTEX " << face->GetNode(pos._index)->GetID() << endl;
+ return distVec.Magnitude();
+ }
+ }
+ return badDistance;
+}
+
//=======================================================================
//function : SimplifyFace
//purpose :
}
}
// BUG 0020185: end
- iRepl[ nbRepl++ ] = iCur;
}
curNodes[ iCur ] = n;
bool isUnique = nodeSet.insert( n ).second;
- if ( isUnique ) {
+ if ( isUnique )
uniqueNodes[ iUnique++ ] = n;
- if ( nbRepl && iRepl[ nbRepl-1 ] == iCur )
- --nbRepl; // n do not stick to a node of the elem
- }
+ else
+ iRepl[ nbRepl++ ] = iCur;
iCur++;
}
//purpose : Return list of group of elements built on the same nodes.
// Search among theElements or in the whole mesh if theElements is empty
//=======================================================================
-void SMESH_MeshEditor::FindEqualElements(set<const SMDS_MeshElement*> & theElements,
- TListOfListOfElementsID & theGroupsOfElementsID)
+
+void SMESH_MeshEditor::FindEqualElements(TIDSortedElemSet & theElements,
+ TListOfListOfElementsID & theGroupsOfElementsID)
{
myLastCreatedElems.Clear();
myLastCreatedNodes.Clear();
- typedef set<const SMDS_MeshElement*> TElemsSet;
typedef map< SortableElement, int > TMapOfNodeSet;
typedef list<int> TGroupOfElems;
- TElemsSet elems;
if ( theElements.empty() )
{ // get all elements in the mesh
SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
while ( eIt->more() )
- elems.insert( elems.end(), eIt->next());
+ theElements.insert( theElements.end(), eIt->next());
}
- else
- elems = theElements;
vector< TGroupOfElems > arrayOfGroups;
TGroupOfElems groupOfElems;
TMapOfNodeSet mapOfNodeSet;
- TElemsSet::iterator elemIt = elems.begin();
- for ( int i = 0, j=0; elemIt != elems.end(); ++elemIt, ++j ) {
+ TIDSortedElemSet::iterator elemIt = theElements.begin();
+ for ( int i = 0, j=0; elemIt != theElements.end(); ++elemIt, ++j ) {
const SMDS_MeshElement* curElem = *elemIt;
SortableElement SE(curElem);
int ind = -1;
void SMESH_MeshEditor::MergeEqualElements()
{
- set<const SMDS_MeshElement*> aMeshElements; /* empty input -
- to merge equal elements in the whole mesh */
+ TIDSortedElemSet aMeshElements; /* empty input ==
+ to merge equal elements in the whole mesh */
TListOfListOfElementsID aGroupsOfElementsID;
FindEqualElements(aMeshElements, aGroupsOfElementsID);
MergeElements(aGroupsOfElementsID);
{
nbElem++;
const SMDS_MeshElement* elem = ElemItr->next();
- if( !elem || elem->IsQuadratic() ) continue;
+ if( !elem ) continue;
+ const SMDSAbs_EntityType aGeomType = elem->GetEntityType();
+ if ( elem->IsQuadratic() )
+ {
+ bool alreadyOK;
+ switch ( aGeomType ) {
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa: alreadyOK = !theHelper.GetIsBiQuadratic(); break;
+ case SMDSEntity_BiQuad_Quadrangle:
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theHelper.GetIsBiQuadratic(); break;
+ default: alreadyOK = true;
+ }
+ if ( alreadyOK ) continue;
+ }
// get elem data needed to re-create it
//
- int id = elem->GetID();
- int nbNodes = elem->NbNodes();
- SMDSAbs_ElementType aType = elem->GetType();
+ const int id = elem->GetID();
+ const int nbNodes = elem->NbCornerNodes();
+ const SMDSAbs_ElementType aType = elem->GetType();
nodes.assign(elem->begin_nodes(), elem->end_nodes());
- if ( elem->GetEntityType() == SMDSEntity_Polyhedra )
+ if ( aGeomType == SMDSEntity_Polyhedra )
nbNodeInFaces = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
- else if ( elem->GetEntityType() == SMDSEntity_Hexagonal_Prism )
+ else if ( aGeomType == SMDSEntity_Hexagonal_Prism )
volumeToPolyhedron( elem, nodes, nbNodeInFaces );
// remove a linear element
}
case SMDSAbs_Volume :
{
- switch( elem->GetEntityType() )
+ switch( aGeomType )
{
case SMDSEntity_Tetra:
NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], id, theForce3d);
break;
case SMDSEntity_Hexa:
+ case SMDSEntity_Quad_Hexa:
+ case SMDSEntity_TriQuad_Hexa:
NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
break;
if( NewElem )
theSm->AddElement( NewElem );
}
-// if (!GetMeshDS()->isCompacted())
-// GetMeshDS()->compactMesh();
return nbElem;
}
-
//=======================================================================
//function : ConvertToQuadratic
//purpose :
//=======================================================================
-void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
+void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theToBiQuad)
{
SMESHDS_Mesh* meshDS = GetMeshDS();
SMESH_MesherHelper aHelper(*myMesh);
+
aHelper.SetIsQuadratic( true );
+ aHelper.SetIsBiQuadratic( theToBiQuad );
+ aHelper.SetElementsOnShape(true);
int nbCheckedElems = 0;
if ( myMesh->HasShapeToMesh() )
while(aFaceItr->more())
{
const SMDS_MeshFace* face = aFaceItr->next();
- if(!face || face->IsQuadratic() ) continue;
+ if ( !face ) continue;
+
+ const SMDSAbs_EntityType type = face->GetEntityType();
+ if (( theToBiQuad && type == SMDSEntity_BiQuad_Quadrangle ) ||
+ ( !theToBiQuad && type == SMDSEntity_Quad_Quadrangle ))
+ continue;
const int id = face->GetID();
- const SMDSAbs_EntityType type = face->GetEntityType();
vector<const SMDS_MeshNode *> nodes ( face->begin_nodes(), face->end_nodes());
meshDS->RemoveFreeElement(face, smDS, /*fromGroups=*/false);
const SMDS_MeshVolume* volume = aVolumeItr->next();
if(!volume || volume->IsQuadratic() ) continue;
- const int id = volume->GetID();
const SMDSAbs_EntityType type = volume->GetEntityType();
+ if (( theToBiQuad && type == SMDSEntity_TriQuad_Hexa ) ||
+ ( !theToBiQuad && type == SMDSEntity_Quad_Hexa ))
+ continue;
+
+ const int id = volume->GetID();
vector<const SMDS_MeshNode *> nodes (volume->begin_nodes(), volume->end_nodes());
if ( type == SMDSEntity_Polyhedra )
nbNodeInFaces = static_cast<const SMDS_VtkVolume* >(volume)->GetQuantities();
switch ( type )
{
case SMDSEntity_Tetra:
- NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
- nodes[3], id, theForce3d );
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d );
break;
case SMDSEntity_Hexa:
+ case SMDSEntity_Quad_Hexa:
+ case SMDSEntity_TriQuad_Hexa:
NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
break;
if ( !theForce3d )
{ // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
aHelper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
- aHelper.FixQuadraticElements();
+ aHelper.FixQuadraticElements(myError);
}
}
/*!
* \brief Makes given elements quadratic
* \param theForce3d - if true, the medium nodes will be placed in the middle of link
- * \param theElements - elements to make quadratic
+ * \param theElements - elements to make quadratic
*/
//================================================================================
void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d,
- TIDSortedElemSet& theElements)
+ TIDSortedElemSet& theElements,
+ const bool theToBiQuad)
{
if ( theElements.empty() ) return;
// we believe that all theElements are of the same type
- SMDSAbs_ElementType elemType = (*theElements.begin())->GetType();
-
+ const SMDSAbs_ElementType elemType = (*theElements.begin())->GetType();
+
// get all nodes shared by theElements
TIDSortedNodeSet allNodes;
TIDSortedElemSet::iterator eIt = theElements.begin();
const SMDS_MeshElement* e = invIt->next();
if ( e->IsQuadratic() )
{
- quadAdjacentElems[ e->GetType() ].insert( e );
- continue;
+ bool alreadyOK;
+ switch ( e->GetEntityType() ) {
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break;
+ case SMDSEntity_BiQuad_Quadrangle:
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break;
+ default: alreadyOK = true;
+ }
+ if ( alreadyOK )
+ {
+ quadAdjacentElems[ e->GetType() ].insert( e );
+ continue;
+ }
}
if ( e->GetType() >= elemType )
{
SMESH_MesherHelper helper(*myMesh);
helper.SetIsQuadratic( true );
+ helper.SetIsBiQuadratic( theToBiQuad );
// add links of quadratic adjacent elements to the helper
helper.AddTLinks( static_cast< const SMDS_MeshVolume*> (*eIt) );
}
- // make quadratic elements instead of linear ones
+ // make quadratic (or bi-tri-quadratic) elements instead of linear ones
- SMESHDS_Mesh* meshDS = GetMeshDS();
+ SMESHDS_Mesh* meshDS = GetMeshDS();
SMESHDS_SubMesh* smDS = 0;
for ( eIt = theElements.begin(); eIt != theElements.end(); ++eIt )
{
const SMDS_MeshElement* elem = *eIt;
- if( elem->IsQuadratic() || elem->NbNodes() < 2 || elem->IsPoly() )
+ if( elem->NbNodes() < 2 || elem->IsPoly() )
continue;
- int id = elem->GetID();
- SMDSAbs_ElementType type = elem->GetType();
+ if ( elem->IsQuadratic() )
+ {
+ bool alreadyOK;
+ switch ( elem->GetEntityType() ) {
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break;
+ case SMDSEntity_BiQuad_Quadrangle:
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break;
+ default: alreadyOK = true;
+ }
+ if ( alreadyOK ) continue;
+ }
+
+ const SMDSAbs_ElementType type = elem->GetType();
+ const int id = elem->GetID();
+ const int nbNodes = elem->NbCornerNodes();
vector<const SMDS_MeshNode *> nodes ( elem->begin_nodes(), elem->end_nodes());
if ( !smDS || !smDS->Contains( elem ))
meshDS->RemoveFreeElement(elem, smDS, /*fromGroups=*/false);
SMDS_MeshElement * newElem = 0;
- switch( nodes.size() )
+ switch( nbNodes )
{
- case 4: // cases for most multiple element types go first (for optimization)
+ case 4: // cases for most frequently used element types go first (for optimization)
if ( type == SMDSAbs_Volume )
newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
else
if ( !theForce3d && !getenv("NO_FixQuadraticElements"))
{ // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
helper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
- helper.FixQuadraticElements();
+ helper.FixQuadraticElements( myError );
}
}
SMESHDS_Mesh* aMesh = GetMeshDS();
// TODO algoritm not OK with vtkUnstructuredGrid: 2 meshes can't share nodes
//SMDS_Mesh aTmpFacesMesh; // try to use the same mesh
- set<const SMDS_MeshElement*> faceSet1, faceSet2;
+ TIDSortedElemSet faceSet1, faceSet2;
set<const SMDS_MeshElement*> volSet1, volSet2;
set<const SMDS_MeshNode*> nodeSet1, nodeSet2;
- set<const SMDS_MeshElement*> * faceSetPtr[] = { &faceSet1, &faceSet2 };
- set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
+ TIDSortedElemSet * faceSetPtr[] = { &faceSet1, &faceSet2 };
+ set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
set<const SMDS_MeshNode*> * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
- TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
+ TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
int iSide, iFace, iNode;
list<const SMDS_MeshElement* > tempFaceList;
for ( iSide = 0; iSide < 2; iSide++ ) {
set<const SMDS_MeshNode*> * nodeSet = nodeSetPtr[ iSide ];
- TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
- set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
+ TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
+ TIDSortedElemSet * faceSet = faceSetPtr[ iSide ];
set<const SMDS_MeshElement*> * volSet = volSetPtr [ iSide ];
set<const SMDS_MeshElement*>::iterator vIt;
TIDSortedElemSet::iterator eIt;
bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
if ( isNewFace ) {
// no such a face is given but it still can exist, check it
- if ( nbNodes == 3 ) {
- aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2] );
- }
- else if ( nbNodes == 4 ) {
- aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
- }
- else {
- vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
- aFreeFace = aMesh->FindFace(poly_nodes);
- }
+ vector<const SMDS_MeshNode *> nodes ( fNodes, fNodes + nbNodes);
+ aFreeFace = aMesh->FindElement( nodes, SMDSAbs_Face, /*noMedium=*/false );
}
if ( !aFreeFace ) {
// create a temporary face
//aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
aFreeFace = aMesh->AddPolygonalFace(poly_nodes);
}
+ if ( aFreeFace )
+ tempFaceList.push_back( aFreeFace );
}
- if ( aFreeFace ) {
+
+ if ( aFreeFace )
freeFaceList.push_back( aFreeFace );
- tempFaceList.push_back( aFreeFace );
- }
} // loop on faces of a volume
- // choose one of several free faces
- // --------------------------------------
+ // choose one of several free faces of a volume
+ // --------------------------------------------
if ( freeFaceList.size() > 1 ) {
// choose a face having max nb of nodes shared by other elems of a side
- int maxNbNodes = -1/*, nbExcludedFaces = 0*/;
+ int maxNbNodes = -1;
list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
while ( fIt != freeFaceList.end() ) { // loop on free faces
int nbSharedNodes = 0;
SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
while ( invElemIt->more() ) {
const SMDS_MeshElement* e = invElemIt->next();
- if ( faceSet->find( e ) != faceSet->end() )
- nbSharedNodes++;
- if ( elemSet->find( e ) != elemSet->end() )
- nbSharedNodes++;
+ nbSharedNodes += faceSet->count( e );
+ nbSharedNodes += elemSet->count( e );
}
}
- if ( nbSharedNodes >= maxNbNodes ) {
+ if ( nbSharedNodes > maxNbNodes ) {
maxNbNodes = nbSharedNodes;
+ freeFaceList.erase( freeFaceList.begin(), fIt++ );
+ }
+ else if ( nbSharedNodes == maxNbNodes ) {
fIt++;
}
- else
+ else {
freeFaceList.erase( fIt++ ); // here fIt++ occurs before erase
+ }
}
if ( freeFaceList.size() > 1 )
{
TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
if ( theFirstNode1 != theFirstNode2 )
- nReplaceMap.insert( TNodeNodeMap::value_type( theFirstNode1, theFirstNode2 ));
+ nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
if ( theSecondNode1 != theSecondNode2 )
- nReplaceMap.insert( TNodeNodeMap::value_type( theSecondNode1, theSecondNode2 ));
+ nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
LinkID_Gen aLinkID_Gen( GetMeshDS() );
set< long > linkIdSet; // links to process
// loop on links in linkList; find faces by links and append links
// of the found faces to linkList
list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
- for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
+ for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ )
+ {
NLink link[] = { *linkIt[0], *linkIt[1] };
long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
- if ( linkIdSet.find( linkID ) == linkIdSet.end() )
+ if ( !linkIdSet.count( linkID ) )
continue;
// by links, find faces in the face sets,
// ---------------------------------------------------------------
const SMDS_MeshElement* face[] = { 0, 0 };
- //const SMDS_MeshNode* faceNodes[ 2 ][ 5 ];
- vector<const SMDS_MeshNode*> fnodes1(9);
- vector<const SMDS_MeshNode*> fnodes2(9);
- //const SMDS_MeshNode* notLinkNodes[ 2 ][ 2 ] = {{ 0, 0 },{ 0, 0 }} ;
- vector<const SMDS_MeshNode*> notLinkNodes1(6);
- vector<const SMDS_MeshNode*> notLinkNodes2(6);
+ vector<const SMDS_MeshNode*> fnodes[2];
int iLinkNode[2][2];
+ TIDSortedElemSet avoidSet;
for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
const SMDS_MeshNode* n1 = link[iSide].first;
const SMDS_MeshNode* n2 = link[iSide].second;
- set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
- set< const SMDS_MeshElement* > fMap;
- for ( int i = 0; i < 2; i++ ) { // loop on 2 nodes of a link
- const SMDS_MeshNode* n = i ? n1 : n2; // a node of a link
- SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
- while ( fIt->more() ) { // loop on faces sharing a node
- const SMDS_MeshElement* f = fIt->next();
- if (faceSet->find( f ) != faceSet->end() && // f is in face set
- ! fMap.insert( f ).second ) // f encounters twice
- {
- if ( face[ iSide ] ) {
- MESSAGE( "2 faces per link " );
- aResult = iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES;
- break;
- }
- face[ iSide ] = f;
- faceSet->erase( f );
- // get face nodes and find ones of a link
- iNode = 0;
- int nbl = -1;
- if(f->IsPoly()) {
- if(iSide==0) {
- fnodes1.resize(f->NbNodes()+1);
- notLinkNodes1.resize(f->NbNodes()-2);
- }
- else {
- fnodes2.resize(f->NbNodes()+1);
- notLinkNodes2.resize(f->NbNodes()-2);
- }
- }
- if(!f->IsQuadratic()) {
- SMDS_ElemIteratorPtr nIt = f->nodesIterator();
- while ( nIt->more() ) {
- const SMDS_MeshNode* n =
- static_cast<const SMDS_MeshNode*>( nIt->next() );
- if ( n == n1 ) {
- iLinkNode[ iSide ][ 0 ] = iNode;
- }
- else if ( n == n2 ) {
- iLinkNode[ iSide ][ 1 ] = iNode;
- }
- //else if ( notLinkNodes[ iSide ][ 0 ] )
- // notLinkNodes[ iSide ][ 1 ] = n;
- //else
- // notLinkNodes[ iSide ][ 0 ] = n;
- else {
- nbl++;
- if(iSide==0)
- notLinkNodes1[nbl] = n;
- //notLinkNodes1.push_back(n);
- else
- notLinkNodes2[nbl] = n;
- //notLinkNodes2.push_back(n);
- }
- //faceNodes[ iSide ][ iNode++ ] = n;
- if(iSide==0) {
- fnodes1[iNode++] = n;
- }
- else {
- fnodes2[iNode++] = n;
- }
- }
- }
- else { // f->IsQuadratic()
- const SMDS_VtkFace* F =
- dynamic_cast<const SMDS_VtkFace*>(f);
- if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
- // use special nodes iterator
- SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
- while ( anIter->more() ) {
- const SMDS_MeshNode* n =
- static_cast<const SMDS_MeshNode*>( anIter->next() );
- if ( n == n1 ) {
- iLinkNode[ iSide ][ 0 ] = iNode;
- }
- else if ( n == n2 ) {
- iLinkNode[ iSide ][ 1 ] = iNode;
- }
- else {
- nbl++;
- if(iSide==0) {
- notLinkNodes1[nbl] = n;
- }
- else {
- notLinkNodes2[nbl] = n;
- }
- }
- if(iSide==0) {
- fnodes1[iNode++] = n;
- }
- else {
- fnodes2[iNode++] = n;
- }
- }
- }
- //faceNodes[ iSide ][ iNode ] = faceNodes[ iSide ][ 0 ];
- if(iSide==0) {
- fnodes1[iNode] = fnodes1[0];
- }
- else {
- fnodes2[iNode] = fnodes1[0];
- }
- }
+ //cout << "Side " << iSide << " ";
+ //cout << "L( " << n1->GetID() << ", " << n2->GetID() << " ) " << endl;
+ // find a face by two link nodes
+ face[ iSide ] = FindFaceInSet( n1, n2, *faceSetPtr[ iSide ], avoidSet,
+ &iLinkNode[iSide][0], &iLinkNode[iSide][1] );
+ if ( face[ iSide ])
+ {
+ //cout << " F " << face[ iSide]->GetID() <<endl;
+ faceSetPtr[ iSide ]->erase( face[ iSide ]);
+ // put face nodes to fnodes
+ if ( face[ iSide ]->IsQuadratic() )
+ {
+ // use interlaced nodes iterator
+ const SMDS_VtkFace* F = dynamic_cast<const SMDS_VtkFace*>( face[ iSide ]);
+ if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
+ SMDS_ElemIteratorPtr nIter = F->interlacedNodesElemIterator();
+ while ( nIter->more() )
+ fnodes[ iSide ].push_back( cast2Node( nIter->next() ));
+ }
+ else
+ {
+ fnodes[ iSide ].assign( face[ iSide ]->begin_nodes(),
+ face[ iSide ]->end_nodes() );
}
+ fnodes[ iSide ].push_back( fnodes[ iSide ].front());
}
}
else {
aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
}
- break; // do not return because it s necessary to remove tmp faces
+ break; // do not return because it's necessary to remove tmp faces
}
// set nodes to merge
// -------------------
if ( face[0] && face[1] ) {
- int nbNodes = face[0]->NbNodes();
+ const int nbNodes = face[0]->NbNodes();
if ( nbNodes != face[1]->NbNodes() ) {
MESSAGE("Diff nb of face nodes");
aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
break; // do not return because it s necessary to remove tmp faces
}
- bool reverse[] = { false, false }; // order of notLinkNodes of quadrangle
- if ( nbNodes == 3 ) {
- //nReplaceMap.insert( TNodeNodeMap::value_type
- // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
- nReplaceMap.insert( TNodeNodeMap::value_type
- ( notLinkNodes1[0], notLinkNodes2[0] ));
+ bool reverse[] = { false, false }; // order of nodes in the link
+ for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
+ // analyse link orientation in faces
+ int i1 = iLinkNode[ iSide ][ 0 ];
+ int i2 = iLinkNode[ iSide ][ 1 ];
+ reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
}
- else {
- for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
- // analyse link orientation in faces
- int i1 = iLinkNode[ iSide ][ 0 ];
- int i2 = iLinkNode[ iSide ][ 1 ];
- reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
- // if notLinkNodes are the first and the last ones, then
- // their order does not correspond to the link orientation
- if (( i1 == 1 && i2 == 2 ) ||
- ( i1 == 2 && i2 == 1 ))
- reverse[ iSide ] = !reverse[ iSide ];
- }
- if ( reverse[0] == reverse[1] ) {
- //nReplaceMap.insert( TNodeNodeMap::value_type
- // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
- //nReplaceMap.insert( TNodeNodeMap::value_type
- // ( notLinkNodes[0][1], notLinkNodes[1][1] ));
- for(int nn=0; nn<nbNodes-2; nn++) {
- nReplaceMap.insert( TNodeNodeMap::value_type
- ( notLinkNodes1[nn], notLinkNodes2[nn] ));
- }
- }
- else {
- //nReplaceMap.insert( TNodeNodeMap::value_type
- // ( notLinkNodes[0][0], notLinkNodes[1][1] ));
- //nReplaceMap.insert( TNodeNodeMap::value_type
- // ( notLinkNodes[0][1], notLinkNodes[1][0] ));
- for(int nn=0; nn<nbNodes-2; nn++) {
- nReplaceMap.insert( TNodeNodeMap::value_type
- ( notLinkNodes1[nn], notLinkNodes2[nbNodes-3-nn] ));
- }
- }
+ int di1 = reverse[0] ? -1 : +1, i1 = iLinkNode[0][1] + di1;
+ int di2 = reverse[1] ? -1 : +1, i2 = iLinkNode[1][1] + di2;
+ for ( int i = nbNodes - 2; i > 0; --i, i1 += di1, i2 += di2 )
+ {
+ nReplaceMap.insert ( make_pair ( fnodes[0][ ( i1 + nbNodes ) % nbNodes ],
+ fnodes[1][ ( i2 + nbNodes ) % nbNodes ]));
}
// add other links of the faces to linkList
// -----------------------------------------
- //const SMDS_MeshNode** nodes = faceNodes[ 0 ];
for ( iNode = 0; iNode < nbNodes; iNode++ ) {
- //linkID = aLinkID_Gen.GetLinkID( nodes[iNode], nodes[iNode+1] );
- linkID = aLinkID_Gen.GetLinkID( fnodes1[iNode], fnodes1[iNode+1] );
+ linkID = aLinkID_Gen.GetLinkID( fnodes[0][iNode], fnodes[0][iNode+1] );
pair< set<long>::iterator, bool > iter_isnew = linkIdSet.insert( linkID );
if ( !iter_isnew.second ) { // already in a set: no need to process
linkIdSet.erase( iter_isnew.first );
}
else // new in set == encountered for the first time: add
{
- //const SMDS_MeshNode* n1 = nodes[ iNode ];
- //const SMDS_MeshNode* n2 = nodes[ iNode + 1];
- const SMDS_MeshNode* n1 = fnodes1[ iNode ];
- const SMDS_MeshNode* n2 = fnodes1[ iNode + 1];
+ const SMDS_MeshNode* n1 = fnodes[0][ iNode ];
+ const SMDS_MeshNode* n2 = fnodes[0][ iNode + 1];
linkList[0].push_back ( NLink( n1, n2 ));
linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
}
}
} // 2 faces found
+
+ if ( faceSetPtr[0]->empty() || faceSetPtr[1]->empty() )
+ break;
+
} // loop on link lists
if ( aResult == SEW_OK &&
- ( linkIt[0] != linkList[0].end() ||
+ ( //linkIt[0] != linkList[0].end() ||
!faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
" " << (faceSetPtr[1]->empty()));
// 3. Replace nodes in elements of the side 1 and remove replaced nodes
// ====================================================================
- // delete temporary faces: they are in reverseElements of actual nodes
+ // delete temporary faces
// SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
// while ( tmpFaceIt->more() )
// aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
-// list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
-// for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
-// aMesh->RemoveElement(*tmpFaceIt);
+ list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
+ for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
+ aMesh->RemoveElement(*tmpFaceIt);
if ( aResult != SEW_OK)
return aResult;
\param theElems - the list of elements (edges or faces) to be replicated
The nodes for duplication could be found from these elements
\param theNodesNot - list of nodes to NOT replicate
- \param theAffectedElems - the list of elements (cells and edges) to which the
+ \param theAffectedElems - the list of elements (cells and edges) to which the
replicated nodes should be associated to.
\return TRUE if operation has been completed successfully, FALSE otherwise
*/
std::vector<const SMDS_MeshNode*> newNodes( anElem->NbNodes() );
SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
int ind = 0;
- while ( anIter->more() )
- {
+ while ( anIter->more() )
+ {
SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
SMDS_MeshNode* aNewNode = aCurrNode;
/*!
\brief Creates a hole in a mesh by doubling the nodes of some particular elements
\param theNodes - identifiers of nodes to be doubled
- \param theModifiedElems - identifiers of elements to be updated by the new (doubled)
- nodes. If list of element identifiers is empty then nodes are doubled but
+ \param theModifiedElems - identifiers of elements to be updated by the new (doubled)
+ nodes. If list of element identifiers is empty then nodes are doubled but
they not assigned to elements
\return TRUE if operation has been completed successfully, FALSE otherwise
*/
//================================================================================
-bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
+bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
const std::list< int >& theListOfModifiedElems )
{
MESSAGE("DoubleNodes");
std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> > anElemToNodes;
std::list< int >::const_iterator anElemIter;
- for ( anElemIter = theListOfModifiedElems.begin();
+ for ( anElemIter = theListOfModifiedElems.begin();
anElemIter != theListOfModifiedElems.end(); ++anElemIter )
{
int aCurr = *anElemIter;
SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
int ind = 0;
- while ( anIter->more() )
- {
+ while ( anIter->more() )
+ {
SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() )
{
anElemToNodes[ anElem ] = aNodeArr;
}
- // Change nodes of elements
+ // Change nodes of elements
std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> >::iterator
anElemToNodesIter = anElemToNodes.begin();
};
}
+//================================================================================
+/*!
+ \brief Identify the elements that will be affected by node duplication (actual duplication is not performed.
+ This method is the first step of DoubleNodeElemGroupsInRegion.
+ \param theElems - list of groups of elements (edges or faces) to be replicated
+ \param theNodesNot - list of groups of nodes not to replicated
+ \param theShape - shape to detect affected elements (element which geometric center
+ located on or inside shape).
+ The replicated nodes should be associated to affected elements.
+ \return groups of affected elements
+ \sa DoubleNodeElemGroupsInRegion()
+ */
+//================================================================================
+
+bool SMESH_MeshEditor::AffectedElemGroupsInRegion( const TIDSortedElemSet& theElems,
+ const TIDSortedElemSet& theNodesNot,
+ const TopoDS_Shape& theShape,
+ TIDSortedElemSet& theAffectedElems)
+{
+ if ( theShape.IsNull() )
+ return false;
+
+ const double aTol = Precision::Confusion();
+ auto_ptr< BRepClass3d_SolidClassifier> bsc3d;
+ auto_ptr<_FaceClassifier> aFaceClassifier;
+ if ( theShape.ShapeType() == TopAbs_SOLID )
+ {
+ bsc3d.reset( new BRepClass3d_SolidClassifier(theShape));;
+ bsc3d->PerformInfinitePoint(aTol);
+ }
+ else if (theShape.ShapeType() == TopAbs_FACE )
+ {
+ aFaceClassifier.reset( new _FaceClassifier(TopoDS::Face(theShape)));
+ }
+
+ // iterates on indicated elements and get elements by back references from their nodes
+ TIDSortedElemSet::const_iterator elemItr = theElems.begin();
+ for ( ; elemItr != theElems.end(); ++elemItr )
+ {
+ SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
+ if (!anElem)
+ continue;
+
+ SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
+ while ( nodeItr->more() )
+ {
+ const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
+ if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
+ continue;
+ SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
+ while ( backElemItr->more() )
+ {
+ const SMDS_MeshElement* curElem = backElemItr->next();
+ if ( curElem && theElems.find(curElem) == theElems.end() &&
+ ( bsc3d.get() ?
+ isInside( curElem, *bsc3d, aTol ) :
+ isInside( curElem, *aFaceClassifier, aTol )))
+ theAffectedElems.insert( curElem );
+ }
+ }
+ }
+ return true;
+}
+
//================================================================================
/*!
\brief Creates a hole in a mesh by doubling the nodes of some particular elements
// --- build a map (face to duplicate --> volume to modify)
// with all the faces shared by 2 domains (group of elements)
// and corresponding volume of this domain, for each shared face.
- // a volume has a face shared by 2 domains if it has a neighbor which is not in is domain.
+ // a volume has a face shared by 2 domains if it has a neighbor which is not in his domain.
+ //MESSAGE("Domain " << idom);
const TIDSortedElemSet& domain = theElems[idom];
TIDSortedElemSet::const_iterator elemItr = domain.begin();
for (; elemItr != domain.end(); ++elemItr)
if (!anElem)
continue;
int vtkId = anElem->getVtkId();
+ //MESSAGE(" vtkId " << vtkId << " smdsId " << anElem->GetID());
int neighborsVtkIds[NBMAXNEIGHBORS];
int downIds[NBMAXNEIGHBORS];
unsigned char downTypes[NBMAXNEIGHBORS];
{
faceDomains[face][idom] = vtkId; // volume associated to face in this domain
celldom[vtkId] = idom;
+ //MESSAGE(" cell with a border " << vtkId << " domain " << idom);
}
}
}
for (int idomain = 0; idomain < theElems.size(); idomain++)
{
+ //MESSAGE("Domain " << idomain);
const TIDSortedElemSet& domain = theElems[idomain];
itface = faceDomains.begin();
for (; itface != faceDomains.end(); ++itface)
{
int oldId = *itn;
//MESSAGE(" node " << oldId);
- std::set<int> cells;
- cells.clear();
vtkCellLinks::Link l = grid->GetCellLinks()->GetLink(oldId);
for (int i=0; i<l.ncells; i++)
{
//no cells created after BuildDownWardConnectivity
}
DownIdType aCell(downId, vtkType);
- if (celldom.count(vtkId))
- continue;
+ if (!cellDomains.count(aCell))
+ cellDomains[aCell] = emptyMap; // create an empty entry for cell
cellDomains[aCell][idomain] = vtkId;
celldom[vtkId] = idomain;
+ //MESSAGE(" cell " << vtkId << " domain " << idomain);
}
}
}
// the value is the ordered domain ids. (more than 4 domains not taken into account)
std::map<std::vector<int>, std::vector<int> > edgesMultiDomains; // nodes of edge --> ordered domains
- std::map<int, std::vector<int> > mutipleNodes; // nodes muti domains with domain order
+ std::map<int, std::vector<int> > mutipleNodes; // nodes multi domains with domain order
+ std::map<int, std::vector<int> > mutipleNodesToFace; // nodes multi domains with domain order to transform in Face (junction between 3 or more 2D domains)
for (int idomain = 0; idomain < theElems.size(); idomain++)
{
std::set<int> oldNodes;
oldNodes.clear();
grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
- bool isMultipleDetected = false;
std::set<int>::iterator itn = oldNodes.begin();
for (; itn != oldNodes.end(); ++itn)
{
int oldId = *itn;
- //MESSAGE(" node " << oldId);
+ //MESSAGE("-+-+-a node " << oldId);
if (!nodeDomains.count(oldId))
nodeDomains[oldId] = emptyMap; // create an empty entry for node
if (nodeDomains[oldId].empty())
- nodeDomains[oldId][idomain] = oldId; // keep the old node in the first domain
+ {
+ nodeDomains[oldId][idomain] = oldId; // keep the old node in the first domain
+ //MESSAGE("-+-+-b oldNode " << oldId << " domain " << idomain);
+ }
std::map<int, int>::iterator itdom = domvol.begin();
for (; itdom != domvol.end(); ++itdom)
{
{
vector<int> orderedDoms;
//MESSAGE("multiple node " << oldId);
- isMultipleDetected =true;
if (mutipleNodes.count(oldId))
orderedDoms = mutipleNodes[oldId];
else
SMDS_MeshNode *newNode = meshDS->AddNode(coords[0], coords[1], coords[2]);
int newId = newNode->getVtkId();
nodeDomains[oldId][idom] = newId; // cloned node for other domains
- //MESSAGE(" newNode " << newId << " oldNode " << oldId << " size=" <<nodeDomains[oldId].size());
- }
- if (nodeDomains[oldId].size() >= 3)
- {
- //MESSAGE("confirm multiple node " << oldId);
- isMultipleDetected =true;
+ //MESSAGE("-+-+-c oldNode " << oldId << " domain " << idomain << " newNode " << newId << " domain " << idom << " size=" <<nodeDomains[oldId].size());
}
}
}
- if (isMultipleDetected) // check if an edge of the face is shared between 3 or more domains
+ }
+ }
+
+ for (int idomain = 0; idomain < theElems.size(); idomain++)
+ {
+ itface = faceDomains.begin();
+ for (; itface != faceDomains.end(); ++itface)
+ {
+ std::map<int, int> domvol = itface->second;
+ if (!domvol.count(idomain))
+ continue;
+ DownIdType face = itface->first;
+ //MESSAGE(" --- face " << face.cellId);
+ std::set<int> oldNodes;
+ oldNodes.clear();
+ grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+ int nbMultipleNodes = 0;
+ std::set<int>::iterator itn = oldNodes.begin();
+ for (; itn != oldNodes.end(); ++itn)
+ {
+ int oldId = *itn;
+ if (mutipleNodes.count(oldId))
+ nbMultipleNodes++;
+ }
+ if (nbMultipleNodes > 1) // check if an edge of the face is shared between 3 or more domains
{
//MESSAGE("multiple Nodes detected on a shared face");
int downId = itface->first.cellId;
unsigned char cellType = itface->first.cellType;
- int nbEdges = grid->getDownArray(cellType)->getNumberOfDownCells(downId);
- const int *downEdgeIds = grid->getDownArray(cellType)->getDownCells(downId);
- const unsigned char* edgeType = grid->getDownArray(cellType)->getDownTypes(downId);
- for (int ie =0; ie < nbEdges; ie++)
+ // --- shared edge or shared face ?
+ if ((cellType == VTK_LINE) || (cellType == VTK_QUADRATIC_EDGE)) // shared edge (between two faces)
{
int nodes[3];
- int nbNodes = grid->getDownArray(edgeType[ie])->getNodes(downEdgeIds[ie], nodes);
- if (mutipleNodes.count(nodes[0]) && mutipleNodes.count(nodes[nbNodes-1]))
+ int nbNodes = grid->getDownArray(cellType)->getNodes(downId, nodes);
+ for (int i=0; i< nbNodes; i=i+nbNodes-1) // i=0 , i=nbNodes-1
+ if (mutipleNodes.count(nodes[i]))
+ if (!mutipleNodesToFace.count(nodes[i]))
+ mutipleNodesToFace[nodes[i]] = mutipleNodes[nodes[i]];
+ }
+ else // shared face (between two volumes)
+ {
+ int nbEdges = grid->getDownArray(cellType)->getNumberOfDownCells(downId);
+ const int* downEdgeIds = grid->getDownArray(cellType)->getDownCells(downId);
+ const unsigned char* edgeType = grid->getDownArray(cellType)->getDownTypes(downId);
+ for (int ie =0; ie < nbEdges; ie++)
{
- vector<int> vn0 = mutipleNodes[nodes[0]];
- vector<int> vn1 = mutipleNodes[nodes[nbNodes - 1]];
- sort( vn0.begin(), vn0.end() );
- sort( vn1.begin(), vn1.end() );
- if (vn0 == vn1)
+ int nodes[3];
+ int nbNodes = grid->getDownArray(edgeType[ie])->getNodes(downEdgeIds[ie], nodes);
+ if (mutipleNodes.count(nodes[0]) && mutipleNodes.count(nodes[nbNodes-1]))
{
- //MESSAGE(" detect edgesMultiDomains " << nodes[0] << " " << nodes[nbNodes - 1]);
- double *coords = grid->GetPoint(nodes[0]);
- gp_Pnt p0(coords[0], coords[1], coords[2]);
- coords = grid->GetPoint(nodes[nbNodes - 1]);
- gp_Pnt p1(coords[0], coords[1], coords[2]);
- gp_Pnt gref;
- int vtkVolIds[1000]; // an edge can belong to a lot of volumes
- map<int, SMDS_VtkVolume*> domvol; // domain --> a volume with the edge
- map<int, double> angleDom; // oriented angles between planes defined by edge and volume centers
- int nbvol = grid->GetParentVolumes(vtkVolIds, downEdgeIds[ie], edgeType[ie]);
- for (int id=0; id < vn0.size(); id++)
+ vector<int> vn0 = mutipleNodes[nodes[0]];
+ vector<int> vn1 = mutipleNodes[nodes[nbNodes - 1]];
+ vector<int> doms;
+ for (int i0 = 0; i0 < vn0.size(); i0++)
+ for (int i1 = 0; i1 < vn1.size(); i1++)
+ if (vn0[i0] == vn1[i1])
+ doms.push_back(vn0[i0]);
+ if (doms.size() >2)
{
- int idom = vn0[id];
- for (int ivol=0; ivol<nbvol; ivol++)
+ //MESSAGE(" detect edgesMultiDomains " << nodes[0] << " " << nodes[nbNodes - 1]);
+ double *coords = grid->GetPoint(nodes[0]);
+ gp_Pnt p0(coords[0], coords[1], coords[2]);
+ coords = grid->GetPoint(nodes[nbNodes - 1]);
+ gp_Pnt p1(coords[0], coords[1], coords[2]);
+ gp_Pnt gref;
+ int vtkVolIds[1000]; // an edge can belong to a lot of volumes
+ map<int, SMDS_VtkVolume*> domvol; // domain --> a volume with the edge
+ map<int, double> angleDom; // oriented angles between planes defined by edge and volume centers
+ int nbvol = grid->GetParentVolumes(vtkVolIds, downEdgeIds[ie], edgeType[ie]);
+ for (int id=0; id < doms.size(); id++)
{
- int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
- SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
- if (theElems[idom].count(elem))
+ int idom = doms[id];
+ for (int ivol=0; ivol<nbvol; ivol++)
{
- SMDS_VtkVolume* svol = dynamic_cast<SMDS_VtkVolume*>(elem);
- domvol[idom] = svol;
- //MESSAGE(" domain " << idom << " volume " << elem->GetID());
- double values[3];
- vtkIdType npts = 0;
- vtkIdType* pts = 0;
- grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
- SMDS_VtkVolume::gravityCenter(grid, pts, npts, values);
- if (id ==0)
- {
- gref.SetXYZ(gp_XYZ(values[0], values[1], values[2]));
- angleDom[idom] = 0;
- }
- else
+ int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
+ SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
+ if (theElems[idom].count(elem))
{
- gp_Pnt g(values[0], values[1], values[2]);
- angleDom[idom] = OrientedAngle(p0, p1, gref, g); // -pi<angle<+pi
- //MESSAGE(" angle=" << angleDom[idom]);
+ SMDS_VtkVolume* svol = dynamic_cast<SMDS_VtkVolume*>(elem);
+ domvol[idom] = svol;
+ //MESSAGE(" domain " << idom << " volume " << elem->GetID());
+ double values[3];
+ vtkIdType npts = 0;
+ vtkIdType* pts = 0;
+ grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
+ SMDS_VtkVolume::gravityCenter(grid, pts, npts, values);
+ if (id ==0)
+ {
+ gref.SetXYZ(gp_XYZ(values[0], values[1], values[2]));
+ angleDom[idom] = 0;
+ }
+ else
+ {
+ gp_Pnt g(values[0], values[1], values[2]);
+ angleDom[idom] = OrientedAngle(p0, p1, gref, g); // -pi<angle<+pi
+ //MESSAGE(" angle=" << angleDom[idom]);
+ }
+ break;
}
- break;
}
}
+ map<double, int> sortedDom; // sort domains by angle
+ for (map<int, double>::iterator ia = angleDom.begin(); ia != angleDom.end(); ++ia)
+ sortedDom[ia->second] = ia->first;
+ vector<int> vnodes;
+ vector<int> vdom;
+ for (map<double, int>::iterator ib = sortedDom.begin(); ib != sortedDom.end(); ++ib)
+ {
+ vdom.push_back(ib->second);
+ //MESSAGE(" ordered domain " << ib->second << " angle " << ib->first);
+ }
+ for (int ino = 0; ino < nbNodes; ino++)
+ vnodes.push_back(nodes[ino]);
+ edgesMultiDomains[vnodes] = vdom; // nodes vector --> ordered domains
}
- map<double, int> sortedDom; // sort domains by angle
- for (map<int, double>::iterator ia = angleDom.begin(); ia != angleDom.end(); ++ia)
- sortedDom[ia->second] = ia->first;
- vector<int> vnodes;
- vector<int> vdom;
- for (map<double, int>::iterator ib = sortedDom.begin(); ib != sortedDom.end(); ++ib)
- {
- vdom.push_back(ib->second);
- //MESSAGE(" ordered domain " << ib->second << " angle " << ib->first);
- }
- for (int ino = 0; ino < nbNodes; ino++)
- vnodes.push_back(nodes[ino]);
- edgesMultiDomains[vnodes] = vdom; // nodes vector --> ordered domains
}
}
}
if (createJointElems)
{
+ int idg;
+ string joints2DName = "joints2D";
+ mapOfJunctionGroups[joints2DName] = this->myMesh->AddGroup(SMDSAbs_Face, joints2DName.c_str(), idg);
+ SMESHDS_Group *joints2DGrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[joints2DName]->GetGroupDS());
+ string joints3DName = "joints3D";
+ mapOfJunctionGroups[joints3DName] = this->myMesh->AddGroup(SMDSAbs_Volume, joints3DName.c_str(), idg);
+ SMESHDS_Group *joints3DGrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[joints3DName]->GetGroupDS());
+
itface = faceDomains.begin();
for (; itface != faceDomains.end(); ++itface)
{
int vtkVolId = itdom->second;
itdom++;
int dom2 = itdom->first;
- SMDS_MeshVolume *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
+ SMDS_MeshCell *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
nodeQuadDomains);
stringstream grpname;
grpname << "j_";
grpname << dom1 << "_" << dom2;
else
grpname << dom2 << "_" << dom1;
- int idg;
string namegrp = grpname.str();
if (!mapOfJunctionGroups.count(namegrp))
- mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+ mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(vol->GetType(), namegrp.c_str(), idg);
SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
if (sgrp)
sgrp->Add(vol->GetID());
+ if (vol->GetType() == SMDSAbs_Volume)
+ joints3DGrp->Add(vol->GetID());
+ else if (vol->GetType() == SMDSAbs_Face)
+ joints2DGrp->Add(vol->GetID());
}
}
// --- create volumes on multiple domain intersection if requested
+ // iterate on mutipleNodesToFace
// iterate on edgesMultiDomains
if (createJointElems)
{
+ // --- iterate on mutipleNodesToFace
+
+ std::map<int, std::vector<int> >::iterator itn = mutipleNodesToFace.begin();
+ for (; itn != mutipleNodesToFace.end(); ++itn)
+ {
+ int node = itn->first;
+ vector<int> orderDom = itn->second;
+ vector<vtkIdType> orderedNodes;
+ for (int idom = 0; idom <orderDom.size(); idom++)
+ orderedNodes.push_back( nodeDomains[node][orderDom[idom]] );
+ SMDS_MeshFace* face = this->GetMeshDS()->AddFaceFromVtkIds(orderedNodes);
+
+ stringstream grpname;
+ grpname << "m2j_";
+ grpname << 0 << "_" << 0;
+ int idg;
+ string namegrp = grpname.str();
+ if (!mapOfJunctionGroups.count(namegrp))
+ mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Face, namegrp.c_str(), idg);
+ SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+ if (sgrp)
+ sgrp->Add(face->GetID());
+ }
+
+ // --- iterate on edgesMultiDomains
+
std::map<std::vector<int>, std::vector<int> >::iterator ite = edgesMultiDomains.begin();
for (; ite != edgesMultiDomains.end(); ++ite)
{
orderedNodes.push_back( nodeDomains[nodes[ino]][orderDom[idom]] );
SMDS_MeshVolume* vol = this->GetMeshDS()->AddVolumeFromVtkIds(orderedNodes);
- stringstream grpname;
- grpname << "mj_";
- grpname << 0 << "_" << 0;
int idg;
- string namegrp = grpname.str();
+ string namegrp = "jointsMultiples";
if (!mapOfJunctionGroups.count(namegrp))
mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
}
else
{
- //MESSAGE("Quadratic multiple joints not implemented");
+ INFOS("Quadratic multiple joints not implemented");
// TODO quadratic nodes
}
}
return true;
}
+/*!
+ * \brief identify all the elements around a geom shape, get the faces delimiting the hole
+ * Build groups of volume to remove, groups of faces to replace on the skin of the object,
+ * groups of faces to remove inside the object, (idem edges).
+ * Build ordered list of nodes at the border of each group of faces to replace (to be used to build a geom subshape)
+ */
+void SMESH_MeshEditor::CreateHoleSkin(double radius,
+ const TopoDS_Shape& theShape,
+ SMESH_NodeSearcher* theNodeSearcher,
+ const char* groupName,
+ std::vector<double>& nodesCoords,
+ std::vector<std::vector<int> >& listOfListOfNodes)
+{
+ MESSAGE("--------------------------------");
+ MESSAGE("SMESH_MeshEditor::CreateHoleSkin");
+ MESSAGE("--------------------------------");
+
+ // --- zone of volumes to remove is given :
+ // 1 either by a geom shape (one or more vertices) and a radius,
+ // 2 either by a group of nodes (representative of the shape)to use with the radius,
+ // 3 either by a group of nodes where all the elements build on one of this nodes are to remove,
+ // In the case 2, the group of nodes is an external group of nodes from another mesh,
+ // In the case 3, the group of nodes is an internal group of the mesh (obtained for instance by a filter),
+ // defined by it's name.
+
+ SMESHDS_GroupBase* groupDS = 0;
+ SMESH_Mesh::GroupIteratorPtr groupIt = this->myMesh->GetGroups();
+ while ( groupIt->more() )
+ {
+ groupDS = 0;
+ SMESH_Group * group = groupIt->next();
+ if ( !group ) continue;
+ groupDS = group->GetGroupDS();
+ if ( !groupDS || groupDS->IsEmpty() ) continue;
+ std::string grpName = group->GetName();
+ //MESSAGE("grpName=" << grpName);
+ if (grpName == groupName)
+ break;
+ else
+ groupDS = 0;
+ }
+
+ bool isNodeGroup = false;
+ bool isNodeCoords = false;
+ if (groupDS)
+ {
+ if (groupDS->GetType() != SMDSAbs_Node)
+ return;
+ isNodeGroup = true; // a group of nodes exists and it is in this mesh
+ }
+
+ if (nodesCoords.size() > 0)
+ isNodeCoords = true; // a list o nodes given by their coordinates
+ //MESSAGE("---" << isNodeGroup << " " << isNodeCoords);
+
+ // --- define groups to build
+
+ int idg; // --- group of SMDS volumes
+ string grpvName = groupName;
+ grpvName += "_vol";
+ SMESH_Group *grp = this->myMesh->AddGroup(SMDSAbs_Volume, grpvName.c_str(), idg);
+ if (!grp)
+ {
+ MESSAGE("group not created " << grpvName);
+ return;
+ }
+ SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(grp->GetGroupDS());
+
+ int idgs; // --- group of SMDS faces on the skin
+ string grpsName = groupName;
+ grpsName += "_skin";
+ SMESH_Group *grps = this->myMesh->AddGroup(SMDSAbs_Face, grpsName.c_str(), idgs);
+ if (!grps)
+ {
+ MESSAGE("group not created " << grpsName);
+ return;
+ }
+ SMESHDS_Group *sgrps = dynamic_cast<SMESHDS_Group*>(grps->GetGroupDS());
+
+ int idgi; // --- group of SMDS faces internal (several shapes)
+ string grpiName = groupName;
+ grpiName += "_internalFaces";
+ SMESH_Group *grpi = this->myMesh->AddGroup(SMDSAbs_Face, grpiName.c_str(), idgi);
+ if (!grpi)
+ {
+ MESSAGE("group not created " << grpiName);
+ return;
+ }
+ SMESHDS_Group *sgrpi = dynamic_cast<SMESHDS_Group*>(grpi->GetGroupDS());
+
+ int idgei; // --- group of SMDS faces internal (several shapes)
+ string grpeiName = groupName;
+ grpeiName += "_internalEdges";
+ SMESH_Group *grpei = this->myMesh->AddGroup(SMDSAbs_Edge, grpeiName.c_str(), idgei);
+ if (!grpei)
+ {
+ MESSAGE("group not created " << grpeiName);
+ return;
+ }
+ SMESHDS_Group *sgrpei = dynamic_cast<SMESHDS_Group*>(grpei->GetGroupDS());
+
+ // --- build downward connectivity
+
+ SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
+ meshDS->BuildDownWardConnectivity(true);
+ SMDS_UnstructuredGrid* grid = meshDS->getGrid();
+
+ // --- set of volumes detected inside
+
+ std::set<int> setOfInsideVol;
+ std::set<int> setOfVolToCheck;
+
+ std::vector<gp_Pnt> gpnts;
+ gpnts.clear();
+
+ if (isNodeGroup) // --- a group of nodes is provided : find all the volumes using one or more of this nodes
+ {
+ MESSAGE("group of nodes provided");
+ SMDS_ElemIteratorPtr elemIt = groupDS->GetElements();
+ while ( elemIt->more() )
+ {
+ const SMDS_MeshElement* elem = elemIt->next();
+ if (!elem)
+ continue;
+ const SMDS_MeshNode* node = dynamic_cast<const SMDS_MeshNode*>(elem);
+ if (!node)
+ continue;
+ SMDS_MeshElement* vol = 0;
+ SMDS_ElemIteratorPtr volItr = node->GetInverseElementIterator(SMDSAbs_Volume);
+ while (volItr->more())
+ {
+ vol = (SMDS_MeshElement*)volItr->next();
+ setOfInsideVol.insert(vol->getVtkId());
+ sgrp->Add(vol->GetID());
+ }
+ }
+ }
+ else if (isNodeCoords)
+ {
+ MESSAGE("list of nodes coordinates provided");
+ int i = 0;
+ int k = 0;
+ while (i < nodesCoords.size()-2)
+ {
+ double x = nodesCoords[i++];
+ double y = nodesCoords[i++];
+ double z = nodesCoords[i++];
+ gp_Pnt p = gp_Pnt(x, y ,z);
+ gpnts.push_back(p);
+ MESSAGE("TopoDS_Vertex " << k++ << " " << p.X() << " " << p.Y() << " " << p.Z());
+ }
+ }
+ else // --- no group, no coordinates : use the vertices of the geom shape provided, and radius
+ {
+ MESSAGE("no group of nodes provided, using vertices from geom shape, and radius");
+ TopTools_IndexedMapOfShape vertexMap;
+ TopExp::MapShapes( theShape, TopAbs_VERTEX, vertexMap );
+ gp_Pnt p = gp_Pnt(0,0,0);
+ if (vertexMap.Extent() < 1)
+ return;
+
+ for ( int i = 1; i <= vertexMap.Extent(); ++i )
+ {
+ const TopoDS_Vertex& vertex = TopoDS::Vertex( vertexMap( i ));
+ p = BRep_Tool::Pnt(vertex);
+ gpnts.push_back(p);
+ MESSAGE("TopoDS_Vertex " << i << " " << p.X() << " " << p.Y() << " " << p.Z());
+ }
+ }
+
+ if (gpnts.size() > 0)
+ {
+ int nodeId = 0;
+ const SMDS_MeshNode* startNode = theNodeSearcher->FindClosestTo(gpnts[0]);
+ if (startNode)
+ nodeId = startNode->GetID();
+ MESSAGE("nodeId " << nodeId);
+
+ double radius2 = radius*radius;
+ MESSAGE("radius2 " << radius2);
+
+ // --- volumes on start node
+
+ setOfVolToCheck.clear();
+ SMDS_MeshElement* startVol = 0;
+ SMDS_ElemIteratorPtr volItr = startNode->GetInverseElementIterator(SMDSAbs_Volume);
+ while (volItr->more())
+ {
+ startVol = (SMDS_MeshElement*)volItr->next();
+ setOfVolToCheck.insert(startVol->getVtkId());
+ }
+ if (setOfVolToCheck.empty())
+ {
+ MESSAGE("No volumes found");
+ return;
+ }
+
+ // --- starting with central volumes then their neighbors, check if they are inside
+ // or outside the domain, until no more new neighbor volume is inside.
+ // Fill the group of inside volumes
+
+ std::map<int, double> mapOfNodeDistance2;
+ mapOfNodeDistance2.clear();
+ std::set<int> setOfOutsideVol;
+ while (!setOfVolToCheck.empty())
+ {
+ std::set<int>::iterator it = setOfVolToCheck.begin();
+ int vtkId = *it;
+ MESSAGE("volume to check, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+ bool volInside = false;
+ vtkIdType npts = 0;
+ vtkIdType* pts = 0;
+ grid->GetCellPoints(vtkId, npts, pts);
+ for (int i=0; i<npts; i++)
+ {
+ double distance2 = 0;
+ if (mapOfNodeDistance2.count(pts[i]))
+ {
+ distance2 = mapOfNodeDistance2[pts[i]];
+ MESSAGE("point " << pts[i] << " distance2 " << distance2);
+ }
+ else
+ {
+ double *coords = grid->GetPoint(pts[i]);
+ gp_Pnt aPoint = gp_Pnt(coords[0], coords[1], coords[2]);
+ distance2 = 1.E40;
+ for (int j=0; j<gpnts.size(); j++)
+ {
+ double d2 = aPoint.SquareDistance(gpnts[j]);
+ if (d2 < distance2)
+ {
+ distance2 = d2;
+ if (distance2 < radius2)
+ break;
+ }
+ }
+ mapOfNodeDistance2[pts[i]] = distance2;
+ MESSAGE(" point " << pts[i] << " distance2 " << distance2 << " coords " << coords[0] << " " << coords[1] << " " << coords[2]);
+ }
+ if (distance2 < radius2)
+ {
+ volInside = true; // one or more nodes inside the domain
+ sgrp->Add(meshDS->fromVtkToSmds(vtkId));
+ break;
+ }
+ }
+ if (volInside)
+ {
+ setOfInsideVol.insert(vtkId);
+ MESSAGE(" volume inside, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+ int neighborsVtkIds[NBMAXNEIGHBORS];
+ int downIds[NBMAXNEIGHBORS];
+ unsigned char downTypes[NBMAXNEIGHBORS];
+ int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId);
+ for (int n = 0; n < nbNeighbors; n++)
+ if (!setOfInsideVol.count(neighborsVtkIds[n]) ||setOfOutsideVol.count(neighborsVtkIds[n]))
+ setOfVolToCheck.insert(neighborsVtkIds[n]);
+ }
+ else
+ {
+ setOfOutsideVol.insert(vtkId);
+ MESSAGE(" volume outside, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+ }
+ setOfVolToCheck.erase(vtkId);
+ }
+ }
+
+ // --- for outside hexahedrons, check if they have more than one neighbor volume inside
+ // If yes, add the volume to the inside set
+
+ bool addedInside = true;
+ std::set<int> setOfVolToReCheck;
+ while (addedInside)
+ {
+ MESSAGE(" --------------------------- re check");
+ addedInside = false;
+ std::set<int>::iterator itv = setOfInsideVol.begin();
+ for (; itv != setOfInsideVol.end(); ++itv)
+ {
+ int vtkId = *itv;
+ int neighborsVtkIds[NBMAXNEIGHBORS];
+ int downIds[NBMAXNEIGHBORS];
+ unsigned char downTypes[NBMAXNEIGHBORS];
+ int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId);
+ for (int n = 0; n < nbNeighbors; n++)
+ if (!setOfInsideVol.count(neighborsVtkIds[n]))
+ setOfVolToReCheck.insert(neighborsVtkIds[n]);
+ }
+ setOfVolToCheck = setOfVolToReCheck;
+ setOfVolToReCheck.clear();
+ while (!setOfVolToCheck.empty())
+ {
+ std::set<int>::iterator it = setOfVolToCheck.begin();
+ int vtkId = *it;
+ if (grid->GetCellType(vtkId) == VTK_HEXAHEDRON)
+ {
+ MESSAGE("volume to recheck, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+ int countInside = 0;
+ int neighborsVtkIds[NBMAXNEIGHBORS];
+ int downIds[NBMAXNEIGHBORS];
+ unsigned char downTypes[NBMAXNEIGHBORS];
+ int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId);
+ for (int n = 0; n < nbNeighbors; n++)
+ if (setOfInsideVol.count(neighborsVtkIds[n]))
+ countInside++;
+ MESSAGE("countInside " << countInside);
+ if (countInside > 1)
+ {
+ MESSAGE(" volume inside, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+ setOfInsideVol.insert(vtkId);
+ sgrp->Add(meshDS->fromVtkToSmds(vtkId));
+ addedInside = true;
+ }
+ else
+ setOfVolToReCheck.insert(vtkId);
+ }
+ setOfVolToCheck.erase(vtkId);
+ }
+ }
+
+ // --- map of Downward faces at the boundary, inside the global volume
+ // map of Downward faces on the skin of the global volume (equivalent to SMDS faces on the skin)
+ // fill group of SMDS faces inside the volume (when several volume shapes)
+ // fill group of SMDS faces on the skin of the global volume (if skin)
+
+ std::map<DownIdType, int, DownIdCompare> boundaryFaces; // boundary faces inside the volume --> corresponding cell
+ std::map<DownIdType, int, DownIdCompare> skinFaces; // faces on the skin of the global volume --> corresponding cell
+ std::set<int>::iterator it = setOfInsideVol.begin();
+ for (; it != setOfInsideVol.end(); ++it)
+ {
+ int vtkId = *it;
+ //MESSAGE(" vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+ int neighborsVtkIds[NBMAXNEIGHBORS];
+ int downIds[NBMAXNEIGHBORS];
+ unsigned char downTypes[NBMAXNEIGHBORS];
+ int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId, true);
+ for (int n = 0; n < nbNeighbors; n++)
+ {
+ int neighborDim = SMDS_Downward::getCellDimension(grid->GetCellType(neighborsVtkIds[n]));
+ if (neighborDim == 3)
+ {
+ if (! setOfInsideVol.count(neighborsVtkIds[n])) // neighbor volume is not inside : face is boundary
+ {
+ DownIdType face(downIds[n], downTypes[n]);
+ boundaryFaces[face] = vtkId;
+ }
+ // if the face between to volumes is in the mesh, get it (internal face between shapes)
+ int vtkFaceId = grid->getDownArray(downTypes[n])->getVtkCellId(downIds[n]);
+ if (vtkFaceId >= 0)
+ {
+ sgrpi->Add(meshDS->fromVtkToSmds(vtkFaceId));
+ // find also the smds edges on this face
+ int nbEdges = grid->getDownArray(downTypes[n])->getNumberOfDownCells(downIds[n]);
+ const int* dEdges = grid->getDownArray(downTypes[n])->getDownCells(downIds[n]);
+ const unsigned char* dTypes = grid->getDownArray(downTypes[n])->getDownTypes(downIds[n]);
+ for (int i = 0; i < nbEdges; i++)
+ {
+ int vtkEdgeId = grid->getDownArray(dTypes[i])->getVtkCellId(dEdges[i]);
+ if (vtkEdgeId >= 0)
+ sgrpei->Add(meshDS->fromVtkToSmds(vtkEdgeId));
+ }
+ }
+ }
+ else if (neighborDim == 2) // skin of the volume
+ {
+ DownIdType face(downIds[n], downTypes[n]);
+ skinFaces[face] = vtkId;
+ int vtkFaceId = grid->getDownArray(downTypes[n])->getVtkCellId(downIds[n]);
+ if (vtkFaceId >= 0)
+ sgrps->Add(meshDS->fromVtkToSmds(vtkFaceId));
+ }
+ }
+ }
+
+ // --- identify the edges constituting the wire of each subshape on the skin
+ // define polylines with the nodes of edges, equivalent to wires
+ // project polylines on subshapes, and partition, to get geom faces
+
+ std::map<int, std::set<int> > shapeIdToVtkIdSet; // shapeId --> set of vtkId on skin
+ std::set<int> emptySet;
+ emptySet.clear();
+ std::set<int> shapeIds;
+
+ SMDS_ElemIteratorPtr itelem = sgrps->GetElements();
+ while (itelem->more())
+ {
+ const SMDS_MeshElement *elem = itelem->next();
+ int shapeId = elem->getshapeId();
+ int vtkId = elem->getVtkId();
+ if (!shapeIdToVtkIdSet.count(shapeId))
+ {
+ shapeIdToVtkIdSet[shapeId] = emptySet;
+ shapeIds.insert(shapeId);
+ }
+ shapeIdToVtkIdSet[shapeId].insert(vtkId);
+ }
+
+ std::map<int, std::set<DownIdType, DownIdCompare> > shapeIdToEdges; // shapeId --> set of downward edges
+ std::set<DownIdType, DownIdCompare> emptyEdges;
+ emptyEdges.clear();
+
+ std::map<int, std::set<int> >::iterator itShape = shapeIdToVtkIdSet.begin();
+ for (; itShape != shapeIdToVtkIdSet.end(); ++itShape)
+ {
+ int shapeId = itShape->first;
+ MESSAGE(" --- Shape ID --- "<< shapeId);
+ shapeIdToEdges[shapeId] = emptyEdges;
+
+ std::vector<int> nodesEdges;
+
+ std::set<int>::iterator its = itShape->second.begin();
+ for (; its != itShape->second.end(); ++its)
+ {
+ int vtkId = *its;
+ MESSAGE(" " << vtkId);
+ int neighborsVtkIds[NBMAXNEIGHBORS];
+ int downIds[NBMAXNEIGHBORS];
+ unsigned char downTypes[NBMAXNEIGHBORS];
+ int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId);
+ for (int n = 0; n < nbNeighbors; n++)
+ {
+ if (neighborsVtkIds[n]<0) // only smds faces are considered as neighbors here
+ continue;
+ int smdsId = meshDS->fromVtkToSmds(neighborsVtkIds[n]);
+ const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
+ if ( shapeIds.count(elem->getshapeId()) && !sgrps->Contains(elem)) // edge : neighbor in the set of shape, not in the group
+ {
+ DownIdType edge(downIds[n], downTypes[n]);
+ if (!shapeIdToEdges[shapeId].count(edge))
+ {
+ shapeIdToEdges[shapeId].insert(edge);
+ int vtkNodeId[3];
+ int nbNodes = grid->getDownArray(downTypes[n])->getNodes(downIds[n],vtkNodeId);
+ nodesEdges.push_back(vtkNodeId[0]);
+ nodesEdges.push_back(vtkNodeId[nbNodes-1]);
+ MESSAGE(" --- nodes " << vtkNodeId[0]+1 << " " << vtkNodeId[nbNodes-1]+1);
+ }
+ }
+ }
+ }
+
+ std::list<int> order;
+ order.clear();
+ if (nodesEdges.size() > 0)
+ {
+ order.push_back(nodesEdges[0]); MESSAGE(" --- back " << order.back()+1); // SMDS id = VTK id + 1;
+ nodesEdges[0] = -1;
+ order.push_back(nodesEdges[1]); MESSAGE(" --- back " << order.back()+1);
+ nodesEdges[1] = -1; // do not reuse this edge
+ bool found = true;
+ while (found)
+ {
+ int nodeTofind = order.back(); // try first to push back
+ int i = 0;
+ for (i = 0; i<nodesEdges.size(); i++)
+ if (nodesEdges[i] == nodeTofind)
+ break;
+ if (i == nodesEdges.size())
+ found = false; // no follower found on back
+ else
+ {
+ if (i%2) // odd ==> use the previous one
+ if (nodesEdges[i-1] < 0)
+ found = false;
+ else
+ {
+ order.push_back(nodesEdges[i-1]); MESSAGE(" --- back " << order.back()+1);
+ nodesEdges[i-1] = -1;
+ }
+ else // even ==> use the next one
+ if (nodesEdges[i+1] < 0)
+ found = false;
+ else
+ {
+ order.push_back(nodesEdges[i+1]); MESSAGE(" --- back " << order.back()+1);
+ nodesEdges[i+1] = -1;
+ }
+ }
+ if (found)
+ continue;
+ // try to push front
+ found = true;
+ nodeTofind = order.front(); // try to push front
+ for (i = 0; i<nodesEdges.size(); i++)
+ if (nodesEdges[i] == nodeTofind)
+ break;
+ if (i == nodesEdges.size())
+ {
+ found = false; // no predecessor found on front
+ continue;
+ }
+ if (i%2) // odd ==> use the previous one
+ if (nodesEdges[i-1] < 0)
+ found = false;
+ else
+ {
+ order.push_front(nodesEdges[i-1]); MESSAGE(" --- front " << order.front()+1);
+ nodesEdges[i-1] = -1;
+ }
+ else // even ==> use the next one
+ if (nodesEdges[i+1] < 0)
+ found = false;
+ else
+ {
+ order.push_front(nodesEdges[i+1]); MESSAGE(" --- front " << order.front()+1);
+ nodesEdges[i+1] = -1;
+ }
+ }
+ }
+
+
+ std::vector<int> nodes;
+ nodes.push_back(shapeId);
+ std::list<int>::iterator itl = order.begin();
+ for (; itl != order.end(); itl++)
+ {
+ nodes.push_back((*itl) + 1); // SMDS id = VTK id + 1;
+ MESSAGE(" ordered node " << nodes[nodes.size()-1]);
+ }
+ listOfListOfNodes.push_back(nodes);
+ }
+
+ // partition geom faces with blocFissure
+ // mesh blocFissure and geom faces of the skin (external wires given, triangle algo to choose)
+ // mesh volume around blocFissure (skin triangles and quadrangle given, tetra algo to choose)
+
+ return;
+}
+
+
//================================================================================
/*!
* \brief Generates skin mesh (containing 2D cells) from 3D mesh
SMESH_MeshEditor tgtEditor2( tgtEditor.GetMesh() );
presentEditor = toAddExistingBondary ? &tgtEditor : &tgtEditor2;
+ SMESH_MesherHelper helper( *myMesh );
+ const TopAbs_ShapeEnum missShapeType = ( missType==SMDSAbs_Face ? TopAbs_FACE : TopAbs_EDGE );
SMDS_VolumeTool vTool;
TIDSortedElemSet avoidSet;
const TIDSortedElemSet emptySet, *elemSet = aroundElements ? &elements : &emptySet;
missType,
/*noMedium=*/false))
continue;
- tgtEditor.AddElement(nodes, missType, !iQuad && nodes.size()/(iQuad+1)>4);
+ SMDS_MeshElement* elem =
+ tgtEditor.AddElement(nodes, missType, !iQuad && nodes.size()/(iQuad+1)>4);
++nbAddedBnd;
+
+ // try to set a new element to a shape
+ if ( myMesh->HasShapeToMesh() )
+ {
+ bool ok = true;
+ set< pair<TopAbs_ShapeEnum, int > > mediumShapes;
+ const int nbN = nodes.size() / (iQuad+1 );
+ for ( inode = 0; inode < nbN && ok; ++inode )
+ {
+ pair<int, TopAbs_ShapeEnum> i_stype =
+ helper.GetMediumPos( nodes[inode], nodes[(inode+1)%nbN]);
+ if (( ok = ( i_stype.first > 0 && i_stype.second >= TopAbs_FACE )))
+ mediumShapes.insert( make_pair ( i_stype.second, i_stype.first ));
+ }
+ if ( ok && mediumShapes.size() > 1 )
+ {
+ set< pair<TopAbs_ShapeEnum, int > >::iterator stype_i = mediumShapes.begin();
+ pair<TopAbs_ShapeEnum, int> stype_i_0 = *stype_i;
+ for ( ++stype_i; stype_i != mediumShapes.end() && ok; ++stype_i )
+ {
+ if (( ok = ( stype_i->first != stype_i_0.first )))
+ ok = helper.IsSubShape( aMesh->IndexToShape( stype_i->second ),
+ aMesh->IndexToShape( stype_i_0.second ));
+ }
+ }
+ if ( ok && mediumShapes.begin()->first == missShapeType )
+ aMesh->SetMeshElementOnShape( elem, mediumShapes.begin()->second );
+ }
}
// ----------------------------------
{
presentEditor->myLastCreatedElems.Append(presentBndElems[i]);
}
-
+
} // loop on given elements
// ---------------------------------------------
