+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::RotationSweep(TIDSortedElemSet theElemSets[2],
+ const gp_Ax1& theAxis,
+ const double theAngle,
+ const int theNbSteps,
+ const double theTol,
+ const bool theMakeGroups,
+ const bool theMakeWalls)
+{
+ ClearLastCreated();
+
+ setElemsFirst( theElemSets );
+ myLastCreatedElems.reserve( theElemSets[0].size() * theNbSteps );
+ myLastCreatedNodes.reserve( theElemSets[1].size() * theNbSteps );
+
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+ srcElems.reserve( theElemSets[0].size() );
+ srcNodes.reserve( theElemSets[1].size() );
+
+ gp_Trsf aTrsf;
+ aTrsf.SetRotation( theAxis, theAngle );
+ gp_Trsf aTrsf2;
+ aTrsf2.SetRotation( theAxis, theAngle/2. );
+
+ gp_Lin aLine( theAxis );
+ double aSqTol = theTol * theTol;
+
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ TNodeOfNodeListMap mapNewNodes;
+ TElemOfVecOfNnlmiMap mapElemNewNodes;
+ TTElemOfElemListMap newElemsMap;
+
+ const bool isQuadraticMesh = bool( myMesh->NbEdges(ORDER_QUADRATIC) +
+ myMesh->NbFaces(ORDER_QUADRATIC) +
+ myMesh->NbVolumes(ORDER_QUADRATIC) );
+ // loop on theElemSets
+ TIDSortedElemSet::iterator itElem;
+ for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+ {
+ TIDSortedElemSet& theElems = theElemSets[ is2ndSet ];
+ for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
+ const SMDS_MeshElement* elem = *itElem;
+ if ( !elem || elem->GetType() == SMDSAbs_Volume )
+ continue;
+ vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
+ newNodesItVec.reserve( elem->NbNodes() );
+
+ // loop on elem nodes
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() )
+ {
+ const SMDS_MeshNode* node = cast2Node( itN->next() );
+
+ gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
+ double coord[3];
+ aXYZ.Coord( coord[0], coord[1], coord[2] );
+ bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
+
+ // check if a node has been already sweeped
+ TNodeOfNodeListMapItr nIt =
+ mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
+ list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
+ if ( listNewNodes.empty() )
+ {
+ // check if we are to create medium nodes between corner ones
+ bool needMediumNodes = false;
+ if ( isQuadraticMesh )
+ {
+ SMDS_ElemIteratorPtr it = node->GetInverseElementIterator();
+ while (it->more() && !needMediumNodes )
+ {
+ const SMDS_MeshElement* invElem = it->next();
+ if ( invElem != elem && !theElems.count( invElem )) continue;
+ needMediumNodes = ( invElem->IsQuadratic() && !invElem->IsMediumNode(node) );
+ if ( !needMediumNodes && invElem->GetEntityType() == SMDSEntity_BiQuad_Quadrangle )
+ needMediumNodes = true;
+ }
+ }
+
+ // make new nodes
+ const SMDS_MeshNode * newNode = node;
+ for ( int i = 0; i < theNbSteps; i++ ) {
+ if ( !isOnAxis ) {
+ if ( needMediumNodes ) // create a medium node
+ {
+ aTrsf2.Transforms( coord[0], coord[1], coord[2] );
+ newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+ myLastCreatedNodes.push_back(newNode);
+ srcNodes.push_back( node );
+ listNewNodes.push_back( newNode );
+ aTrsf2.Transforms( coord[0], coord[1], coord[2] );
+ }
+ else {
+ aTrsf.Transforms( coord[0], coord[1], coord[2] );
+ }
+ // create a corner node
+ newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+ myLastCreatedNodes.push_back(newNode);
+ srcNodes.push_back( node );
+ listNewNodes.push_back( newNode );
+ }
+ else {
+ listNewNodes.push_back( newNode );
+ // if ( needMediumNodes )
+ // listNewNodes.push_back( newNode );
+ }
+ }
+ }
+ newNodesItVec.push_back( nIt );
+ }
+ // make new elements
+ sweepElement( elem, newNodesItVec, newElemsMap[elem], theNbSteps, srcElems );
+ }
+ }
+
+ if ( theMakeWalls )
+ makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElemSets[0], theNbSteps, srcElems );
+
+ PGroupIDs newGroupIDs;
+ if ( theMakeGroups )
+ newGroupIDs = generateGroups( srcNodes, srcElems, "rotated");
+
+ return newGroupIDs;
+}
+
+//=======================================================================
+//function : ExtrusParam
+//purpose : standard construction
+//=======================================================================
+
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Vec& theStep,
+ const int theNbSteps,
+ const std::list<double>& theScales,
+ const std::list<double>& theAngles,
+ const gp_XYZ* theBasePoint,
+ const int theFlags,
+ const double theTolerance):
+ myDir( theStep ),
+ myBaseP( Precision::Infinite(), 0, 0 ),
+ myFlags( theFlags ),
+ myTolerance( theTolerance ),
+ myElemsToUse( NULL )
+{
+ mySteps = new TColStd_HSequenceOfReal;
+ const double stepSize = theStep.Magnitude();
+ for (int i=1; i<=theNbSteps; i++ )
+ mySteps->Append( stepSize );
+
+ if ( !theScales.empty() )
+ {
+ if ( IsScaleVariation() && (int)theScales.size() < theNbSteps )
+ linearScaleVariation( theNbSteps, const_cast< std::list<double>& >( theScales ));
+
+ // add medium scales
+ std::list<double>::const_iterator s2 = theScales.begin(), s1 = s2++;
+ myScales.reserve( theNbSteps * 2 );
+ myScales.push_back( 0.5 * ( *s1 + 1. ));
+ myScales.push_back( *s1 );
+ for ( ; s2 != theScales.end(); s1 = s2++ )
+ {
+ myScales.push_back( 0.5 * ( *s1 + *s2 ));
+ myScales.push_back( *s2 );
+ }
+ }
+
+ if ( !theAngles.empty() )
+ {
+ std::list<double>& angles = const_cast< std::list<double>& >( theAngles );
+ if ( IsAngleVariation() && (int)theAngles.size() < theNbSteps )
+ linearAngleVariation( theNbSteps, angles );
+
+ // accumulate angles
+ double angle = 0;
+ int nbAngles = 0;
+ std::list<double>::iterator a1 = angles.begin(), a2;
+ for ( ; a1 != angles.end(); ++a1, ++nbAngles )
+ {
+ angle += *a1;
+ *a1 = angle;
+ }
+ while ( nbAngles++ < theNbSteps )
+ angles.push_back( angles.back() );
+
+ // add medium angles
+ a2 = angles.begin(), a1 = a2++;
+ myAngles.push_back( 0.5 * *a1 );
+ myAngles.push_back( *a1 );
+ for ( ; a2 != angles.end(); a1 = a2++ )
+ {
+ myAngles.push_back( 0.5 * ( *a1 + *a2 ));
+ myAngles.push_back( *a2 );
+ }
+ }
+
+ if ( theBasePoint )
+ {
+ myBaseP = *theBasePoint;
+ }
+
+ if (( theFlags & EXTRUSION_FLAG_SEW ) &&
+ ( theTolerance > 0 ))
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDirAndSew;
+ }
+ else
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDir;
+ }
+}
+
+//=======================================================================
+//function : ExtrusParam
+//purpose : steps are given explicitly
+//=======================================================================
+
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Dir& theDir,
+ Handle(TColStd_HSequenceOfReal) theSteps,
+ const int theFlags,
+ const double theTolerance):
+ myDir( theDir ),
+ mySteps( theSteps ),
+ myFlags( theFlags ),
+ myTolerance( theTolerance ),
+ myElemsToUse( NULL )
+{
+ if (( theFlags & EXTRUSION_FLAG_SEW ) &&
+ ( theTolerance > 0 ))
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDirAndSew;
+ }
+ else
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDir;
+ }
+}
+
+//=======================================================================
+//function : ExtrusParam
+//purpose : for extrusion by normal
+//=======================================================================
+
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const double theStepSize,
+ const int theNbSteps,
+ const int theFlags,
+ const int theDim ):
+ myDir( 1,0,0 ),
+ mySteps( new TColStd_HSequenceOfReal ),
+ myFlags( theFlags ),
+ myTolerance( 0 ),
+ myElemsToUse( NULL )
+{
+ for (int i = 0; i < theNbSteps; i++ )
+ mySteps->Append( theStepSize );
+
+ if ( theDim == 1 )
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByNormal1D;
+ }
+ else
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByNormal2D;
+ }
+}
+
+//=======================================================================
+//function : ExtrusParam
+//purpose : for extrusion along path
+//=======================================================================
+
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const std::vector< PathPoint >& thePoints,
+ const gp_Pnt* theBasePoint,
+ const std::list<double>& theScales,
+ const bool theMakeGroups )
+ : myBaseP( Precision::Infinite(), 0, 0 ),
+ myFlags( EXTRUSION_FLAG_BOUNDARY | ( theMakeGroups ? EXTRUSION_FLAG_GROUPS : 0 )),
+ myPathPoints( thePoints )
+{
+ if ( theBasePoint )
+ {
+ myBaseP = theBasePoint->XYZ();
+ }
+
+ if ( !theScales.empty() )
+ {
+ // add medium scales
+ std::list<double>::const_iterator s2 = theScales.begin(), s1 = s2++;
+ myScales.reserve( thePoints.size() * 2 );
+ myScales.push_back( 0.5 * ( 1. + *s1 ));
+ myScales.push_back( *s1 );
+ for ( ; s2 != theScales.end(); s1 = s2++ )
+ {
+ myScales.push_back( 0.5 * ( *s1 + *s2 ));
+ myScales.push_back( *s2 );
+ }
+ }
+
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesAlongTrack;
+}
+
+//=======================================================================
+//function : ExtrusParam::SetElementsToUse
+//purpose : stores elements to use for extrusion by normal, depending on
+// state of EXTRUSION_FLAG_USE_INPUT_ELEMS_ONLY flag;
+// define myBaseP for scaling
+//=======================================================================
+
+void SMESH_MeshEditor::ExtrusParam::SetElementsToUse( const TIDSortedElemSet& elems,
+ const TIDSortedElemSet& nodes )
+{
+ myElemsToUse = ToUseInpElemsOnly() ? & elems : 0;
+
+ if ( Precision::IsInfinite( myBaseP.X() )) // myBaseP not defined
+ {
+ myBaseP.SetCoord( 0.,0.,0. );
+ TIDSortedElemSet newNodes;
+
+ const TIDSortedElemSet* elemSets[] = { &elems, &nodes };
+ for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+ {
+ const TIDSortedElemSet& elements = *( elemSets[ is2ndSet ]);
+ TIDSortedElemSet::const_iterator itElem = elements.begin();
+ for ( ; itElem != elements.end(); itElem++ )
+ {
+ const SMDS_MeshElement* elem = *itElem;
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() ) {
+ const SMDS_MeshElement* node = itN->next();
+ if ( newNodes.insert( node ).second )
+ myBaseP += SMESH_NodeXYZ( node );
+ }
+ }
+ }
+ myBaseP /= newNodes.size();
+ }
+}
+
+//=======================================================================
+//function : ExtrusParam::beginStepIter
+//purpose : prepare iteration on steps
+//=======================================================================
+
+void SMESH_MeshEditor::ExtrusParam::beginStepIter( bool withMediumNodes )
+{
+ myWithMediumNodes = withMediumNodes;
+ myNextStep = 1;
+ myCurSteps.clear();
+}
+//=======================================================================
+//function : ExtrusParam::moreSteps
+//purpose : are there more steps?
+//=======================================================================
+
+bool SMESH_MeshEditor::ExtrusParam::moreSteps()
+{
+ return myNextStep <= mySteps->Length() || !myCurSteps.empty();
+}
+//=======================================================================
+//function : ExtrusParam::nextStep
+//purpose : returns the next step
+//=======================================================================
+
+double SMESH_MeshEditor::ExtrusParam::nextStep()
+{
+ double res = 0;
+ if ( !myCurSteps.empty() )
+ {
+ res = myCurSteps.back();
+ myCurSteps.pop_back();
+ }
+ else if ( myNextStep <= mySteps->Length() )
+ {
+ myCurSteps.push_back( mySteps->Value( myNextStep ));
+ ++myNextStep;
+ if ( myWithMediumNodes )
+ {
+ myCurSteps.back() /= 2.;
+ myCurSteps.push_back( myCurSteps.back() );
+ }
+ res = nextStep();
+ }
+ return res;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesByDir
+//purpose : create nodes for standard extrusion
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByDir( SMESHDS_Mesh* mesh,
+ const SMDS_MeshNode* srcNode,
+ std::list<const SMDS_MeshNode*> & newNodes,
+ const bool makeMediumNodes)
+{
+ gp_XYZ p = SMESH_NodeXYZ( srcNode );
+
+ int nbNodes = 0;
+ for ( beginStepIter( makeMediumNodes ); moreSteps(); ++nbNodes ) // loop on steps
+ {
+ p += myDir.XYZ() * nextStep();
+ const SMDS_MeshNode * newNode = mesh->AddNode( p.X(), p.Y(), p.Z() );
+ newNodes.push_back( newNode );
+ }
+
+ if ( !myScales.empty() || !myAngles.empty() )
+ {
+ gp_XYZ center = myBaseP;
+ gp_Ax1 ratationAxis( center, myDir );
+ gp_Trsf rotation;
+
+ std::list<const SMDS_MeshNode*>::iterator nIt = newNodes.begin();
+ size_t i = !makeMediumNodes;
+ for ( beginStepIter( makeMediumNodes );
+ moreSteps();
+ ++nIt, i += 1 + !makeMediumNodes )
+ {
+ center += myDir.XYZ() * nextStep();
+
+ gp_XYZ xyz = SMESH_NodeXYZ( *nIt );
+ bool moved = false;
+ if ( i < myScales.size() )
+ {
+ xyz = ( myScales[i] * ( xyz - center )) + center;
+ moved = true;
+ }
+ if ( !myAngles.empty() )
+ {
+ rotation.SetRotation( ratationAxis, myAngles[i] );
+ rotation.Transforms( xyz );
+ moved = true;
+ }
+ if ( moved )
+ mesh->MoveNode( *nIt, xyz.X(), xyz.Y(), xyz.Z() );
+ else
+ break;
+ }
+ }
+ return nbNodes;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesByDirAndSew
+//purpose : create nodes for standard extrusion with sewing
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByDirAndSew( SMESHDS_Mesh* mesh,
+ const SMDS_MeshNode* srcNode,
+ std::list<const SMDS_MeshNode*> & newNodes,
+ const bool makeMediumNodes)
+{
+ gp_XYZ P1 = SMESH_NodeXYZ( srcNode );
+
+ int nbNodes = 0;
+ for ( beginStepIter( makeMediumNodes ); moreSteps(); ++nbNodes ) // loop on steps
+ {
+ P1 += myDir.XYZ() * nextStep();
+
+ // try to search in sequence of existing nodes
+ // if myNodes.size()>0 we 'nave to use given sequence
+ // else - use all nodes of mesh
+ const SMDS_MeshNode * node = 0;
+ if ( myNodes.Length() > 0 )
+ {
+ for ( int i = 1; i <= myNodes.Length(); i++ )
+ {
+ SMESH_NodeXYZ P2 = myNodes.Value(i);
+ if (( P1 - P2 ).SquareModulus() < myTolerance * myTolerance )
+ {
+ node = myNodes.Value(i);
+ break;
+ }
+ }
+ }
+ else
+ {
+ SMDS_NodeIteratorPtr itn = mesh->nodesIterator();
+ while(itn->more())
+ {
+ SMESH_NodeXYZ P2 = itn->next();
+ if (( P1 - P2 ).SquareModulus() < myTolerance * myTolerance )
+ {
+ node = P2._node;
+ break;
+ }
+ }
+ }
+
+ if ( !node )
+ node = mesh->AddNode( P1.X(), P1.Y(), P1.Z() );
+
+ newNodes.push_back( node );
+
+ } // loop on steps
+
+ return nbNodes;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesByNormal2D
+//purpose : create nodes for extrusion using normals of faces
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByNormal2D( SMESHDS_Mesh* mesh,
+ const SMDS_MeshNode* srcNode,
+ std::list<const SMDS_MeshNode*> & newNodes,
+ const bool makeMediumNodes)
+{
+ const bool alongAvgNorm = ( myFlags & EXTRUSION_FLAG_BY_AVG_NORMAL );
+
+ gp_XYZ p = SMESH_NodeXYZ( srcNode );
+
+ // get normals to faces sharing srcNode
+ vector< gp_XYZ > norms, baryCenters;
+ gp_XYZ norm, avgNorm( 0,0,0 );
+ SMDS_ElemIteratorPtr faceIt = srcNode->GetInverseElementIterator( SMDSAbs_Face );
+ while ( faceIt->more() )
+ {
+ const SMDS_MeshElement* face = faceIt->next();
+ if ( myElemsToUse && !myElemsToUse->count( face ))
+ continue;
+ if ( SMESH_MeshAlgos::FaceNormal( face, norm, /*normalized=*/true ))
+ {
+ norms.push_back( norm );
+ avgNorm += norm;
+ if ( !alongAvgNorm )
+ {
+ gp_XYZ bc(0,0,0);
+ int nbN = 0;
+ for ( SMDS_ElemIteratorPtr nIt = face->nodesIterator(); nIt->more(); ++nbN )
+ bc += SMESH_NodeXYZ( nIt->next() );
+ baryCenters.push_back( bc / nbN );
+ }
+ }
+ }
+
+ if ( norms.empty() ) return 0;
+
+ double normSize = avgNorm.Modulus();
+ if ( normSize < std::numeric_limits<double>::min() )
+ return 0;
+
+ if ( myFlags & EXTRUSION_FLAG_BY_AVG_NORMAL ) // extrude along avgNorm
+ {
+ myDir = avgNorm;
+ return makeNodesByDir( mesh, srcNode, newNodes, makeMediumNodes );
+ }
+
+ avgNorm /= normSize;
+
+ int nbNodes = 0;
+ for ( beginStepIter( makeMediumNodes ); moreSteps(); ++nbNodes ) // loop on steps
+ {
+ gp_XYZ pNew = p;
+ double stepSize = nextStep();
+
+ if ( norms.size() > 1 )
+ {
+ for ( size_t iF = 0; iF < norms.size(); ++iF ) // loop on faces
+ {
+ // translate plane of a face
+ baryCenters[ iF ] += norms[ iF ] * stepSize;
+
+ // find point of intersection of the face plane located at baryCenters[ iF ]
+ // and avgNorm located at pNew
+ double d = -( norms[ iF ] * baryCenters[ iF ]); // d of plane equation ax+by+cz+d=0
+ double dot = ( norms[ iF ] * avgNorm );
+ if ( dot < std::numeric_limits<double>::min() )
+ dot = stepSize * 1e-3;
+ double step = -( norms[ iF ] * pNew + d ) / dot;
+ pNew += step * avgNorm;
+ }
+ }
+ else
+ {
+ pNew += stepSize * avgNorm;
+ }
+ p = pNew;
+
+ const SMDS_MeshNode * newNode = mesh->AddNode( p.X(), p.Y(), p.Z() );
+ newNodes.push_back( newNode );
+ }
+ return nbNodes;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesByNormal1D
+//purpose : create nodes for extrusion using normals of edges
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByNormal1D( SMESHDS_Mesh* mesh,
+ const SMDS_MeshNode* srcNode,
+ std::list<const SMDS_MeshNode*> & newNodes,
+ const bool makeMediumNodes)
+{
+ throw SALOME_Exception("Extrusion 1D by Normal not implemented");
+ return 0;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesAlongTrack
+//purpose : create nodes for extrusion along path
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesAlongTrack( SMESHDS_Mesh* mesh,
+ const SMDS_MeshNode* srcNode,
+ std::list<const SMDS_MeshNode*> & newNodes,
+ const bool makeMediumNodes)
+{
+ const Standard_Real aTolAng=1.e-4;
+
+ gp_Pnt aV0x = myBaseP;
+ gp_Pnt aPN0 = SMESH_NodeXYZ( srcNode );
+
+ const PathPoint& aPP0 = myPathPoints[0];
+ gp_Pnt aP0x = aPP0.myPnt;
+ gp_Dir aDT0x= aPP0.myTgt;
+
+ std::vector< gp_Pnt > centers;
+ centers.reserve( NbSteps() * 2 );
+
+ gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
+
+ for ( size_t j = 1; j < myPathPoints.size(); ++j )
+ {
+ const PathPoint& aPP = myPathPoints[j];
+ const gp_Pnt& aP1x = aPP.myPnt;
+ const gp_Dir& aDT1x = aPP.myTgt;
+
+ // Translation
+ gp_Vec aV01x( aP0x, aP1x );
+ aTrsf.SetTranslation( aV01x );
+ gp_Pnt aV1x = aV0x.Transformed( aTrsf );
+ gp_Pnt aPN1 = aPN0.Transformed( aTrsf );
+
+ // rotation 1 [ T1,T0 ]
+ Standard_Real aAngleT1T0 = -aDT1x.Angle( aDT0x );
+ if ( fabs( aAngleT1T0 ) > aTolAng )
+ {
+ gp_Dir aDT1T0 = aDT1x ^ aDT0x;
+ aTrsfRotT1T0.SetRotation( gp_Ax1( aV1x, aDT1T0 ), aAngleT1T0 );
+
+ aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
+ }
+
+ // rotation 2
+ if ( aPP.myAngle != 0. )
+ {
+ aTrsfRot.SetRotation( gp_Ax1( aV1x, aDT1x ), aPP.myAngle );
+ aPN1 = aPN1.Transformed( aTrsfRot );
+ }
+
+ // make new node
+ if ( makeMediumNodes )
+ {
+ // create additional node
+ gp_XYZ midP = 0.5 * ( aPN1.XYZ() + aPN0.XYZ() );
+ const SMDS_MeshNode* newNode = mesh->AddNode( midP.X(), midP.Y(), midP.Z() );
+ newNodes.push_back( newNode );
+
+ }
+ const SMDS_MeshNode* newNode = mesh->AddNode( aPN1.X(), aPN1.Y(), aPN1.Z() );
+ newNodes.push_back( newNode );
+
+ centers.push_back( 0.5 * ( aV0x.XYZ() + aV1x.XYZ() ));
+ centers.push_back( aV1x );
+
+ aPN0 = aPN1;
+ aP0x = aP1x;
+ aV0x = aV1x;
+ aDT0x = aDT1x;
+ }
+
+ // scale
+ if ( !myScales.empty() )
+ {
+ gp_Trsf aTrsfScale;
+ std::list<const SMDS_MeshNode*>::iterator node = newNodes.begin();
+ for ( size_t i = !makeMediumNodes;
+ i < myScales.size() && node != newNodes.end();
+ i += ( 1 + !makeMediumNodes ), ++node )
+ {
+ aTrsfScale.SetScale( centers[ i ], myScales[ i ] );
+ gp_Pnt aN = SMESH_NodeXYZ( *node );
+ gp_Pnt aP = aN.Transformed( aTrsfScale );
+ mesh->MoveNode( *node, aP.X(), aP.Y(), aP.Z() );
+ }
+ }
+
+ return myPathPoints.size() + makeMediumNodes * ( myPathPoints.size() - 2 );
+}
+
+//=======================================================================
+//function : ExtrusionSweep
+//purpose :
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet theElems[2],
+ const gp_Vec& theStep,
+ const int theNbSteps,
+ TTElemOfElemListMap& newElemsMap,
+ const int theFlags,
+ const double theTolerance)
+{
+ std::list<double> dummy;
+ ExtrusParam aParams( theStep, theNbSteps, dummy, dummy, 0,
+ theFlags, theTolerance );
+ return ExtrusionSweep( theElems, aParams, newElemsMap );
+}
+
+namespace
+{
+
+//=======================================================================
+//function : getOriFactor
+//purpose : Return -1 or 1 depending on if order of given nodes corresponds to
+// edge curve orientation
+//=======================================================================
+
+ double getOriFactor( const TopoDS_Edge& edge,
+ const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ SMESH_MesherHelper& helper)
+ {
+ double u1 = helper.GetNodeU( edge, n1, n2 );
+ double u2 = helper.GetNodeU( edge, n2, n1 );
+ return u1 < u2 ? 1. : -1.;
+ }
+}
+
+//=======================================================================
+//function : ExtrusionSweep
+//purpose :
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet theElemSets[2],
+ ExtrusParam& theParams,
+ TTElemOfElemListMap& newElemsMap)
+{
+ ClearLastCreated();
+
+ setElemsFirst( theElemSets );
+ myLastCreatedElems.reserve( theElemSets[0].size() * theParams.NbSteps() );
+ myLastCreatedNodes.reserve( theElemSets[1].size() * theParams.NbSteps() );
+
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+ srcElems.reserve( theElemSets[0].size() );
+ srcNodes.reserve( theElemSets[1].size() );
+
+ const int nbSteps = theParams.NbSteps();
+ theParams.SetElementsToUse( theElemSets[0], theElemSets[1] );
+
+ TNodeOfNodeListMap mapNewNodes;
+ TElemOfVecOfNnlmiMap mapElemNewNodes;
+
+ const bool isQuadraticMesh = bool( myMesh->NbEdges(ORDER_QUADRATIC) +
+ myMesh->NbFaces(ORDER_QUADRATIC) +
+ myMesh->NbVolumes(ORDER_QUADRATIC) );
+ // loop on theElems
+ TIDSortedElemSet::iterator itElem;
+ for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+ {
+ TIDSortedElemSet& theElems = theElemSets[ is2ndSet ];
+ for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
+ {
+ // check element type
+ const SMDS_MeshElement* elem = *itElem;
+ if ( !elem || elem->GetType() == SMDSAbs_Volume )
+ continue;
+
+ const size_t nbNodes = elem->NbNodes();
+ vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
+ newNodesItVec.reserve( nbNodes );
+
+ // loop on elem nodes
+ SMDS_NodeIteratorPtr itN = elem->nodeIterator();
+ while ( itN->more() )
+ {
+ // check if a node has been already sweeped
+ const SMDS_MeshNode* node = itN->next();
+ TNodeOfNodeListMap::iterator nIt =
+ mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
+ list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
+ if ( listNewNodes.empty() )
+ {
+ // make new nodes
+
+ // check if we are to create medium nodes between corner ones
+ bool needMediumNodes = false;
+ if ( isQuadraticMesh )
+ {
+ SMDS_ElemIteratorPtr it = node->GetInverseElementIterator();
+ while (it->more() && !needMediumNodes )
+ {
+ const SMDS_MeshElement* invElem = it->next();
+ if ( invElem != elem && !theElems.count( invElem )) continue;
+ needMediumNodes = ( invElem->IsQuadratic() && !invElem->IsMediumNode(node) );
+ if ( !needMediumNodes && invElem->GetEntityType() == SMDSEntity_BiQuad_Quadrangle )
+ needMediumNodes = true;
+ }
+ }
+ // create nodes for all steps
+ if ( theParams.MakeNodes( GetMeshDS(), node, listNewNodes, needMediumNodes ))
+ {
+ list<const SMDS_MeshNode*>::iterator newNodesIt = listNewNodes.begin();
+ for ( ; newNodesIt != listNewNodes.end(); ++newNodesIt )
+ {
+ myLastCreatedNodes.push_back( *newNodesIt );
+ srcNodes.push_back( node );
+ }
+ }
+ else
+ {
+ if ( theParams.ToMakeBoundary() )
+ {
+ GetMeshDS()->Modified();
+ throw SALOME_Exception( SMESH_Comment("Can't extrude node #") << node->GetID() );
+ }
+ break; // newNodesItVec will be shorter than nbNodes
+ }
+ }
+ newNodesItVec.push_back( nIt );
+ }
+ // make new elements
+ if ( newNodesItVec.size() == nbNodes )
+ sweepElement( elem, newNodesItVec, newElemsMap[elem], nbSteps, srcElems );
+ }
+ }
+
+ if ( theParams.ToMakeBoundary() ) {
+ makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElemSets[0], nbSteps, srcElems );
+ }
+ PGroupIDs newGroupIDs;
+ if ( theParams.ToMakeGroups() )
+ newGroupIDs = generateGroups( srcNodes, srcElems, "extruded");
+
+ return newGroupIDs;
+}
+
+//=======================================================================
+//function : ExtrusionAlongTrack
+//purpose :
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet theElements[2],
+ SMESH_Mesh* theTrackMesh,
+ SMDS_ElemIteratorPtr theTrackIterator,
+ const SMDS_MeshNode* theN1,
+ std::list<double>& theAngles,
+ const bool theAngleVariation,
+ std::list<double>& theScales,
+ const bool theScaleVariation,
+ const gp_Pnt* theRefPoint,
+ const bool theMakeGroups)
+{
+ ClearLastCreated();
+
+ // 1. Check data
+ if ( theElements[0].empty() && theElements[1].empty() )
+ return EXTR_NO_ELEMENTS;
+
+ ASSERT( theTrackMesh );
+ if ( ! theTrackIterator || !theTrackIterator->more() )
+ return EXTR_NO_ELEMENTS;
+
+ // 2. Get ordered nodes
+ SMESH_MeshAlgos::TElemGroupVector branchEdges;
+ SMESH_MeshAlgos::TNodeGroupVector branchNods;
+ SMESH_MeshAlgos::Get1DBranches( theTrackIterator, branchEdges, branchNods, theN1 );
+ if ( branchEdges.empty() )
+ return EXTR_PATH_NOT_EDGE;
+
+ if ( branchEdges.size() > 1 )
+ return EXTR_BAD_PATH_SHAPE;
+
+ std::vector< const SMDS_MeshNode* >& pathNodes = branchNods[0];
+ std::vector< const SMDS_MeshElement* >& pathEdges = branchEdges[0];
+ if ( pathNodes[0] != theN1 && pathNodes[1] != theN1 )
+ return EXTR_BAD_STARTING_NODE;
+
+ if ( theTrackMesh->NbEdges( ORDER_QUADRATIC ) > 0 )
+ {
+ // add medium nodes to pathNodes
+ std::vector< const SMDS_MeshNode* > pathNodes2;
+ std::vector< const SMDS_MeshElement* > pathEdges2;
+ pathNodes2.reserve( pathNodes.size() * 2 );
+ pathEdges2.reserve( pathEdges.size() * 2 );
+ for ( size_t i = 0; i < pathEdges.size(); ++i )
+ {
+ pathNodes2.push_back( pathNodes[i] );
+ pathEdges2.push_back( pathEdges[i] );
+ if ( pathEdges[i]->IsQuadratic() )
+ {
+ pathNodes2.push_back( pathEdges[i]->GetNode(2) );
+ pathEdges2.push_back( pathEdges[i] );
+ }
+ }
+ pathNodes2.push_back( pathNodes.back() );
+ pathEdges.swap( pathEdges2 );
+ pathNodes.swap( pathNodes2 );
+ }
+
+ // 3. Get path data at pathNodes
+
+ std::vector< ExtrusParam::PathPoint > points( pathNodes.size() );
+
+ if ( theAngleVariation )
+ linearAngleVariation( points.size()-1, theAngles );
+ if ( theScaleVariation )
+ linearScaleVariation( points.size()-1, theScales );
+
+ theAngles.push_front( 0 ); // for the 1st point that is not transformed
+ std::list<double>::iterator angle = theAngles.begin();
+
+ SMESHDS_Mesh* pathMeshDS = theTrackMesh->GetMeshDS();
+
+ std::map< int, double > edgeID2OriFactor; // orientation of EDGEs
+ std::map< int, double >::iterator id2factor;
+ SMESH_MesherHelper pathHelper( *theTrackMesh );
+ gp_Pnt p; gp_Vec tangent;
+ const double tol2 = gp::Resolution() * gp::Resolution();
+
+ for ( size_t i = 0; i < pathNodes.size(); ++i )
+ {
+ ExtrusParam::PathPoint & point = points[ i ];
+
+ point.myPnt = SMESH_NodeXYZ( pathNodes[ i ]);
+
+ if ( angle != theAngles.end() )
+ point.myAngle = *angle++;
+
+ tangent.SetCoord( 0,0,0 );
+ const int shapeID = pathNodes[ i ]->GetShapeID();
+ const TopoDS_Shape& shape = pathMeshDS->IndexToShape( shapeID );
+ TopAbs_ShapeEnum shapeType = shape.IsNull() ? TopAbs_SHAPE : shape.ShapeType();
+ switch ( shapeType )
+ {
+ case TopAbs_EDGE:
+ {
+ TopoDS_Edge edge = TopoDS::Edge( shape );
+ id2factor = edgeID2OriFactor.insert( std::make_pair( shapeID, 0 )).first;
+ if ( id2factor->second == 0 )
+ {
+ if ( i ) id2factor->second = getOriFactor( edge, pathNodes[i-1], pathNodes[i], pathHelper );
+ else id2factor->second = getOriFactor( edge, pathNodes[i], pathNodes[i+1], pathHelper );
+ }
+ double u = pathHelper.GetNodeU( edge, pathNodes[i] ), u0, u1;
+ Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, u0, u1 );
+ curve->D1( u, p, tangent );
+ tangent *= id2factor->second;
+ break;
+ }
+ case TopAbs_VERTEX:
+ {
+ int nbEdges = 0;
+ PShapeIteratorPtr shapeIt = pathHelper.GetAncestors( shape, *theTrackMesh, TopAbs_EDGE );
+ while ( const TopoDS_Shape* edgePtr = shapeIt->next() )
+ {
+ int edgeID = pathMeshDS->ShapeToIndex( *edgePtr );
+ for ( int di = -1; di <= 0; ++di )
+ {
+ size_t j = i + di;
+ if ( j < pathEdges.size() && edgeID == pathEdges[ j ]->GetShapeID() )
+ {
+ TopoDS_Edge edge = TopoDS::Edge( *edgePtr );
+ id2factor = edgeID2OriFactor.insert( std::make_pair( edgeID, 0 )).first;
+ if ( id2factor->second == 0 )
+ {
+ if ( j < i )
+ id2factor->second = getOriFactor( edge, pathNodes[i-1], pathNodes[i], pathHelper );
+ else
+ id2factor->second = getOriFactor( edge, pathNodes[i], pathNodes[i+1], pathHelper );
+ }
+ double u = pathHelper.GetNodeU( edge, pathNodes[i] ), u0, u1;
+ Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, u0, u1 );
+ gp_Vec du;
+ curve->D1( u, p, du );
+ double size2 = du.SquareMagnitude();
+ if ( du.SquareMagnitude() > tol2 )
+ {
+ tangent += du.Divided( Sqrt( size2 )) * id2factor->second;
+ nbEdges++;
+ }
+ break;
+ }
+ }
+ }
+ if ( nbEdges > 0 )
+ break;
+ }
+ default:
+ {
+ for ( int di = -1; di <= 1; di += 2 )
+ {
+ size_t j = i + di;
+ if ( j < pathNodes.size() )
+ {
+ gp_Vec dir( point.myPnt, SMESH_NodeXYZ( pathNodes[ j ]));
+ double size2 = dir.SquareMagnitude();
+ if ( size2 > tol2 )
+ tangent += dir.Divided( Sqrt( size2 )) * di;
+ }
+ }
+ }
+ } // switch ( shapeType )
+
+ if ( tangent.SquareMagnitude() < tol2 )
+ return EXTR_CANT_GET_TANGENT;
+
+ point.myTgt = tangent;
+
+ } // loop on pathNodes
+
+
+ ExtrusParam nodeMaker( points, theRefPoint, theScales, theMakeGroups );
+ TTElemOfElemListMap newElemsMap;
+
+ ExtrusionSweep( theElements, nodeMaker, newElemsMap );
+
+ return EXTR_OK;
+}
+
+//=======================================================================
+//function : linearAngleVariation
+//purpose : spread values over nbSteps
+//=======================================================================
+
+void SMESH_MeshEditor::linearAngleVariation(const int nbSteps,
+ list<double>& Angles)
+{
+ int nbAngles = Angles.size();
+ if( nbSteps > nbAngles && nbAngles > 0 )
+ {
+ vector<double> theAngles(nbAngles);
+ theAngles.assign( Angles.begin(), Angles.end() );
+
+ list<double> res;
+ double rAn2St = double( nbAngles ) / double( nbSteps );
+ double angPrev = 0, angle;
+ for ( int iSt = 0; iSt < nbSteps; ++iSt )
+ {
+ double angCur = rAn2St * ( iSt+1 );
+ double angCurFloor = floor( angCur );
+ double angPrevFloor = floor( angPrev );
+ if ( angPrevFloor == angCurFloor )
+ angle = rAn2St * theAngles[ int( angCurFloor ) ];
+ else {
+ int iP = int( angPrevFloor );
+ double angPrevCeil = ceil(angPrev);
+ angle = ( angPrevCeil - angPrev ) * theAngles[ iP ];
+
+ int iC = int( angCurFloor );
+ if ( iC < nbAngles )
+ angle += ( angCur - angCurFloor ) * theAngles[ iC ];
+
+ iP = int( angPrevCeil );
+ while ( iC-- > iP )
+ angle += theAngles[ iC ];
+ }
+ res.push_back(angle);
+ angPrev = angCur;
+ }
+ Angles.swap( res );
+ }
+}
+
+//=======================================================================
+//function : linearScaleVariation
+//purpose : spread values over nbSteps
+//=======================================================================
+
+void SMESH_MeshEditor::linearScaleVariation(const int theNbSteps,
+ std::list<double>& theScales)
+{
+ int nbScales = theScales.size();
+ std::vector<double> myScales;
+ myScales.reserve( theNbSteps );
+ std::list<double>::const_iterator scale = theScales.begin();
+ double prevScale = 1.0;
+ for ( int iSc = 1; scale != theScales.end(); ++scale, ++iSc )
+ {
+ int iStep = int( iSc / double( nbScales ) * theNbSteps + 0.5 );
+ int stDelta = Max( 1, iStep - myScales.size());
+ double scDelta = ( *scale - prevScale ) / stDelta;
+ for ( int iStep = 0; iStep < stDelta; ++iStep )
+ {
+ myScales.push_back( prevScale + scDelta );
+ prevScale = myScales.back();
+ }
+ prevScale = *scale;
+ }
+ theScales.assign( myScales.begin(), myScales.end() );
+}
+
+//================================================================================
+/*!
+ * \brief Move or copy theElements applying theTrsf to their nodes
+ * \param theElems - elements to transform, if theElems is empty then apply to all mesh nodes
+ * \param theTrsf - transformation to apply
+ * \param theCopy - if true, create translated copies of theElems
+ * \param theMakeGroups - if true and theCopy, create translated groups
+ * \param theTargetMesh - mesh to copy translated elements into
+ * \return SMESH_MeshEditor::PGroupIDs - list of ids of created groups
+ */
+//================================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
+ const gp_Trsf& theTrsf,
+ const bool theCopy,
+ const bool theMakeGroups,
+ SMESH_Mesh* theTargetMesh)
+{
+ ClearLastCreated();
+ myLastCreatedElems.reserve( theElems.size() );
+
+ bool needReverse = false;
+ string groupPostfix;
+ switch ( theTrsf.Form() ) {
+ case gp_PntMirror:
+ needReverse = true;
+ groupPostfix = "mirrored";
+ break;
+ case gp_Ax1Mirror:
+ groupPostfix = "mirrored";
+ break;
+ case gp_Ax2Mirror:
+ needReverse = true;
+ groupPostfix = "mirrored";
+ break;
+ case gp_Rotation:
+ groupPostfix = "rotated";
+ break;
+ case gp_Translation:
+ groupPostfix = "translated";
+ break;
+ case gp_Scale:
+ groupPostfix = "scaled";
+ break;
+ case gp_CompoundTrsf: // different scale by axis
+ groupPostfix = "scaled";
+ break;
+ default:
+ needReverse = false;
+ groupPostfix = "transformed";
+ }
+
+ SMESHDS_Mesh* aTgtMesh = theTargetMesh ? theTargetMesh->GetMeshDS() : 0;
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ SMESH_MeshEditor targetMeshEditor( theTargetMesh );
+ SMESH_MeshEditor* editor = theTargetMesh ? & targetMeshEditor : theCopy ? this : 0;
+ SMESH_MeshEditor::ElemFeatures elemType;
+
+ // map old node to new one
+ TNodeNodeMap nodeMap;
+
+ // elements sharing moved nodes; those of them which have all
+ // nodes mirrored but are not in theElems are to be reversed
+ TIDSortedElemSet inverseElemSet;
+
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+
+ // issue 021015: EDF 1578 SMESH: Free nodes are removed when translating a mesh
+ TIDSortedElemSet orphanNode;
+
+ if ( theElems.empty() ) // transform the whole mesh
+ {
+ // add all elements
+ SMDS_ElemIteratorPtr eIt = aMesh->elementsIterator();
+ while ( eIt->more() ) theElems.insert( eIt->next() );
+ // add orphan nodes
+ SMDS_NodeIteratorPtr nIt = aMesh->nodesIterator();
+ while ( nIt->more() )
+ {
+ const SMDS_MeshNode* node = nIt->next();
+ if ( node->NbInverseElements() == 0)
+ orphanNode.insert( node );
+ }
+ }
+
+ // loop on elements to transform nodes : first orphan nodes then elems
+ TIDSortedElemSet::iterator itElem;
+ TIDSortedElemSet *elements[] = { &orphanNode, &theElems };
+ for (int i=0; i<2; i++)
+ for ( itElem = elements[i]->begin(); itElem != elements[i]->end(); itElem++ )
+ {
+ const SMDS_MeshElement* elem = *itElem;
+ if ( !elem )
+ continue;
+
+ // loop on elem nodes
+ double coord[3];
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() )
+ {
+ const SMDS_MeshNode* node = cast2Node( itN->next() );
+ // check if a node has been already transformed
+ pair<TNodeNodeMap::iterator,bool> n2n_isnew =
+ nodeMap.insert( make_pair ( node, node ));
+ if ( !n2n_isnew.second )
+ continue;
+
+ node->GetXYZ( coord );
+ theTrsf.Transforms( coord[0], coord[1], coord[2] );
+ if ( theTargetMesh ) {
+ const SMDS_MeshNode * newNode = aTgtMesh->AddNode( coord[0], coord[1], coord[2] );
+ n2n_isnew.first->second = newNode;
+ myLastCreatedNodes.push_back(newNode);
+ srcNodes.push_back( node );
+ }
+ else if ( theCopy ) {
+ const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+ n2n_isnew.first->second = newNode;
+ myLastCreatedNodes.push_back(newNode);
+ srcNodes.push_back( node );
+ }
+ else {
+ aMesh->MoveNode( node, coord[0], coord[1], coord[2] );
+ // node position on shape becomes invalid
+ const_cast< SMDS_MeshNode* > ( node )->SetPosition
+ ( SMDS_SpacePosition::originSpacePosition() );
+ }
+
+ // keep inverse elements
+ if ( !theCopy && !theTargetMesh && needReverse ) {
+ SMDS_ElemIteratorPtr invElemIt = node->GetInverseElementIterator();
+ while ( invElemIt->more() ) {
+ const SMDS_MeshElement* iel = invElemIt->next();
+ inverseElemSet.insert( iel );
+ }
+ }
+ }
+ } // loop on elems in { &orphanNode, &theElems };
+
+ // either create new elements or reverse mirrored ones
+ if ( !theCopy && !needReverse && !theTargetMesh )
+ return PGroupIDs();
+
+ theElems.insert( inverseElemSet.begin(),inverseElemSet.end() );
+
+ // Replicate or reverse elements
+
+ std::vector<int> iForw;
+ vector<const SMDS_MeshNode*> nodes;
+ for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
+ {
+ const SMDS_MeshElement* elem = *itElem;
+ if ( !elem ) continue;
+
+ SMDSAbs_GeometryType geomType = elem->GetGeomType();
+ size_t nbNodes = elem->NbNodes();
+ if ( geomType == SMDSGeom_NONE ) continue; // node
+
+ nodes.resize( nbNodes );
+
+ if ( geomType == SMDSGeom_POLYHEDRA ) // ------------------ polyhedral volume
+ {
+ const SMDS_MeshVolume* aPolyedre = SMDS_Mesh::DownCast< SMDS_MeshVolume >( elem );
+ if ( !aPolyedre )
+ continue;
+ nodes.clear();
+ 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
+ nodes.push_back((*nodeMapIt).second);
+ }
+ if ( needReverse && allTransformed )
+ std::reverse( nodes.end() - nbFaceNodes, nodes.end() );
+ }
+ if ( !allTransformed )
+ continue; // not all nodes transformed
+ }
+ else // ----------------------- the rest element types
+ {
+ while ( iForw.size() < nbNodes ) iForw.push_back( iForw.size() );
+ const vector<int>& iRev = SMDS_MeshCell::reverseSmdsOrder( elem->GetEntityType(), nbNodes );
+ const vector<int>& i = needReverse ? iRev : iForw;
+
+ // find transformed nodes
+ size_t 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 ( editor ) {
+ // copy in this or a new mesh
+ if ( editor->AddElement( nodes, elemType.Init( elem, /*basicOnly=*/false )))
+ srcElems.push_back( elem );
+ }
+ else {
+ // reverse element as it was reversed by transformation
+ if ( nbNodes > 2 )
+ aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes );
+ }
+
+ } // loop on elements
+
+ if ( editor && editor != this )
+ myLastCreatedElems.swap( editor->myLastCreatedElems );
+
+ PGroupIDs newGroupIDs;
+
+ if ( ( theMakeGroups && theCopy ) ||
+ ( theMakeGroups && theTargetMesh ) )
+ newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh, false );
+
+ return newGroupIDs;
+}
+
+//================================================================================
+/*!
+ * \brief Make an offset mesh from a source 2D mesh
+ * \param [in] theElements - source faces
+ * \param [in] theValue - offset value
+ * \param [out] theTgtMesh - a mesh to add offset elements to
+ * \param [in] theMakeGroups - to generate groups
+ * \return PGroupIDs - IDs of created groups. NULL means failure
+ */
+//================================================================================
+
+SMESH_MeshEditor::PGroupIDs SMESH_MeshEditor::Offset( TIDSortedElemSet & theElements,
+ const double theValue,
+ SMESH_Mesh* theTgtMesh,
+ const bool theMakeGroups,
+ const bool theCopyElements,
+ const bool theFixSelfIntersection)
+{
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+ SMESHDS_Mesh* tgtMeshDS = theTgtMesh->GetMeshDS();
+ SMESH_MeshEditor tgtEditor( theTgtMesh );
+
+ SMDS_ElemIteratorPtr eIt;
+ if ( theElements.empty() ) eIt = meshDS->elementsIterator( SMDSAbs_Face );
+ else eIt = SMESHUtils::elemSetIterator( theElements );
+
+ SMESH_MeshAlgos::TElemIntPairVec new2OldFaces;
+ SMESH_MeshAlgos::TNodeIntPairVec new2OldNodes;
+ std::unique_ptr< SMDS_Mesh > offsetMesh
+ ( SMESH_MeshAlgos::MakeOffset( eIt, *meshDS, theValue,
+ theFixSelfIntersection,
+ new2OldFaces, new2OldNodes ));
+ if ( offsetMesh->NbElements() == 0 )
+ return PGroupIDs(); // MakeOffset() failed
+
+
+ if ( theTgtMesh == myMesh && !theCopyElements )
+ {
+ // clear the source elements
+ if ( theElements.empty() ) eIt = meshDS->elementsIterator( SMDSAbs_Face );
+ else eIt = SMESHUtils::elemSetIterator( theElements );
+ while ( eIt->more() )
+ meshDS->RemoveFreeElement( eIt->next(), 0 );
+ }
+
+ // offsetMesh->Modified();
+ // offsetMesh->CompactMesh(); // make IDs start from 1
+
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+ srcElems.reserve( new2OldFaces.size() );
+ srcNodes.reserve( new2OldNodes.size() );
+
+ ClearLastCreated();
+ myLastCreatedElems.reserve( new2OldFaces.size() );
+ myLastCreatedNodes.reserve( new2OldNodes.size() );
+
+ // copy offsetMesh to theTgtMesh
+
+ int idShift = meshDS->MaxNodeID();
+ for ( size_t i = 0; i < new2OldNodes.size(); ++i )
+ if ( const SMDS_MeshNode* n = new2OldNodes[ i ].first )
+ {
+#ifndef _DEBUG_
+ if ( n->NbInverseElements() > 0 )
+#endif
+ {
+ const SMDS_MeshNode* n2 =
+ tgtMeshDS->AddNodeWithID( n->X(), n->Y(), n->Z(), idShift + n->GetID() );
+ myLastCreatedNodes.push_back( n2 );
+ srcNodes.push_back( meshDS->FindNode( new2OldNodes[ i ].second ));
+ }
+ }
+
+ ElemFeatures elemType;
+ for ( size_t i = 0; i < new2OldFaces.size(); ++i )
+ if ( const SMDS_MeshElement* f = new2OldFaces[ i ].first )
+ {
+ elemType.Init( f );
+ elemType.myNodes.clear();
+ for ( SMDS_NodeIteratorPtr nIt = f->nodeIterator(); nIt->more(); )
+ {
+ const SMDS_MeshNode* n2 = nIt->next();
+ elemType.myNodes.push_back( tgtMeshDS->FindNode( idShift + n2->GetID() ));
+ }
+ tgtEditor.AddElement( elemType.myNodes, elemType );
+ srcElems.push_back( meshDS->FindElement( new2OldFaces[ i ].second ));
+ }
+
+ myLastCreatedElems.swap( tgtEditor.myLastCreatedElems );
+
+ PGroupIDs newGroupIDs;
+ if ( theMakeGroups )
+ newGroupIDs = generateGroups( srcNodes, srcElems, "offset", theTgtMesh, false );
+ else
+ newGroupIDs.reset( new std::list< int > );
+
+ return newGroupIDs;
+}
+
+//=======================================================================
+/*!
+ * \brief Create groups of elements made during transformation
+ * \param nodeGens - nodes making corresponding myLastCreatedNodes
+ * \param elemGens - elements making corresponding myLastCreatedElems
+ * \param postfix - to push_back to names of new groups
+ * \param targetMesh - mesh to create groups in
+ * \param topPresent - is there are "top" elements that are created by sweeping
+ */
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
+ const SMESH_SequenceOfElemPtr& elemGens,
+ const std::string& postfix,
+ SMESH_Mesh* targetMesh,
+ const bool topPresent)
+{
+ PGroupIDs newGroupIDs( new list<int> );
+ SMESH_Mesh* mesh = targetMesh ? targetMesh : GetMesh();
+
+ // Sort existing groups by types and collect their names
+
+ // containers to store an old group and generated new ones;
+ // 1st new group is for result elems of different type than a source one;
+ // 2nd new group is for same type result elems ("top" group at extrusion)
+ 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;
+
+ SMESH_Mesh::GroupIteratorPtr groupIt = GetMesh()->GetGroups();
+ 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() );
+ groupDS->SetStoreName( group->GetName() );
+ 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[ type ].push_back( make_tuple( groupDS, newGroup, newTopGroup ));
+ orderedOldNewGroups.push_back( & groupsByType[ type ].back() );
+ }
+
+ // Loop on nodes and elements to add them in new groups
+
+ vector< const SMDS_MeshElement* > resultElems;
+ for ( int isNodes = 0; isNodes < 2; ++isNodes )
+ {
+ const SMESH_SequenceOfElemPtr& gens = isNodes ? nodeGens : elemGens;
+ const SMESH_SequenceOfElemPtr& elems = isNodes ? myLastCreatedNodes : myLastCreatedElems;
+ if ( gens.size() != elems.size() )
+ throw SALOME_Exception("SMESH_MeshEditor::generateGroups(): invalid args");
+
+ // loop on created elements
+ for (size_t iElem = 0; iElem < elems.size(); ++iElem )
+ {
+ const SMDS_MeshElement* sourceElem = gens[ iElem ];
+ if ( !sourceElem ) {
+ MESSAGE("generateGroups(): NULL source element");
+ continue;
+ }
+ list< TOldNewGroup > & groupsOldNew = groupsByType[ sourceElem->GetType() ];
+ if ( groupsOldNew.empty() ) { // no groups of this type at all
+ while ( iElem+1 < gens.size() && gens[ iElem+1 ] == sourceElem )
+ ++iElem; // skip all elements made by sourceElem
+ continue;
+ }
+ // collect all elements made by the iElem-th sourceElem
+ resultElems.clear();
+ if ( const SMDS_MeshElement* resElem = elems[ iElem ])
+ if ( resElem != sourceElem )
+ resultElems.push_back( resElem );
+ while ( iElem+1 < gens.size() && gens[ iElem+1 ] == sourceElem )
+ if ( const SMDS_MeshElement* resElem = elems[ ++iElem ])
+ if ( resElem != sourceElem )
+ resultElems.push_back( resElem );
+
+ const SMDS_MeshElement* topElem = 0;
+ if ( isNodes ) // there must be a top element
+ {
+ topElem = resultElems.back();
+ resultElems.pop_back();
+ }
+ else
+ {
+ vector< const SMDS_MeshElement* >::reverse_iterator resElemIt = resultElems.rbegin();
+ for ( ; resElemIt != resultElems.rend() ; ++resElemIt )
+ if ( (*resElemIt)->GetType() == sourceElem->GetType() )
+ {
+ topElem = *resElemIt;
+ *resElemIt = 0; // erase *resElemIt
+ break;
+ }
+ }
+ // add resultElems to groups originted from ones the sourceElem belongs to
+ list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end();
+ for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew )
+ {
+ SMESHDS_GroupBase* oldGroup = gOldNew->get<0>();
+ if ( oldGroup->Contains( sourceElem )) // sourceElem is in oldGroup
+ {
+ // fill in a new group
+ SMDS_MeshGroup & newGroup = gOldNew->get<1>()->SMDSGroup();
+ vector< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
+ for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt )
+ if ( *resElemIt )
+ newGroup.Add( *resElemIt );
+
+ // fill a "top" group
+ if ( topElem )
+ {
+ SMDS_MeshGroup & newTopGroup = gOldNew->get<2>()->SMDSGroup();
+ newTopGroup.Add( topElem );
+ }
+ }
+ }
+ } // 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>() };
+ 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 = ( topPresent &&
+ newGroupDS->GetType() == oldGroupDS->GetType() &&
+ is2nd );
+
+ string name = oldGroupDS->GetStoreName();
+ { // remove trailing whitespaces (issue 22599)
+ size_t size = name.size();
+ while ( size > 1 && isspace( name[ size-1 ]))
+ --size;
+ if ( size != name.size() )
+ {
+ name.resize( size );
+ oldGroupDS->SetStoreName( name.c_str() );
+ }
+ }
+ 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;
+}
+
+//================================================================================
+/*!
+ * * \brief Return list of group of nodes close to each other within theTolerance
+ * * Search among theNodes or in the whole mesh if theNodes is empty using
+ * * an Octree algorithm
+ * \param [in,out] theNodes - the nodes to treat
+ * \param [in] theTolerance - the tolerance
+ * \param [out] theGroupsOfNodes - the result groups of coincident nodes
+ * \param [in] theSeparateCornersAndMedium - if \c true, in quadratic mesh puts
+ * corner and medium nodes in separate groups
+ */
+//================================================================================
+
+void SMESH_MeshEditor::FindCoincidentNodes (TIDSortedNodeSet & theNodes,
+ const double theTolerance,
+ TListOfListOfNodes & theGroupsOfNodes,
+ bool theSeparateCornersAndMedium)
+{
+ ClearLastCreated();
+
+ if ( myMesh->NbEdges ( ORDER_QUADRATIC ) +
+ myMesh->NbFaces ( ORDER_QUADRATIC ) +
+ myMesh->NbVolumes( ORDER_QUADRATIC ) == 0 )
+ theSeparateCornersAndMedium = false;
+
+ TIDSortedNodeSet& corners = theNodes;
+ TIDSortedNodeSet medium;
+
+ if ( theNodes.empty() ) // get all nodes in the mesh
+ {
+ TIDSortedNodeSet* nodes[2] = { &corners, &medium };
+ SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator();
+ if ( theSeparateCornersAndMedium )
+ while ( nIt->more() )
+ {
+ const SMDS_MeshNode* n = nIt->next();
+ TIDSortedNodeSet* & nodeSet = nodes[ SMESH_MesherHelper::IsMedium( n )];
+ nodeSet->insert( nodeSet->end(), n );
+ }
+ else
+ while ( nIt->more() )
+ theNodes.insert( theNodes.end(), nIt->next() );
+ }
+ else if ( theSeparateCornersAndMedium ) // separate corners from medium nodes
+ {
+ TIDSortedNodeSet::iterator nIt = corners.begin();
+ while ( nIt != corners.end() )
+ if ( SMESH_MesherHelper::IsMedium( *nIt ))
+ {
+ medium.insert( medium.end(), *nIt );
+ corners.erase( nIt++ );
+ }
+ else
+ {
+ ++nIt;
+ }
+ }
+
+ if ( !corners.empty() )
+ SMESH_OctreeNode::FindCoincidentNodes ( corners, &theGroupsOfNodes, theTolerance );
+ if ( !medium.empty() )
+ SMESH_OctreeNode::FindCoincidentNodes ( medium, &theGroupsOfNodes, theTolerance );
+}
+
+//=======================================================================
+//function : SimplifyFace
+//purpose : split a chain of nodes into several closed chains
+//=======================================================================
+
+int SMESH_MeshEditor::SimplifyFace (const vector<const SMDS_MeshNode *>& faceNodes,
+ vector<const SMDS_MeshNode *>& poly_nodes,
+ vector<int>& quantities) const
+{
+ int nbNodes = faceNodes.size();
+ while ( faceNodes[ 0 ] == faceNodes[ nbNodes-1 ] && nbNodes > 2 )
+ --nbNodes;
+ if ( nbNodes < 3 )
+ return 0;
+ size_t prevNbQuant = quantities.size();
+
+ vector< const SMDS_MeshNode* > simpleNodes; simpleNodes.reserve( nbNodes );
+ map< const SMDS_MeshNode*, int > nodeIndices; // indices within simpleNodes
+ map< const SMDS_MeshNode*, int >::iterator nInd;
+
+ nodeIndices.insert( make_pair( faceNodes[0], 0 ));
+ simpleNodes.push_back( faceNodes[0] );
+ for ( int iCur = 1; iCur < nbNodes; iCur++ )
+ {
+ if ( faceNodes[ iCur ] != simpleNodes.back() )
+ {
+ int index = simpleNodes.size();
+ nInd = nodeIndices.insert( make_pair( faceNodes[ iCur ], index )).first;
+ int prevIndex = nInd->second;
+ if ( prevIndex < index )
+ {
+ // a sub-loop found
+ int loopLen = index - prevIndex;
+ if ( loopLen > 2 )
+ {
+ // store the sub-loop
+ quantities.push_back( loopLen );
+ for ( int i = prevIndex; i < index; i++ )
+ poly_nodes.push_back( simpleNodes[ i ]);
+ }
+ simpleNodes.resize( prevIndex+1 );
+ }
+ else
+ {
+ simpleNodes.push_back( faceNodes[ iCur ]);
+ }
+ }
+ }
+
+ if ( simpleNodes.size() > 2 )
+ {
+ quantities.push_back( simpleNodes.size() );
+ poly_nodes.insert ( poly_nodes.end(), simpleNodes.begin(), simpleNodes.end() );
+ }
+
+ return quantities.size() - prevNbQuant;
+}
+
+//=======================================================================
+//function : MergeNodes
+//purpose : In each group, the cdr of nodes are substituted by the first one
+// in all elements.
+//=======================================================================
+
+void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes,
+ const bool theAvoidMakingHoles)
+{
+ ClearLastCreated();
+
+ SMESHDS_Mesh* mesh = GetMeshDS();
+
+ TNodeNodeMap nodeNodeMap; // node to replace - new node
+ set<const SMDS_MeshElement*> elems; // all elements with changed nodes
+ list< int > rmElemIds, rmNodeIds;
+ vector< ElemFeatures > newElemDefs;
+
+ // Fill nodeNodeMap and elems
+
+ TListOfListOfNodes::iterator grIt = theGroupsOfNodes.begin();
+ for ( ; grIt != theGroupsOfNodes.end(); grIt++ )
+ {
+ list<const SMDS_MeshNode*>& nodes = *grIt;
+ list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
+ const SMDS_MeshNode* nToKeep = *nIt;
+ for ( ++nIt; nIt != nodes.end(); nIt++ )
+ {
+ const SMDS_MeshNode* nToRemove = *nIt;
+ nodeNodeMap.insert( make_pair( nToRemove, nToKeep ));
+ SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
+ while ( invElemIt->more() ) {
+ const SMDS_MeshElement* elem = invElemIt->next();
+ elems.insert(elem);
+ }
+ }
+ }
+
+ // Apply recursive replacements (BUG 0020185)
+ TNodeNodeMap::iterator nnIt = nodeNodeMap.begin();
+ for ( ; nnIt != nodeNodeMap.end(); ++nnIt )
+ {
+ const SMDS_MeshNode* nToKeep = nnIt->second;
+ TNodeNodeMap::iterator nnIt_i = nodeNodeMap.find( nToKeep );
+ while ( nnIt_i != nodeNodeMap.end() && nnIt_i->second != nnIt->second )
+ {
+ nToKeep = nnIt_i->second;
+ nnIt->second = nToKeep;
+ nnIt_i = nodeNodeMap.find( nToKeep );
+ }
+ }
+
+ if ( theAvoidMakingHoles )
+ {
+ // find elements whose topology changes
+
+ vector<const SMDS_MeshElement*> pbElems;
+ set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
+ for ( ; eIt != elems.end(); ++eIt )
+ {
+ const SMDS_MeshElement* elem = *eIt;
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() )
+ {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( itN->next() );
+ TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
+ if ( nnIt != nodeNodeMap.end() && elem->GetNodeIndex( nnIt->second ) >= 0 )
+ {
+ // several nodes of elem stick
+ pbElems.push_back( elem );
+ break;
+ }
+ }
+ }
+ // exclude from merge nodes causing spoiling element
+ for ( size_t iLoop = 0; iLoop < pbElems.size(); ++iLoop ) // avoid infinite cycle
+ {
+ bool nodesExcluded = false;
+ for ( size_t i = 0; i < pbElems.size(); ++i )
+ {
+ size_t prevNbMergeNodes = nodeNodeMap.size();
+ if ( !applyMerge( pbElems[i], newElemDefs, nodeNodeMap, /*noHoles=*/true ) &&
+ prevNbMergeNodes < nodeNodeMap.size() )
+ nodesExcluded = true;
+ }
+ if ( !nodesExcluded )
+ break;
+ }
+ }
+
+ for ( nnIt = nodeNodeMap.begin(); nnIt != nodeNodeMap.end(); ++nnIt )
+ {
+ const SMDS_MeshNode* nToRemove = nnIt->first;
+ const SMDS_MeshNode* nToKeep = nnIt->second;
+ if ( nToRemove != nToKeep )
+ {
+ rmNodeIds.push_back( nToRemove->GetID() );
+ AddToSameGroups( nToKeep, nToRemove, mesh );
+ // set _alwaysComputed to a sub-mesh of VERTEX to enable further mesh computing
+ // w/o creating node in place of merged ones.
+ SMDS_PositionPtr pos = nToRemove->GetPosition();
+ if ( pos && pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
+ if ( SMESH_subMesh* sm = myMesh->GetSubMeshContaining( nToRemove->getshapeId() ))
+ sm->SetIsAlwaysComputed( true );
+ }
+ }
+
+ // Change element nodes or remove an element
+
+ set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
+ for ( ; eIt != elems.end(); eIt++ )
+ {
+ const SMDS_MeshElement* elem = *eIt;
+ SMESHDS_SubMesh* sm = mesh->MeshElements( elem->getshapeId() );
+
+ bool keepElem = applyMerge( elem, newElemDefs, nodeNodeMap, /*noHoles=*/false );
+ if ( !keepElem )
+ rmElemIds.push_back( elem->GetID() );
+
+ for ( size_t i = 0; i < newElemDefs.size(); ++i )
+ {
+ if ( i > 0 || !mesh->ChangeElementNodes( elem,
+ & newElemDefs[i].myNodes[0],
+ newElemDefs[i].myNodes.size() ))
+ {
+ if ( i == 0 )
+ {
+ newElemDefs[i].SetID( elem->GetID() );
+ mesh->RemoveFreeElement(elem, sm, /*fromGroups=*/false);
+ if ( !keepElem ) rmElemIds.pop_back();
+ }
+ else
+ {
+ newElemDefs[i].SetID( -1 );
+ }
+ SMDS_MeshElement* newElem = this->AddElement( newElemDefs[i].myNodes, newElemDefs[i] );
+ if ( sm && newElem )
+ sm->AddElement( newElem );
+ if ( elem != newElem )
+ ReplaceElemInGroups( elem, newElem, mesh );
+ }
+ }
+ }
+
+ // Remove bad elements, then equal nodes (order important)
+ Remove( rmElemIds, /*isNodes=*/false );
+ Remove( rmNodeIds, /*isNodes=*/true );
+
+ return;
+}
+
+//=======================================================================
+//function : applyMerge
+//purpose : Compute new connectivity of an element after merging nodes
+// \param [in] elems - the element
+// \param [out] newElemDefs - definition(s) of result element(s)
+// \param [inout] nodeNodeMap - nodes to merge
+// \param [in] avoidMakingHoles - if true and and the element becomes invalid
+// after merging (but not degenerated), removes nodes causing
+// the invalidity from \a nodeNodeMap.
+// \return bool - true if the element should be removed
+//=======================================================================
+
+bool SMESH_MeshEditor::applyMerge( const SMDS_MeshElement* elem,
+ vector< ElemFeatures >& newElemDefs,
+ TNodeNodeMap& nodeNodeMap,
+ const bool avoidMakingHoles )
+{
+ bool toRemove = false; // to remove elem
+ int nbResElems = 1; // nb new elements
+
+ newElemDefs.resize(nbResElems);
+ newElemDefs[0].Init( elem );
+ newElemDefs[0].myNodes.clear();
+
+ set<const SMDS_MeshNode*> nodeSet;
+ vector< const SMDS_MeshNode*> curNodes;
+ vector< const SMDS_MeshNode*> & uniqueNodes = newElemDefs[0].myNodes;
+ vector<int> iRepl;
+
+ const int nbNodes = elem->NbNodes();
+ SMDSAbs_EntityType entity = elem->GetEntityType();
+
+ curNodes.resize( nbNodes );
+ uniqueNodes.resize( nbNodes );
+ iRepl.resize( nbNodes );
+ int iUnique = 0, iCur = 0, nbRepl = 0;
+
+ // Get new seq of nodes
+
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() )
+ {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( itN->next() );
+
+ TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
+ if ( nnIt != nodeNodeMap.end() ) {
+ n = (*nnIt).second;
+ }
+ curNodes[ iCur ] = n;
+ bool isUnique = nodeSet.insert( n ).second;
+ if ( isUnique )
+ uniqueNodes[ iUnique++ ] = n;
+ else
+ iRepl[ nbRepl++ ] = iCur;
+ iCur++;
+ }
+
+ // Analyse element topology after replacement
+
+ int nbUniqueNodes = nodeSet.size();
+ if ( nbNodes != nbUniqueNodes ) // some nodes stick
+ {
+ toRemove = true;
+ nbResElems = 0;
+
+ if ( newElemDefs[0].myIsQuad && newElemDefs[0].myType == SMDSAbs_Face && nbNodes > 6 )
+ {
+ // if corner nodes stick, remove medium nodes between them from uniqueNodes
+ int nbCorners = nbNodes / 2;
+ for ( int iCur = 0; iCur < nbCorners; ++iCur )
+ {
+ int iNext = ( iCur + 1 ) % nbCorners;
+ if ( curNodes[ iCur ] == curNodes[ iNext ] ) // corners stick
+ {
+ int iMedium = iCur + nbCorners;
+ vector< const SMDS_MeshNode* >::iterator i =
+ std::find( uniqueNodes.begin() + nbCorners - nbRepl,
+ uniqueNodes.end(),
+ curNodes[ iMedium ]);
+ if ( i != uniqueNodes.end() )
+ {
+ --nbUniqueNodes;
+ for ( ; i+1 != uniqueNodes.end(); ++i )
+ *i = *(i+1);
+ }
+ }
+ }
+ }
+
+ switch ( entity )
+ {
+ case SMDSEntity_Polygon:
+ case SMDSEntity_Quad_Polygon: // Polygon
+ {
+ ElemFeatures* elemType = & newElemDefs[0];
+ const bool isQuad = elemType->myIsQuad;
+ if ( isQuad )
+ SMDS_MeshCell::applyInterlace // interlace medium and corner nodes
+ ( SMDS_MeshCell::interlacedSmdsOrder( SMDSEntity_Quad_Polygon, nbNodes ), curNodes );
+
+ // a polygon can divide into several elements
+ vector<const SMDS_MeshNode *> polygons_nodes;
+ vector<int> quantities;
+ nbResElems = SimplifyFace( curNodes, polygons_nodes, quantities );
+ newElemDefs.resize( nbResElems );
+ for ( int inode = 0, iface = 0; iface < nbResElems; iface++ )
+ {
+ ElemFeatures* elemType = & newElemDefs[iface];
+ if ( iface ) elemType->Init( elem );
+
+ vector<const SMDS_MeshNode *>& face_nodes = elemType->myNodes;
+ int nbNewNodes = quantities[iface];
+ face_nodes.assign( polygons_nodes.begin() + inode,
+ polygons_nodes.begin() + inode + nbNewNodes );
+ inode += nbNewNodes;
+ if ( isQuad ) // check if a result elem is a valid quadratic polygon
+ {
+ bool isValid = ( nbNewNodes % 2 == 0 );
+ for ( int i = 0; i < nbNewNodes && isValid; ++i )
+ isValid = ( elem->IsMediumNode( face_nodes[i]) == bool( i % 2 ));
+ elemType->SetQuad( isValid );
+ if ( isValid ) // put medium nodes after corners
+ SMDS_MeshCell::applyInterlaceRev
+ ( SMDS_MeshCell::interlacedSmdsOrder( SMDSEntity_Quad_Polygon,
+ nbNewNodes ), face_nodes );
+ }
+ elemType->SetPoly(( nbNewNodes / ( elemType->myIsQuad + 1 ) > 4 ));
+ }
+ nbUniqueNodes = newElemDefs[0].myNodes.size();
+ break;
+ } // Polygon
+
+ case SMDSEntity_Polyhedra: // Polyhedral volume
+ {
+ if ( nbUniqueNodes >= 4 )
+ {
+ // each face has to be analyzed in order to check volume validity
+ if ( const SMDS_MeshVolume* aPolyedre = SMDS_Mesh::DownCast< SMDS_MeshVolume >( elem ))
+ {
+ int nbFaces = aPolyedre->NbFaces();
+
+ vector<const SMDS_MeshNode *>& poly_nodes = newElemDefs[0].myNodes;
+ vector<int> & quantities = newElemDefs[0].myPolyhedQuantities;
+ vector<const SMDS_MeshNode *> faceNodes;
+ poly_nodes.clear();
+ quantities.clear();
+
+ for (int iface = 1; iface <= nbFaces; iface++)
+ {
+ int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
+ faceNodes.resize( nbFaceNodes );
+ for (int inode = 1; inode <= nbFaceNodes; inode++)
+ {
+ const SMDS_MeshNode * faceNode = aPolyedre->GetFaceNode(iface, inode);
+ TNodeNodeMap::iterator nnIt = nodeNodeMap.find(faceNode);
+ if ( nnIt != nodeNodeMap.end() ) // faceNode sticks
+ faceNode = (*nnIt).second;
+ faceNodes[inode - 1] = faceNode;
+ }
+ SimplifyFace(faceNodes, poly_nodes, quantities);
+ }
+
+ if ( quantities.size() > 3 )
+ {
+ // TODO: remove coincident faces
+ nbResElems = 1;
+ nbUniqueNodes = newElemDefs[0].myNodes.size();
+ }
+ }
+ }
+ }
+ break;
+
+ // Regular elements
+ // TODO not all the possible cases are solved. Find something more generic?
+ case SMDSEntity_Edge: //////// EDGE
+ case SMDSEntity_Triangle: //// TRIANGLE
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Tetra:
+ case SMDSEntity_Quad_Tetra: // TETRAHEDRON
+ {
+ break;
+ }
+ case SMDSEntity_Quad_Edge:
+ {
+ break;
+ }
+ case SMDSEntity_Quadrangle: //////////////////////////////////// QUADRANGLE
+ {
+ if ( nbUniqueNodes < 3 )
+ toRemove = true;
+ else if ( nbRepl == 1 && curNodes[ iRepl[0]] == curNodes[( iRepl[0]+2 )%4 ])
+ toRemove = true; // opposite nodes stick
+ else
+ toRemove = false;
+ break;
+ }
+ case SMDSEntity_Quad_Quadrangle: // Quadratic QUADRANGLE
+ {
+ // 1 5 2
+ // +---+---+
+ // | |
+ // 4+ +6
+ // | |
+ // +---+---+
+ // 0 7 3
+ if ( nbUniqueNodes == 6 &&
+ iRepl[0] < 4 &&
+ ( nbRepl == 1 || iRepl[1] >= 4 ))
+ {
+ toRemove = false;
+ }
+ break;
+ }
+ case SMDSEntity_BiQuad_Quadrangle: // Bi-Quadratic QUADRANGLE
+ {
+ // 1 5 2
+ // +---+---+
+ // | |
+ // 4+ 8+ +6
+ // | |
+ // +---+---+
+ // 0 7 3
+ if ( nbUniqueNodes == 7 &&
+ iRepl[0] < 4 &&
+ ( nbRepl == 1 || iRepl[1] != 8 ))
+ {
+ toRemove = false;
+ }
+ break;
+ }
+ case SMDSEntity_Penta: ///////////////////////////////////// PENTAHEDRON
+ {
+ if ( nbUniqueNodes == 4 ) {
+ // ---------------------------------> tetrahedron
+ if ( curNodes[3] == curNodes[4] &&
+ curNodes[3] == curNodes[5] ) {
+ // top nodes stick
+ toRemove = false;
+ }
+ else if ( curNodes[0] == curNodes[1] &&
+ curNodes[0] == curNodes[2] ) {
+ // bottom nodes stick: set a top before
+ uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
+ uniqueNodes[ 0 ] = curNodes [ 5 ];
+ uniqueNodes[ 1 ] = curNodes [ 4 ];
+ uniqueNodes[ 2 ] = curNodes [ 3 ];
+ toRemove = false;
+ }
+ else if (( curNodes[0] == curNodes[3] ) +
+ ( curNodes[1] == curNodes[4] ) +
+ ( curNodes[2] == curNodes[5] ) == 2 ) {
+ // a lateral face turns into a line
+ toRemove = false;
+ }
+ }
+ else if ( nbUniqueNodes == 5 ) {
+ // PENTAHEDRON --------------------> pyramid
+ if ( curNodes[0] == curNodes[3] )
+ {
+ uniqueNodes[ 0 ] = curNodes[ 1 ];
+ uniqueNodes[ 1 ] = curNodes[ 4 ];
+ uniqueNodes[ 2 ] = curNodes[ 5 ];
+ uniqueNodes[ 3 ] = curNodes[ 2 ];
+ uniqueNodes[ 4 ] = curNodes[ 0 ];
+ toRemove = false;
+ }
+ if ( curNodes[1] == curNodes[4] )
+ {
+ uniqueNodes[ 0 ] = curNodes[ 0 ];
+ uniqueNodes[ 1 ] = curNodes[ 2 ];
+ uniqueNodes[ 2 ] = curNodes[ 5 ];
+ uniqueNodes[ 3 ] = curNodes[ 3 ];
+ uniqueNodes[ 4 ] = curNodes[ 1 ];
+ toRemove = false;
+ }
+ if ( curNodes[2] == curNodes[5] )
+ {
+ uniqueNodes[ 0 ] = curNodes[ 0 ];
+ uniqueNodes[ 1 ] = curNodes[ 3 ];
+ uniqueNodes[ 2 ] = curNodes[ 4 ];
+ uniqueNodes[ 3 ] = curNodes[ 1 ];
+ uniqueNodes[ 4 ] = curNodes[ 2 ];
+ toRemove = false;
+ }
+ }
+ break;
+ }
+ case SMDSEntity_Hexa:
+ {
+ //////////////////////////////////// HEXAHEDRON
+ SMDS_VolumeTool hexa (elem);
+ hexa.SetExternalNormal();
+ if ( nbUniqueNodes == 4 && nbRepl == 4 ) {
+ //////////////////////// HEX ---> tetrahedron
+ for ( int iFace = 0; iFace < 6; iFace++ ) {
+ const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
+ if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
+ curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
+ curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
+ // one face turns into a point ...
+ int pickInd = ind[ 0 ];
+ int iOppFace = hexa.GetOppFaceIndex( iFace );
+ ind = hexa.GetFaceNodesIndices( iOppFace );
+ int nbStick = 0;
+ uniqueNodes.clear();
+ for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
+ if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
+ nbStick++;
+ else
+ uniqueNodes.push_back( curNodes[ind[ iCur ]]);
+ }
+ if ( nbStick == 1 ) {
+ // ... and the opposite one - into a triangle.
+ // set a top node
+ uniqueNodes.push_back( curNodes[ pickInd ]);
+ toRemove = false;
+ }
+ break;
+ }
+ }
+ }
+ else if ( nbUniqueNodes == 6 && nbRepl == 2 ) {
+ //////////////////////// HEX ---> prism
+ int nbTria = 0, iTria[3];
+ const int *ind; // indices of face nodes
+ // look for triangular faces
+ for ( int iFace = 0; iFace < 6 && nbTria < 3; iFace++ ) {
+ ind = hexa.GetFaceNodesIndices( iFace );
+ TIDSortedNodeSet faceNodes;
+ for ( iCur = 0; iCur < 4; iCur++ )
+ faceNodes.insert( curNodes[ind[iCur]] );
+ if ( faceNodes.size() == 3 )
+ iTria[ nbTria++ ] = iFace;
+ }
+ // check if triangles are opposite
+ if ( nbTria == 2 && iTria[0] == hexa.GetOppFaceIndex( iTria[1] ))
+ {
+ // set nodes of the bottom triangle
+ ind = hexa.GetFaceNodesIndices( iTria[ 0 ]);
+ vector<int> indB;
+ for ( iCur = 0; iCur < 4; iCur++ )
+ if ( ind[iCur] != iRepl[0] && ind[iCur] != iRepl[1])
+ indB.push_back( ind[iCur] );
+ if ( !hexa.IsForward() )
+ std::swap( indB[0], indB[2] );
+ for ( iCur = 0; iCur < 3; iCur++ )
+ uniqueNodes[ iCur ] = curNodes[indB[iCur]];
+ // set nodes of the top triangle
+ const int *indT = hexa.GetFaceNodesIndices( iTria[ 1 ]);
+ for ( iCur = 0; iCur < 3; ++iCur )
+ for ( int j = 0; j < 4; ++j )
+ if ( hexa.IsLinked( indB[ iCur ], indT[ j ] ))
+ {
+ uniqueNodes[ iCur + 3 ] = curNodes[ indT[ j ]];
+ break;
+ }
+ toRemove = false;
+ break;
+ }
+ }
+ else if (nbUniqueNodes == 5 && nbRepl == 3 ) {
+ //////////////////// HEXAHEDRON ---> pyramid
+ for ( int iFace = 0; iFace < 6; iFace++ ) {
+ const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
+ if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
+ curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
+ curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
+ // one face turns into a point ...
+ int iOppFace = hexa.GetOppFaceIndex( iFace );
+ ind = hexa.GetFaceNodesIndices( iOppFace );
+ uniqueNodes.clear();
+ for ( iCur = 0; iCur < 4; iCur++ ) {
+ if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
+ break;
+ else
+ uniqueNodes.push_back( curNodes[ind[ iCur ]]);
+ }
+ if ( uniqueNodes.size() == 4 ) {
+ // ... and the opposite one is a quadrangle
+ // set a top node
+ const int* indTop = hexa.GetFaceNodesIndices( iFace );
+ uniqueNodes.push_back( curNodes[indTop[ 0 ]]);
+ toRemove = false;
+ }
+ break;
+ }
+ }
+ }
+
+ if ( toRemove && nbUniqueNodes > 4 ) {
+ ////////////////// HEXAHEDRON ---> polyhedron
+ hexa.SetExternalNormal();
+ vector<const SMDS_MeshNode *>& poly_nodes = newElemDefs[0].myNodes;
+ vector<int> & quantities = newElemDefs[0].myPolyhedQuantities;
+ poly_nodes.reserve( 6 * 4 ); poly_nodes.clear();
+ quantities.reserve( 6 ); quantities.clear();
+ for ( int iFace = 0; iFace < 6; iFace++ )
+ {
+ const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
+ if ( curNodes[ind[0]] == curNodes[ind[2]] ||
+ curNodes[ind[1]] == curNodes[ind[3]] )
+ {
+ quantities.clear();
+ break; // opposite nodes stick
+ }
+ nodeSet.clear();
+ for ( iCur = 0; iCur < 4; iCur++ )
+ {
+ if ( nodeSet.insert( curNodes[ind[ iCur ]] ).second )
+ poly_nodes.push_back( curNodes[ind[ iCur ]]);
+ }
+ if ( nodeSet.size() < 3 )
+ poly_nodes.resize( poly_nodes.size() - nodeSet.size() );
+ else
+ quantities.push_back( nodeSet.size() );
+ }
+ if ( quantities.size() >= 4 )
+ {
+ nbResElems = 1;
+ nbUniqueNodes = poly_nodes.size();
+ newElemDefs[0].SetPoly(true);
+ }
+ }
+ break;
+ } // case HEXAHEDRON
+
+ default:
+ toRemove = true;
+
+ } // switch ( entity )
+
+ if ( toRemove && nbResElems == 0 && avoidMakingHoles )
+ {
+ // erase from nodeNodeMap nodes whose merge spoils elem
+ vector< const SMDS_MeshNode* > noMergeNodes;
+ SMESH_MeshAlgos::DeMerge( elem, curNodes, noMergeNodes );
+ for ( size_t i = 0; i < noMergeNodes.size(); ++i )
+ nodeNodeMap.erase( noMergeNodes[i] );
+ }
+
+ } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
+
+ uniqueNodes.resize( nbUniqueNodes );
+
+ if ( !toRemove && nbResElems == 0 )
+ nbResElems = 1;
+
+ newElemDefs.resize( nbResElems );
+
+ return !toRemove;
+}
+
+
+// ========================================================
+// class : ComparableElement
+// purpose : allow comparing elements basing on their nodes
+// ========================================================
+
+class ComparableElement : public boost::container::flat_set< int >
+{
+ typedef boost::container::flat_set< int > int_set;
+
+ const SMDS_MeshElement* myElem;
+ int mySumID;
+ mutable int myGroupID;
+
+public:
+
+ ComparableElement( const SMDS_MeshElement* theElem ):
+ myElem ( theElem ), mySumID( 0 ), myGroupID( -1 )
+ {
+ this->reserve( theElem->NbNodes() );
+ for ( SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator(); nodeIt->more(); )
+ {
+ int id = nodeIt->next()->GetID();
+ mySumID += id;
+ this->insert( id );
+ }
+ }
+
+ const SMDS_MeshElement* GetElem() const { return myElem; }
+
+ int& GroupID() const { return myGroupID; }
+ //int& GroupID() const { return const_cast< int& >( myGroupID ); }
+
+ ComparableElement( const ComparableElement& theSource ) // move copy
+ {
+ ComparableElement& src = const_cast< ComparableElement& >( theSource );
+ (int_set&) (*this ) = boost::move( src );
+ myElem = src.myElem;
+ mySumID = src.mySumID;
+ myGroupID = src.myGroupID;
+ }
+
+ static int HashCode(const ComparableElement& se, int limit )
+ {
+ return ::HashCode( se.mySumID, limit );
+ }
+ static Standard_Boolean IsEqual(const ComparableElement& se1, const ComparableElement& se2 )
+ {
+ return ( se1 == se2 );
+ }
+
+};
+
+//=======================================================================
+//function : FindEqualElements
+//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( TIDSortedElemSet & theElements,
+ TListOfListOfElementsID & theGroupsOfElementsID )
+{
+ ClearLastCreated();
+
+ SMDS_ElemIteratorPtr elemIt;
+ if ( theElements.empty() ) elemIt = GetMeshDS()->elementsIterator();
+ else elemIt = SMESHUtils::elemSetIterator( theElements );
+
+ typedef NCollection_Map< ComparableElement, ComparableElement > TMapOfElements;
+ typedef std::list<int> TGroupOfElems;
+ TMapOfElements mapOfElements;
+ std::vector< TGroupOfElems > arrayOfGroups;
+ TGroupOfElems groupOfElems;
+
+ while ( elemIt->more() )
+ {
+ const SMDS_MeshElement* curElem = elemIt->next();
+ if ( curElem->IsNull() )
+ continue;
+ ComparableElement compElem = curElem;
+ // check uniqueness
+ const ComparableElement& elemInSet = mapOfElements.Added( compElem );
+ if ( elemInSet.GetElem() != curElem ) // coincident elem
+ {
+ int& iG = elemInSet.GroupID();
+ if ( iG < 0 )
+ {
+ iG = arrayOfGroups.size();
+ arrayOfGroups.push_back( groupOfElems );
+ arrayOfGroups[ iG ].push_back( elemInSet.GetElem()->GetID() );
+ }
+ arrayOfGroups[ iG ].push_back( curElem->GetID() );
+ }
+ }
+
+ groupOfElems.clear();
+ std::vector< TGroupOfElems >::iterator groupIt = arrayOfGroups.begin();
+ for ( ; groupIt != arrayOfGroups.end(); ++groupIt )
+ {
+ if ( groupIt->size() > 1 ) {
+ //groupOfElems.sort(); -- theElements are sorted already
+ theGroupsOfElementsID.emplace_back( *groupIt );
+ }
+ }
+}
+
+//=======================================================================
+//function : MergeElements
+//purpose : In each given group, substitute all elements by the first one.
+//=======================================================================
+
+void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElementsID)
+{
+ ClearLastCreated();
+
+ typedef list<int> TListOfIDs;
+ TListOfIDs rmElemIds; // IDs of elems to remove
+
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ TListOfListOfElementsID::iterator groupsIt = theGroupsOfElementsID.begin();
+ while ( groupsIt != theGroupsOfElementsID.end() ) {
+ TListOfIDs& aGroupOfElemID = *groupsIt;
+ aGroupOfElemID.sort();
+ int elemIDToKeep = aGroupOfElemID.front();
+ const SMDS_MeshElement* elemToKeep = aMesh->FindElement(elemIDToKeep);
+ aGroupOfElemID.pop_front();
+ TListOfIDs::iterator idIt = aGroupOfElemID.begin();
+ while ( idIt != aGroupOfElemID.end() ) {
+ int elemIDToRemove = *idIt;
+ const SMDS_MeshElement* elemToRemove = aMesh->FindElement(elemIDToRemove);
+ // add the kept element in groups of removed one (PAL15188)
+ AddToSameGroups( elemToKeep, elemToRemove, aMesh );
+ rmElemIds.push_back( elemIDToRemove );
+ ++idIt;
+ }
+ ++groupsIt;
+ }
+
+ Remove( rmElemIds, false );
+}
+
+//=======================================================================
+//function : MergeEqualElements
+//purpose : Remove all but one of elements built on the same nodes.
+//=======================================================================
+
+void SMESH_MeshEditor::MergeEqualElements()
+{
+ TIDSortedElemSet aMeshElements; /* empty input ==
+ to merge equal elements in the whole mesh */
+ TListOfListOfElementsID aGroupsOfElementsID;
+ FindEqualElements( aMeshElements, aGroupsOfElementsID );
+ MergeElements( aGroupsOfElementsID );
+}
+
+//=======================================================================
+//function : findAdjacentFace
+//purpose :
+//=======================================================================
+
+static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ const SMDS_MeshElement* elem)
+{
+ TIDSortedElemSet elemSet, avoidSet;
+ if ( elem )
+ avoidSet.insert ( elem );
+ return SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet );
+}
+
+//=======================================================================
+//function : findSegment
+//purpose : Return a mesh segment by two nodes one of which can be medium
+//=======================================================================
+
+static const SMDS_MeshElement* findSegment(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2)
+{
+ SMDS_ElemIteratorPtr it = n1->GetInverseElementIterator( SMDSAbs_Edge );
+ while ( it->more() )
+ {
+ const SMDS_MeshElement* seg = it->next();
+ if ( seg->GetNodeIndex( n2 ) >= 0 )
+ return seg;
+ }
+ return 0;
+}
+
+//=======================================================================
+//function : FindFreeBorder
+//purpose :
+//=======================================================================
+
+#define ControlFreeBorder SMESH::Controls::FreeEdges::IsFreeEdge
+
+bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode* theFirstNode,
+ const SMDS_MeshNode* theSecondNode,
+ const SMDS_MeshNode* theLastNode,
+ list< const SMDS_MeshNode* > & theNodes,
+ list< const SMDS_MeshElement* >& theFaces)
+{
+ if ( !theFirstNode || !theSecondNode )
+ return false;
+ // find border face between theFirstNode and theSecondNode
+ const SMDS_MeshElement* curElem = findAdjacentFace( theFirstNode, theSecondNode, 0 );
+ if ( !curElem )
+ return false;
+
+ theFaces.push_back( curElem );
+ theNodes.push_back( theFirstNode );
+ theNodes.push_back( theSecondNode );
+
+ const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
+ //TIDSortedElemSet foundElems;
+ bool needTheLast = ( theLastNode != 0 );
+
+ vector<const SMDS_MeshNode*> nodes;
+
+ while ( nStart != theLastNode ) {
+ if ( nStart == theFirstNode )
+ return !needTheLast;
+
+ // find all free border faces sharing nStart
+
+ list< const SMDS_MeshElement* > curElemList;
+ list< const SMDS_MeshNode* > nStartList;
+ SMDS_ElemIteratorPtr invElemIt = nStart->GetInverseElementIterator(SMDSAbs_Face);
+ while ( invElemIt->more() ) {
+ const SMDS_MeshElement* e = invElemIt->next();
+ //if ( e == curElem || foundElems.insert( e ).second ) // e can encounter twice in border
+ {
+ // get nodes
+ nodes.assign( SMDS_MeshElement::iterator( e->interlacedNodesIterator() ),
+ SMDS_MeshElement::iterator() );
+ nodes.push_back( nodes[ 0 ]);
+
+ // check 2 links
+ int iNode = 0, nbNodes = nodes.size() - 1;
+ for ( iNode = 0; iNode < nbNodes; iNode++ )
+ if ((( nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
+ ( nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
+ ( ControlFreeBorder( &nodes[ iNode ], e->GetID() )))
+ {
+ nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart )]);
+ curElemList.push_back( e );
+ }
+ }
+ }
+ // analyse the found
+
+ int nbNewBorders = curElemList.size();
+ if ( nbNewBorders == 0 ) {
+ // no free border furthermore
+ return !needTheLast;
+ }
+ else if ( nbNewBorders == 1 ) {
+ // one more element found
+ nIgnore = nStart;
+ nStart = nStartList.front();
+ curElem = curElemList.front();
+ theFaces.push_back( curElem );
+ theNodes.push_back( nStart );
+ }
+ else {
+ // several continuations found
+ list< const SMDS_MeshElement* >::iterator curElemIt;
+ list< const SMDS_MeshNode* >::iterator nStartIt;
+ // check if one of them reached the last node
+ if ( needTheLast ) {
+ for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
+ curElemIt!= curElemList.end();
+ curElemIt++, nStartIt++ )
+ if ( *nStartIt == theLastNode ) {
+ theFaces.push_back( *curElemIt );
+ theNodes.push_back( *nStartIt );
+ return true;
+ }
+ }
+ // find the best free border by the continuations
+ list<const SMDS_MeshNode*> contNodes[ 2 ], *cNL;
+ list<const SMDS_MeshElement*> contFaces[ 2 ], *cFL;
+ for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
+ curElemIt!= curElemList.end();
+ curElemIt++, nStartIt++ )
+ {
+ cNL = & contNodes[ contNodes[0].empty() ? 0 : 1 ];
+ cFL = & contFaces[ contFaces[0].empty() ? 0 : 1 ];
+ // find one more free border
+ if ( ! SMESH_MeshEditor::FindFreeBorder( nStart, *nStartIt, theLastNode, *cNL, *cFL )) {
+ cNL->clear();
+ cFL->clear();
+ }
+ else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
+ // choice: clear a worse one
+ int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
+ int iWorse = ( needTheLast ? 1 - iLongest : iLongest );
+ contNodes[ iWorse ].clear();
+ contFaces[ iWorse ].clear();
+ }
+ }
+ if ( contNodes[0].empty() && contNodes[1].empty() )
+ return false;
+
+ // push_back the best free border
+ cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
+ cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
+ //theNodes.pop_back(); // remove nIgnore
+ theNodes.pop_back(); // remove nStart
+ //theFaces.pop_back(); // remove curElem
+ theNodes.splice( theNodes.end(), *cNL );
+ theFaces.splice( theFaces.end(), *cFL );
+ return true;
+
+ } // several continuations found
+ } // while ( nStart != theLastNode )
+
+ return true;
+}
+
+//=======================================================================
+//function : CheckFreeBorderNodes
+//purpose : Return true if the tree nodes are on a free border
+//=======================================================================
+
+bool SMESH_MeshEditor::CheckFreeBorderNodes(const SMDS_MeshNode* theNode1,
+ const SMDS_MeshNode* theNode2,
+ const SMDS_MeshNode* theNode3)
+{
+ list< const SMDS_MeshNode* > nodes;
+ list< const SMDS_MeshElement* > faces;
+ return FindFreeBorder( theNode1, theNode2, theNode3, nodes, faces);
+}
+
+//=======================================================================
+//function : SewFreeBorder
+//purpose :
+//warning : for border-to-side sewing theSideSecondNode is considered as
+// the last side node and theSideThirdNode is not used
+//=======================================================================
+
+SMESH_MeshEditor::Sew_Error
+SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
+ const SMDS_MeshNode* theBordSecondNode,
+ const SMDS_MeshNode* theBordLastNode,
+ const SMDS_MeshNode* theSideFirstNode,
+ const SMDS_MeshNode* theSideSecondNode,
+ const SMDS_MeshNode* theSideThirdNode,
+ const bool theSideIsFreeBorder,
+ const bool toCreatePolygons,
+ const bool toCreatePolyedrs)
+{
+ ClearLastCreated();
+
+ Sew_Error aResult = SEW_OK;
+
+ // ====================================
+ // find side nodes and elements
+ // ====================================
+
+ list< const SMDS_MeshNode* > nSide[ 2 ];
+ list< const SMDS_MeshElement* > eSide[ 2 ];
+ list< const SMDS_MeshNode* >::iterator nIt[ 2 ];
+ list< const SMDS_MeshElement* >::iterator eIt[ 2 ];
+
+ // Free border 1
+ // --------------
+ if (!FindFreeBorder(theBordFirstNode,theBordSecondNode,theBordLastNode,
+ nSide[0], eSide[0])) {
+ MESSAGE(" Free Border 1 not found " );
+ aResult = SEW_BORDER1_NOT_FOUND;
+ }
+ if (theSideIsFreeBorder) {
+ // Free border 2
+ // --------------
+ if (!FindFreeBorder(theSideFirstNode, theSideSecondNode, theSideThirdNode,
+ nSide[1], eSide[1])) {
+ MESSAGE(" Free Border 2 not found " );
+ aResult = ( aResult != SEW_OK ? SEW_BOTH_BORDERS_NOT_FOUND : SEW_BORDER2_NOT_FOUND );
+ }
+ }
+ if ( aResult != SEW_OK )
+ return aResult;
+
+ if (!theSideIsFreeBorder) {
+ // Side 2
+ // --------------
+
+ // -------------------------------------------------------------------------
+ // Algo:
+ // 1. If nodes to merge are not coincident, move nodes of the free border
+ // from the coord sys defined by the direction from the first to last
+ // nodes of the border to the correspondent sys of the side 2
+ // 2. On the side 2, find the links most co-directed with the correspondent
+ // links of the free border
+ // -------------------------------------------------------------------------
+
+ // 1. Since sewing may break if there are volumes to split on the side 2,
+ // we won't move nodes but just compute new coordinates for them
+ typedef map<const SMDS_MeshNode*, gp_XYZ> TNodeXYZMap;
+ TNodeXYZMap nBordXYZ;
+ list< const SMDS_MeshNode* >& bordNodes = nSide[ 0 ];
+ list< const SMDS_MeshNode* >::iterator nBordIt;
+
+ gp_XYZ Pb1( theBordFirstNode->X(), theBordFirstNode->Y(), theBordFirstNode->Z() );
+ gp_XYZ Pb2( theBordLastNode->X(), theBordLastNode->Y(), theBordLastNode->Z() );
+ gp_XYZ Ps1( theSideFirstNode->X(), theSideFirstNode->Y(), theSideFirstNode->Z() );
+ gp_XYZ Ps2( theSideSecondNode->X(), theSideSecondNode->Y(), theSideSecondNode->Z() );
+ double tol2 = 1.e-8;
+ gp_Vec Vbs1( Pb1 - Ps1 ),Vbs2( Pb2 - Ps2 );
+ if ( Vbs1.SquareMagnitude() > tol2 || Vbs2.SquareMagnitude() > tol2 ) {
+ // Need node movement.
+
+ // find X and Z axes to create trsf
+ gp_Vec Zb( Pb1 - Pb2 ), Zs( Ps1 - Ps2 );
+ gp_Vec X = Zs ^ Zb;
+ if ( X.SquareMagnitude() <= gp::Resolution() * gp::Resolution() )
+ // Zb || Zs
+ X = gp_Ax2( gp::Origin(), Zb ).XDirection();
+
+ // coord systems
+ gp_Ax3 toBordAx( Pb1, Zb, X );
+ gp_Ax3 fromSideAx( Ps1, Zs, X );
+ gp_Ax3 toGlobalAx( gp::Origin(), gp::DZ(), gp::DX() );
+ // set trsf
+ gp_Trsf toBordSys, fromSide2Sys;
+ toBordSys.SetTransformation( toBordAx );
+ fromSide2Sys.SetTransformation( fromSideAx, toGlobalAx );
+ fromSide2Sys.SetScaleFactor( Zs.Magnitude() / Zb.Magnitude() );
+
+ // move
+ for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
+ const SMDS_MeshNode* n = *nBordIt;
+ gp_XYZ xyz( n->X(),n->Y(),n->Z() );
+ toBordSys.Transforms( xyz );
+ fromSide2Sys.Transforms( xyz );
+ nBordXYZ.insert( TNodeXYZMap::value_type( n, xyz ));
+ }
+ }
+ else {
+ // just insert nodes XYZ in the nBordXYZ map
+ for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
+ const SMDS_MeshNode* n = *nBordIt;
+ nBordXYZ.insert( TNodeXYZMap::value_type( n, gp_XYZ( n->X(),n->Y(),n->Z() )));
+ }
+ }
+
+ // 2. On the side 2, find the links most co-directed with the correspondent
+ // links of the free border
+
+ list< const SMDS_MeshElement* >& sideElems = eSide[ 1 ];
+ list< const SMDS_MeshNode* >& sideNodes = nSide[ 1 ];
+ sideNodes.push_back( theSideFirstNode );
+
+ bool hasVolumes = false;
+ LinkID_Gen aLinkID_Gen( GetMeshDS() );
+ set<long> foundSideLinkIDs, checkedLinkIDs;
+ SMDS_VolumeTool volume;
+ //const SMDS_MeshNode* faceNodes[ 4 ];
+
+ const SMDS_MeshNode* sideNode;
+ const SMDS_MeshElement* sideElem = 0;
+ const SMDS_MeshNode* prevSideNode = theSideFirstNode;
+ const SMDS_MeshNode* prevBordNode = theBordFirstNode;
+ nBordIt = bordNodes.begin();
+ nBordIt++;
+ // border node position and border link direction to compare with
+ gp_XYZ bordPos = nBordXYZ[ *nBordIt ];
+ gp_XYZ bordDir = bordPos - nBordXYZ[ prevBordNode ];
+ // choose next side node by link direction or by closeness to
+ // the current border node:
+ bool searchByDir = ( *nBordIt != theBordLastNode );
+ do {
+ // find the next node on the Side 2
+ sideNode = 0;
+ double maxDot = -DBL_MAX, minDist = DBL_MAX;
+ long linkID;
+ checkedLinkIDs.clear();
+ gp_XYZ prevXYZ( prevSideNode->X(), prevSideNode->Y(), prevSideNode->Z() );
+
+ // loop on inverse elements of current node (prevSideNode) on the Side 2
+ SMDS_ElemIteratorPtr invElemIt = prevSideNode->GetInverseElementIterator();
+ while ( invElemIt->more() )
+ {
+ const SMDS_MeshElement* elem = invElemIt->next();
+ // prepare data for a loop on links coming to prevSideNode, of a face or a volume
+ int iPrevNode = 0, iNode = 0, nbNodes = elem->NbNodes();
+ vector< const SMDS_MeshNode* > faceNodes( nbNodes, (const SMDS_MeshNode*)0 );
+ bool isVolume = volume.Set( elem );
+ const SMDS_MeshNode** nodes = isVolume ? volume.GetNodes() : & faceNodes[0];
+ if ( isVolume ) // --volume
+ hasVolumes = true;
+ else if ( elem->GetType() == SMDSAbs_Face ) { // --face
+ // retrieve all face nodes and find iPrevNode - an index of the prevSideNode
+ SMDS_NodeIteratorPtr nIt = elem->interlacedNodesIterator();
+ while ( nIt->more() ) {
+ nodes[ iNode ] = cast2Node( nIt->next() );
+ if ( nodes[ iNode++ ] == prevSideNode )
+ iPrevNode = iNode - 1;
+ }
+ // there are 2 links to check
+ nbNodes = 2;
+ }
+ else // --edge
+ continue;
+ // loop on links, to be precise, on the second node of links
+ for ( iNode = 0; iNode < nbNodes; iNode++ ) {
+ const SMDS_MeshNode* n = nodes[ iNode ];
+ if ( isVolume ) {
+ if ( !volume.IsLinked( n, prevSideNode ))
+ continue;
+ }
+ else {
+ if ( iNode ) // a node before prevSideNode
+ n = nodes[ iPrevNode == 0 ? elem->NbNodes() - 1 : iPrevNode - 1 ];
+ else // a node after prevSideNode
+ n = nodes[ iPrevNode + 1 == elem->NbNodes() ? 0 : iPrevNode + 1 ];
+ }
+ // check if this link was already used
+ long iLink = aLinkID_Gen.GetLinkID( prevSideNode, n );
+ bool isJustChecked = !checkedLinkIDs.insert( iLink ).second;
+ if (!isJustChecked &&
+ foundSideLinkIDs.find( iLink ) == foundSideLinkIDs.end() )
+ {
+ // test a link geometrically
+ gp_XYZ nextXYZ ( n->X(), n->Y(), n->Z() );
+ bool linkIsBetter = false;
+ double dot = 0.0, dist = 0.0;
+ if ( searchByDir ) { // choose most co-directed link
+ dot = bordDir * ( nextXYZ - prevXYZ ).Normalized();
+ linkIsBetter = ( dot > maxDot );
+ }
+ else { // choose link with the node closest to bordPos
+ dist = ( nextXYZ - bordPos ).SquareModulus();
+ linkIsBetter = ( dist < minDist );
+ }
+ if ( linkIsBetter ) {
+ maxDot = dot;
+ minDist = dist;
+ linkID = iLink;
+ sideNode = n;
+ sideElem = elem;
+ }
+ }
+ }
+ } // loop on inverse elements of prevSideNode
+
+ if ( !sideNode ) {
+ MESSAGE(" Can't find path by links of the Side 2 ");
+ return SEW_BAD_SIDE_NODES;
+ }
+ sideNodes.push_back( sideNode );
+ sideElems.push_back( sideElem );
+ foundSideLinkIDs.insert ( linkID );
+ prevSideNode = sideNode;
+
+ if ( *nBordIt == theBordLastNode )
+ searchByDir = false;
+ else {
+ // find the next border link to compare with
+ gp_XYZ sidePos( sideNode->X(), sideNode->Y(), sideNode->Z() );
+ searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
+ // move to next border node if sideNode is before forward border node (bordPos)
+ while ( *nBordIt != theBordLastNode && !searchByDir ) {
+ prevBordNode = *nBordIt;
+ nBordIt++;
+ bordPos = nBordXYZ[ *nBordIt ];
+ bordDir = bordPos - nBordXYZ[ prevBordNode ];
+ searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
+ }
+ }
+ }
+ while ( sideNode != theSideSecondNode );
+
+ if ( hasVolumes && sideNodes.size () != bordNodes.size() && !toCreatePolyedrs) {
+ MESSAGE("VOLUME SPLITTING IS FORBIDDEN");
+ return SEW_VOLUMES_TO_SPLIT; // volume splitting is forbidden
+ }
+ } // end nodes search on the side 2
+
+ // ============================
+ // sew the border to the side 2
+ // ============================
+
+ int nbNodes[] = { (int)nSide[0].size(), (int)nSide[1].size() };
+ int maxNbNodes = Max( nbNodes[0], nbNodes[1] );
+
+ bool toMergeConformal = ( nbNodes[0] == nbNodes[1] );
+ if ( toMergeConformal && toCreatePolygons )
+ {
+ // do not merge quadrangles if polygons are OK (IPAL0052824)
+ eIt[0] = eSide[0].begin();
+ eIt[1] = eSide[1].begin();
+ bool allQuads[2] = { true, true };
+ for ( int iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
+ for ( ; allQuads[iBord] && eIt[iBord] != eSide[iBord].end(); ++eIt[iBord] )
+ allQuads[iBord] = ( (*eIt[iBord])->NbCornerNodes() == 4 );
+ }
+ toMergeConformal = ( !allQuads[0] && !allQuads[1] );
+ }
+
+ TListOfListOfNodes nodeGroupsToMerge;
+ if (( toMergeConformal ) ||
+ ( theSideIsFreeBorder && !theSideThirdNode )) {
+
+ // all nodes are to be merged
+
+ for (nIt[0] = nSide[0].begin(), nIt[1] = nSide[1].begin();
+ nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end();
+ nIt[0]++, nIt[1]++ )
+ {
+ nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
+ nodeGroupsToMerge.back().push_back( *nIt[1] ); // to keep
+ nodeGroupsToMerge.back().push_back( *nIt[0] ); // to remove
+ }
+ }
+ else {
+
+ // insert new nodes into the border and the side to get equal nb of segments
+
+ // get normalized parameters of nodes on the borders
+ vector< double > param[ 2 ];
+ param[0].resize( maxNbNodes );
+ param[1].resize( maxNbNodes );
+ int iNode, iBord;
+ for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
+ list< const SMDS_MeshNode* >& nodes = nSide[ iBord ];
+ list< const SMDS_MeshNode* >::iterator nIt = nodes.begin();
+ const SMDS_MeshNode* nPrev = *nIt;
+ double bordLength = 0;
+ for ( iNode = 0; nIt != nodes.end(); nIt++, iNode++ ) { // loop on border nodes
+ const SMDS_MeshNode* nCur = *nIt;
+ gp_XYZ segment (nCur->X() - nPrev->X(),
+ nCur->Y() - nPrev->Y(),
+ nCur->Z() - nPrev->Z());
+ double segmentLen = segment.Modulus();
+ bordLength += segmentLen;
+ param[ iBord ][ iNode ] = bordLength;
+ nPrev = nCur;
+ }
+ // normalize within [0,1]
+ for ( iNode = 0; iNode < nbNodes[ iBord ]; iNode++ ) {
+ param[ iBord ][ iNode ] /= bordLength;
+ }
+ }
+
+ // loop on border segments
+ const SMDS_MeshNode *nPrev[ 2 ] = { 0, 0 };
+ int i[ 2 ] = { 0, 0 };
+ nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
+ nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
+
+ // element can be split while iterating on border if it has two edges in the border
+ std::map< const SMDS_MeshElement* , const SMDS_MeshElement* > elemReplaceMap;
+ std::map< const SMDS_MeshElement* , const SMDS_MeshElement* >::iterator elemReplaceMapIt;
+
+ TElemOfNodeListMap insertMap;
+ TElemOfNodeListMap::iterator insertMapIt;
+ // insertMap is
+ // key: elem to insert nodes into
+ // value: 2 nodes to insert between + nodes to be inserted
+ do {
+ bool next[ 2 ] = { false, false };
+
+ // find min adjacent segment length after sewing
+ double nextParam = 10., prevParam = 0;
+ for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
+ if ( i[ iBord ] + 1 < nbNodes[ iBord ])
+ nextParam = Min( nextParam, param[iBord][ i[iBord] + 1 ]);
+ if ( i[ iBord ] > 0 )
+ prevParam = Max( prevParam, param[iBord][ i[iBord] - 1 ]);
+ }
+ double minParam = Min( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
+ double maxParam = Max( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
+ double minSegLen = Min( nextParam - minParam, maxParam - prevParam );
+
+ // choose to insert or to merge nodes
+ double du = param[ 1 ][ i[1] ] - param[ 0 ][ i[0] ];
+ if ( Abs( du ) <= minSegLen * 0.2 ) {
+ // merge
+ // ------
+ nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
+ const SMDS_MeshNode* n0 = *nIt[0];
+ const SMDS_MeshNode* n1 = *nIt[1];
+ nodeGroupsToMerge.back().push_back( n1 );
+ nodeGroupsToMerge.back().push_back( n0 );
+ // position of node of the border changes due to merge
+ param[ 0 ][ i[0] ] += du;
+ // move n1 for the sake of elem shape evaluation during insertion.
+ // n1 will be removed by MergeNodes() anyway
+ const_cast<SMDS_MeshNode*>( n0 )->setXYZ( n1->X(), n1->Y(), n1->Z() );
+ next[0] = next[1] = true;
+ }
+ else {
+ // insert
+ // ------
+ int intoBord = ( du < 0 ) ? 0 : 1;
+ const SMDS_MeshElement* elem = *eIt [ intoBord ];
+ const SMDS_MeshNode* n1 = nPrev[ intoBord ];
+ const SMDS_MeshNode* n2 = *nIt [ intoBord ];
+ const SMDS_MeshNode* nIns = *nIt [ 1 - intoBord ];
+ if ( intoBord == 1 ) {
+ // move node of the border to be on a link of elem of the side
+ SMESH_NodeXYZ p1( n1 ), p2( n2 );
+ double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]);
+ gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio;
+ GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() );
+ }
+ elemReplaceMapIt = elemReplaceMap.find( elem );
+ if ( elemReplaceMapIt != elemReplaceMap.end() )
+ elem = elemReplaceMapIt->second;
+
+ insertMapIt = insertMap.find( elem );
+ bool notFound = ( insertMapIt == insertMap.end() );
+ bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 );
+ if ( otherLink ) {
+ // insert into another link of the same element:
+ // 1. perform insertion into the other link of the elem
+ list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
+ const SMDS_MeshNode* n12 = nodeList.front(); nodeList.pop_front();
+ const SMDS_MeshNode* n22 = nodeList.front(); nodeList.pop_front();
+ InsertNodesIntoLink( elem, n12, n22, nodeList, toCreatePolygons );
+ // 2. perform insertion into the link of adjacent faces
+ while ( const SMDS_MeshElement* adjElem = findAdjacentFace( n12, n22, elem )) {
+ InsertNodesIntoLink( adjElem, n12, n22, nodeList, toCreatePolygons );
+ }
+ while ( const SMDS_MeshElement* seg = findSegment( n12, n22 )) {
+ InsertNodesIntoLink( seg, n12, n22, nodeList );
+ }
+ if (toCreatePolyedrs) {
+ // perform insertion into the links of adjacent volumes
+ UpdateVolumes(n12, n22, nodeList);
+ }
+ // 3. find an element appeared on n1 and n2 after the insertion
+ insertMap.erase( insertMapIt );
+ const SMDS_MeshElement* elem2 = findAdjacentFace( n1, n2, 0 );
+ elemReplaceMap.insert( std::make_pair( elem, elem2 ));
+ elem = elem2;
+ }
+ if ( notFound || otherLink ) {
+ // add element and nodes of the side into the insertMap
+ insertMapIt = insertMap.insert( make_pair( elem, list<const SMDS_MeshNode*>() )).first;
+ (*insertMapIt).second.push_back( n1 );
+ (*insertMapIt).second.push_back( n2 );
+ }
+ // add node to be inserted into elem
+ (*insertMapIt).second.push_back( nIns );
+ next[ 1 - intoBord ] = true;
+ }
+
+ // go to the next segment
+ for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
+ if ( next[ iBord ] ) {
+ if ( i[ iBord ] != 0 && eIt[ iBord ] != eSide[ iBord ].end())
+ eIt[ iBord ]++;
+ nPrev[ iBord ] = *nIt[ iBord ];
+ nIt[ iBord ]++; i[ iBord ]++;
+ }
+ }
+ }
+ while ( nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end());
+
+ // perform insertion of nodes into elements
+
+ for (insertMapIt = insertMap.begin();
+ insertMapIt != insertMap.end();
+ insertMapIt++ )
+ {
+ const SMDS_MeshElement* elem = (*insertMapIt).first;
+ list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
+ if ( nodeList.size() < 3 ) continue;
+ const SMDS_MeshNode* n1 = nodeList.front(); nodeList.pop_front();
+ const SMDS_MeshNode* n2 = nodeList.front(); nodeList.pop_front();
+
+ InsertNodesIntoLink( elem, n1, n2, nodeList, toCreatePolygons );
+
+ while ( const SMDS_MeshElement* seg = findSegment( n1, n2 )) {
+ InsertNodesIntoLink( seg, n1, n2, nodeList );
+ }
+
+ if ( !theSideIsFreeBorder ) {
+ // look for and insert nodes into the faces adjacent to elem
+ while ( const SMDS_MeshElement* adjElem = findAdjacentFace( n1, n2, elem )) {
+ InsertNodesIntoLink( adjElem, n1, n2, nodeList, toCreatePolygons );
+ }
+ }
+ if (toCreatePolyedrs) {
+ // perform insertion into the links of adjacent volumes
+ UpdateVolumes(n1, n2, nodeList);
+ }
+ }
+ } // end: insert new nodes
+
+ MergeNodes ( nodeGroupsToMerge );
+
+
+ // Remove coincident segments
+
+ // get new segments
+ TIDSortedElemSet segments;
+ SMESH_SequenceOfElemPtr newFaces;
+ for ( size_t i = 0; i < myLastCreatedElems.size(); ++i )
+ {
+ if ( !myLastCreatedElems[i] ) continue;
+ if ( myLastCreatedElems[i]->GetType() == SMDSAbs_Edge )
+ segments.insert( segments.end(), myLastCreatedElems[i] );
+ else
+ newFaces.push_back( myLastCreatedElems[i] );
+ }
+ // get segments adjacent to merged nodes
+ TListOfListOfNodes::iterator groupIt = nodeGroupsToMerge.begin();
+ for ( ; groupIt != nodeGroupsToMerge.end(); groupIt++ )
+ {
+ const list<const SMDS_MeshNode*>& nodes = *groupIt;
+ if ( nodes.front()->IsNull() ) continue;
+ SMDS_ElemIteratorPtr segIt = nodes.front()->GetInverseElementIterator( SMDSAbs_Edge );
+ while ( segIt->more() )
+ segments.insert( segIt->next() );
+ }
+
+ // find coincident
+ TListOfListOfElementsID equalGroups;
+ if ( !segments.empty() )
+ FindEqualElements( segments, equalGroups );
+ if ( !equalGroups.empty() )
+ {
+ // remove from segments those that will be removed
+ TListOfListOfElementsID::iterator itGroups = equalGroups.begin();
+ for ( ; itGroups != equalGroups.end(); ++itGroups )
+ {
+ list< int >& group = *itGroups;
+ list< int >::iterator id = group.begin();
+ for ( ++id; id != group.end(); ++id )
+ if ( const SMDS_MeshElement* seg = GetMeshDS()->FindElement( *id ))
+ segments.erase( seg );
+ }
+ // remove equal segments
+ MergeElements( equalGroups );
+
+ // restore myLastCreatedElems
+ myLastCreatedElems = newFaces;
+ TIDSortedElemSet::iterator seg = segments.begin();
+ for ( ; seg != segments.end(); ++seg )
+ myLastCreatedElems.push_back( *seg );
+ }
+
+ return aResult;
+}
+
+//=======================================================================
+//function : InsertNodesIntoLink
+//purpose : insert theNodesToInsert into theElement between theBetweenNode1
+// and theBetweenNode2 and split theElement
+//=======================================================================
+
+void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theElement,
+ const SMDS_MeshNode* theBetweenNode1,
+ const SMDS_MeshNode* theBetweenNode2,
+ list<const SMDS_MeshNode*>& theNodesToInsert,
+ const bool toCreatePoly)
+{
+ if ( !theElement ) return;
+
+ SMESHDS_Mesh *aMesh = GetMeshDS();
+ vector<const SMDS_MeshElement*> newElems;
+
+ if ( theElement->GetType() == SMDSAbs_Edge )
+ {
+ theNodesToInsert.push_front( theBetweenNode1 );
+ theNodesToInsert.push_back ( theBetweenNode2 );
+ list<const SMDS_MeshNode*>::iterator n = theNodesToInsert.begin();
+ const SMDS_MeshNode* n1 = *n;
+ for ( ++n; n != theNodesToInsert.end(); ++n )
+ {
+ const SMDS_MeshNode* n2 = *n;
+ if ( const SMDS_MeshElement* seg = aMesh->FindEdge( n1, n2 ))
+ AddToSameGroups( seg, theElement, aMesh );
+ else
+ newElems.push_back( aMesh->AddEdge ( n1, n2 ));
+ n1 = n2;
+ }
+ theNodesToInsert.pop_front();
+ theNodesToInsert.pop_back();
+
+ if ( theElement->IsQuadratic() ) // add a not split part
+ {
+ vector<const SMDS_MeshNode*> nodes( theElement->begin_nodes(),
+ theElement->end_nodes() );
+ int iOther = 0, nbN = nodes.size();
+ for ( ; iOther < nbN; ++iOther )
+ if ( nodes[iOther] != theBetweenNode1 &&
+ nodes[iOther] != theBetweenNode2 )
+ break;
+ if ( iOther == 0 )
+ {
+ if ( const SMDS_MeshElement* seg = aMesh->FindEdge( nodes[0], nodes[1] ))
+ AddToSameGroups( seg, theElement, aMesh );
+ else
+ newElems.push_back( aMesh->AddEdge ( nodes[0], nodes[1] ));
+ }
+ else if ( iOther == 2 )
+ {
+ if ( const SMDS_MeshElement* seg = aMesh->FindEdge( nodes[1], nodes[2] ))
+ AddToSameGroups( seg, theElement, aMesh );
+ else
+ newElems.push_back( aMesh->AddEdge ( nodes[1], nodes[2] ));
+ }
+ }
+ // treat new elements
+ for ( size_t i = 0; i < newElems.size(); ++i )
+ if ( newElems[i] )
+ {
+ aMesh->SetMeshElementOnShape( newElems[i], theElement->getshapeId() );
+ myLastCreatedElems.push_back( newElems[i] );
+ }
+ ReplaceElemInGroups( theElement, newElems, aMesh );
+ aMesh->RemoveElement( theElement );
+ return;
+
+ } // if ( theElement->GetType() == SMDSAbs_Edge )
+
+ const SMDS_MeshElement* theFace = theElement;
+ if ( theFace->GetType() != SMDSAbs_Face ) return;
+
+ // find indices of 2 link nodes and of the rest nodes
+ int iNode = 0, il1, il2, i3, i4;
+ il1 = il2 = i3 = i4 = -1;
+ vector<const SMDS_MeshNode*> nodes( theFace->NbNodes() );
+
+ SMDS_NodeIteratorPtr nodeIt = theFace->interlacedNodesIterator();
+ while ( nodeIt->more() ) {
+ const SMDS_MeshNode* n = nodeIt->next();
+ if ( n == theBetweenNode1 )
+ il1 = iNode;
+ else if ( n == theBetweenNode2 )
+ il2 = iNode;
+ else if ( i3 < 0 )
+ i3 = iNode;
+ else
+ i4 = iNode;
+ nodes[ iNode++ ] = n;
+ }
+ if ( il1 < 0 || il2 < 0 || i3 < 0 )
+ return ;
+
+ // arrange link nodes to go one after another regarding the face orientation
+ bool reverse = ( Abs( il2 - il1 ) == 1 ? il2 < il1 : il1 < il2 );
+ list<const SMDS_MeshNode *> aNodesToInsert = theNodesToInsert;
+ if ( reverse ) {
+ iNode = il1;
+ il1 = il2;
+ il2 = iNode;
+ aNodesToInsert.reverse();
+ }
+ // check that not link nodes of a quadrangles are in good order
+ int nbFaceNodes = theFace->NbNodes();
+ if ( nbFaceNodes == 4 && i4 - i3 != 1 ) {
+ iNode = i3;
+ i3 = i4;
+ i4 = iNode;
+ }
+
+ if (toCreatePoly || theFace->IsPoly()) {
+
+ iNode = 0;
+ vector<const SMDS_MeshNode *> poly_nodes (nbFaceNodes + aNodesToInsert.size());
+
+ // add nodes of face up to first node of link
+ bool isFLN = false;
+ SMDS_NodeIteratorPtr nodeIt = theFace->interlacedNodesIterator();
+ while ( nodeIt->more() && !isFLN ) {
+ const SMDS_MeshNode* n = nodeIt->next();
+ poly_nodes[iNode++] = n;
+ isFLN = ( n == nodes[il1] );
+ }
+ // add nodes to insert
+ list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
+ for (; nIt != aNodesToInsert.end(); nIt++) {
+ poly_nodes[iNode++] = *nIt;
+ }
+ // add nodes of face starting from last node of link
+ while ( nodeIt->more() ) {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ poly_nodes[iNode++] = n;
+ }
+
+ // make a new face
+ newElems.push_back( aMesh->AddPolygonalFace( poly_nodes ));
+ }
+
+ else if ( !theFace->IsQuadratic() )
+ {
+ // put aNodesToInsert between theBetweenNode1 and theBetweenNode2
+ int nbLinkNodes = 2 + aNodesToInsert.size();
+ //const SMDS_MeshNode* linkNodes[ nbLinkNodes ];
+ vector<const SMDS_MeshNode*> linkNodes( nbLinkNodes );
+ linkNodes[ 0 ] = nodes[ il1 ];
+ linkNodes[ nbLinkNodes - 1 ] = nodes[ il2 ];
+ list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
+ for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
+ linkNodes[ iNode++ ] = *nIt;
+ }
+ // decide how to split a quadrangle: compare possible variants
+ // and choose which of splits to be a quadrangle
+ int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad = 0;
+ if ( nbFaceNodes == 3 ) {
+ iBestQuad = nbSplits;
+ i4 = i3;
+ }
+ else if ( nbFaceNodes == 4 ) {
+ SMESH::Controls::NumericalFunctorPtr aCrit( new SMESH::Controls::AspectRatio);
+ double aBestRate = DBL_MAX;
+ for ( int iQuad = 0; iQuad < nbSplits; iQuad++ ) {
+ i1 = 0; i2 = 1;
+ double aBadRate = 0;
+ // evaluate elements quality
+ for ( iSplit = 0; iSplit < nbSplits; iSplit++ ) {
+ if ( iSplit == iQuad ) {
+ SMDS_FaceOfNodes quad (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ i3 ],
+ nodes[ i4 ]);
+ aBadRate += getBadRate( &quad, aCrit );
+ }
+ else {
+ SMDS_FaceOfNodes tria (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ iSplit < iQuad ? i4 : i3 ]);
+ aBadRate += getBadRate( &tria, aCrit );
+ }
+ }
+ // choice
+ if ( aBadRate < aBestRate ) {
+ iBestQuad = iQuad;
+ aBestRate = aBadRate;
+ }
+ }
+ }
+
+ // create new elements
+ i1 = 0; i2 = 1;
+ for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ )
+ {
+ if ( iSplit == iBestQuad )
+ newElems.push_back( aMesh->AddFace (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ i3 ],
+ nodes[ i4 ]));
+ else
+ newElems.push_back( aMesh->AddFace (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ iSplit < iBestQuad ? i4 : i3 ]));
+ }
+
+ const SMDS_MeshNode* newNodes[ 4 ];
+ newNodes[ 0 ] = linkNodes[ i1 ];
+ newNodes[ 1 ] = linkNodes[ i2 ];
+ newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
+ newNodes[ 3 ] = nodes[ i4 ];
+ if (iSplit == iBestQuad)
+ newElems.push_back( aMesh->AddFace( newNodes[0], newNodes[1], newNodes[2], newNodes[3] ));
+ else
+ newElems.push_back( aMesh->AddFace( newNodes[0], newNodes[1], newNodes[2] ));
+
+ } // end if(!theFace->IsQuadratic())
+
+ else { // theFace is quadratic
+ // we have to split theFace on simple triangles and one simple quadrangle
+ int tmp = il1/2;
+ int nbshift = tmp*2;
+ // shift nodes in nodes[] by nbshift
+ int i,j;
+ for(i=0; i<nbshift; i++) {
+ const SMDS_MeshNode* n = nodes[0];
+ for(j=0; j<nbFaceNodes-1; j++) {
+ nodes[j] = nodes[j+1];
+ }
+ nodes[nbFaceNodes-1] = n;
+ }
+ il1 = il1 - nbshift;
+ // now have to insert nodes between n0 and n1 or n1 and n2 (see below)
+ // n0 n1 n2 n0 n1 n2
+ // +-----+-----+ +-----+-----+
+ // \ / | |
+ // \ / | |
+ // n5+ +n3 n7+ +n3
+ // \ / | |
+ // \ / | |
+ // + +-----+-----+
+ // n4 n6 n5 n4
+
+ // create new elements
+ int n1,n2,n3;
+ if ( nbFaceNodes == 6 ) { // quadratic triangle
+ newElems.push_back( aMesh->AddFace( nodes[3], nodes[4], nodes[5] ));
+ if ( theFace->IsMediumNode(nodes[il1]) ) {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[0], nodes[1], nodes[3], nodes[5] ));
+ n1 = 1;
+ n2 = 2;
+ n3 = 3;
+ }
+ else {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[1], nodes[2], nodes[3], nodes[5] ));
+ n1 = 0;
+ n2 = 1;
+ n3 = 5;
+ }
+ }
+ else { // nbFaceNodes==8 - quadratic quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[3], nodes[4], nodes[5] ));
+ newElems.push_back( aMesh->AddFace( nodes[5], nodes[6], nodes[7] ));
+ newElems.push_back( aMesh->AddFace( nodes[5], nodes[7], nodes[3] ));
+ if ( theFace->IsMediumNode( nodes[ il1 ])) {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[0], nodes[1], nodes[3], nodes[7] ));
+ n1 = 1;
+ n2 = 2;
+ n3 = 3;
+ }
+ else {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[1], nodes[2], nodes[3], nodes[7] ));
+ n1 = 0;
+ n2 = 1;
+ n3 = 7;
+ }
+ }
+ // create needed triangles using n1,n2,n3 and inserted nodes
+ int nbn = 2 + aNodesToInsert.size();
+ vector<const SMDS_MeshNode*> aNodes(nbn);
+ aNodes[0 ] = nodes[n1];
+ aNodes[nbn-1] = nodes[n2];
+ list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
+ for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
+ aNodes[iNode++] = *nIt;
+ }
+ for ( i = 1; i < nbn; i++ )
+ newElems.push_back( aMesh->AddFace( aNodes[i-1], aNodes[i], nodes[n3] ));
+ }
+
+ // remove the old face
+ for ( size_t i = 0; i < newElems.size(); ++i )
+ if ( newElems[i] )
+ {
+ aMesh->SetMeshElementOnShape( newElems[i], theFace->getshapeId() );
+ myLastCreatedElems.push_back( newElems[i] );
+ }
+ ReplaceElemInGroups( theFace, newElems, aMesh );
+ aMesh->RemoveElement(theFace);
+
+} // InsertNodesIntoLink()
+
+//=======================================================================
+//function : UpdateVolumes
+//purpose :
+//=======================================================================
+
+void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode* theBetweenNode1,
+ const SMDS_MeshNode* theBetweenNode2,
+ list<const SMDS_MeshNode*>& theNodesToInsert)
+{
+ ClearLastCreated();
+
+ SMDS_ElemIteratorPtr invElemIt = theBetweenNode1->GetInverseElementIterator(SMDSAbs_Volume);
+ while (invElemIt->more()) { // loop on inverse elements of theBetweenNode1
+ const SMDS_MeshElement* elem = invElemIt->next();
+
+ // check, if current volume has link theBetweenNode1 - theBetweenNode2
+ SMDS_VolumeTool aVolume (elem);
+ if (!aVolume.IsLinked(theBetweenNode1, theBetweenNode2))
+ continue;
+
+ // insert new nodes in all faces of the volume, sharing link theBetweenNode1 - theBetweenNode2
+ int iface, nbFaces = aVolume.NbFaces();
+ vector<const SMDS_MeshNode *> poly_nodes;
+ vector<int> quantities (nbFaces);
+
+ for (iface = 0; iface < nbFaces; iface++) {
+ int nbFaceNodes = aVolume.NbFaceNodes(iface), nbInserted = 0;
+ // faceNodes will contain (nbFaceNodes + 1) nodes, last = first
+ const SMDS_MeshNode** faceNodes = aVolume.GetFaceNodes(iface);
+
+ for (int inode = 0; inode < nbFaceNodes; inode++) {
+ poly_nodes.push_back(faceNodes[inode]);
+
+ if (nbInserted == 0) {
+ if (faceNodes[inode] == theBetweenNode1) {
+ if (faceNodes[inode + 1] == theBetweenNode2) {
+ nbInserted = theNodesToInsert.size();
+
+ // add nodes to insert
+ list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.begin();
+ for (; nIt != theNodesToInsert.end(); nIt++) {
+ poly_nodes.push_back(*nIt);
+ }
+ }
+ }
+ else if (faceNodes[inode] == theBetweenNode2) {
+ if (faceNodes[inode + 1] == theBetweenNode1) {
+ nbInserted = theNodesToInsert.size();
+
+ // add nodes to insert in reversed order
+ list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.end();
+ nIt--;
+ for (; nIt != theNodesToInsert.begin(); nIt--) {
+ poly_nodes.push_back(*nIt);
+ }
+ poly_nodes.push_back(*nIt);
+ }
+ }
+ else {
+ }
+ }
+ }
+ quantities[iface] = nbFaceNodes + nbInserted;
+ }
+
+ // Replace the volume
+ SMESHDS_Mesh *aMesh = GetMeshDS();
+
+ if ( SMDS_MeshElement* newElem = aMesh->AddPolyhedralVolume( poly_nodes, quantities ))
+ {
+ aMesh->SetMeshElementOnShape( newElem, elem->getshapeId() );
+ myLastCreatedElems.push_back( newElem );
+ ReplaceElemInGroups( elem, newElem, aMesh );
+ }
+ aMesh->RemoveElement( elem );
+ }
+}
+
+namespace
+{
+ //================================================================================
+ /*!
+ * \brief Transform any volume into data of SMDSEntity_Polyhedra
+ */
+ //================================================================================
+
+ void volumeToPolyhedron( const SMDS_MeshElement* elem,
+ vector<const SMDS_MeshNode *> & nodes,
+ vector<int> & nbNodeInFaces )
+ {
+ nodes.clear();
+ nbNodeInFaces.clear();
+ SMDS_VolumeTool vTool ( elem );
+ for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+ {
+ const SMDS_MeshNode** fNodes = vTool.GetFaceNodes( iF );
+ nodes.insert( nodes.end(), fNodes, fNodes + vTool.NbFaceNodes( iF ));
+ nbNodeInFaces.push_back( vTool.NbFaceNodes( iF ));
+ }
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Convert elements contained in a sub-mesh to quadratic
+ * \return int - nb of checked elements
+ */
+//=======================================================================
+
+int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm,
+ SMESH_MesherHelper& theHelper,
+ const bool theForce3d)
+{
+ //MESSAGE("convertElemToQuadratic");
+ int nbElem = 0;
+ if( !theSm ) return nbElem;
+
+ vector<int> nbNodeInFaces;
+ vector<const SMDS_MeshNode *> nodes;
+ SMDS_ElemIteratorPtr ElemItr = theSm->GetElements();
+ while(ElemItr->more())
+ {
+ nbElem++;
+ const SMDS_MeshElement* elem = ElemItr->next();
+ if( !elem ) continue;
+
+ // analyse a necessity of conversion
+ const SMDSAbs_ElementType aType = elem->GetType();
+ if ( aType < SMDSAbs_Edge || aType > SMDSAbs_Volume )
+ continue;
+ const SMDSAbs_EntityType aGeomType = elem->GetEntityType();
+ bool hasCentralNodes = false;
+ if ( elem->IsQuadratic() )
+ {
+ bool alreadyOK;
+ switch ( aGeomType ) {
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa:
+ case SMDSEntity_Quad_Penta:
+ alreadyOK = !theHelper.GetIsBiQuadratic(); break;
+
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ case SMDSEntity_TriQuad_Hexa:
+ case SMDSEntity_BiQuad_Penta:
+ alreadyOK = theHelper.GetIsBiQuadratic();
+ hasCentralNodes = true;
+ break;
+ default:
+ alreadyOK = true;
+ }
+ // take into account already present medium nodes
+ switch ( aType ) {
+ case SMDSAbs_Volume:
+ theHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( elem )); break;
+ case SMDSAbs_Face:
+ theHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( elem )); break;
+ case SMDSAbs_Edge:
+ theHelper.AddTLinks( static_cast< const SMDS_MeshEdge* >( elem )); break;
+ default:;
+ }
+ if ( alreadyOK )
+ continue;
+ }
+ // get elem data needed to re-create it
+ //
+ const int id = elem->GetID();
+ const int nbNodes = elem->NbCornerNodes();
+ nodes.assign(elem->begin_nodes(), elem->end_nodes());
+ if ( aGeomType == SMDSEntity_Polyhedra )
+ nbNodeInFaces = static_cast<const SMDS_MeshVolume* >( elem )->GetQuantities();
+ else if ( aGeomType == SMDSEntity_Hexagonal_Prism )
+ volumeToPolyhedron( elem, nodes, nbNodeInFaces );
+
+ // remove a linear element
+ GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
+
+ // remove central nodes of biquadratic elements (biquad->quad conversion)
+ if ( hasCentralNodes )
+ for ( size_t i = nbNodes * 2; i < nodes.size(); ++i )
+ if ( nodes[i]->NbInverseElements() == 0 )
+ GetMeshDS()->RemoveFreeNode( nodes[i], theSm, /*fromGroups=*/true );
+
+ const SMDS_MeshElement* NewElem = 0;
+
+ switch( aType )
+ {
+ case SMDSAbs_Edge :
+ {
+ NewElem = theHelper.AddEdge(nodes[0], nodes[1], id, theForce3d);
+ break;
+ }
+ case SMDSAbs_Face :
+ {
+ switch(nbNodes)
+ {
+ case 3:
+ NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
+ break;
+ case 4:
+ NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ default:
+ NewElem = theHelper.AddPolygonalFace(nodes, id, theForce3d);
+ }
+ break;
+ }
+ case SMDSAbs_Volume :
+ {
+ switch( aGeomType )
+ {
+ case SMDSEntity_Tetra:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ case SMDSEntity_Pyramid:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], id, theForce3d);
+ break;
+ case SMDSEntity_Penta:
+ case SMDSEntity_Quad_Penta:
+ case SMDSEntity_BiQuad_Penta:
+ 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;
+ case SMDSEntity_Hexagonal_Prism:
+ default:
+ NewElem = theHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
+ }
+ break;
+ }
+ default :
+ continue;
+ }
+ ReplaceElemInGroups( elem, NewElem, GetMeshDS());
+ if( NewElem && NewElem->getshapeId() < 1 )
+ theSm->AddElement( NewElem );
+ }
+ return nbElem;
+}
+//=======================================================================
+//function : ConvertToQuadratic
+//purpose :
+//=======================================================================
+
+void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theToBiQuad)
+{
+ //MESSAGE("ConvertToQuadratic "<< theForce3d << " " << theToBiQuad);
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+
+ SMESH_MesherHelper aHelper(*myMesh);
+
+ aHelper.SetIsQuadratic( true );
+ aHelper.SetIsBiQuadratic( theToBiQuad );
+ aHelper.SetElementsOnShape(true);
+ aHelper.ToFixNodeParameters( true );
+
+ // convert elements assigned to sub-meshes
+ int nbCheckedElems = 0;
+ if ( myMesh->HasShapeToMesh() )
+ {
+ if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
+ {
+ SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
+ while ( smIt->more() ) {
+ SMESH_subMesh* sm = smIt->next();
+ if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() ) {
+ aHelper.SetSubShape( sm->GetSubShape() );
+ nbCheckedElems += convertElemToQuadratic(smDS, aHelper, theForce3d);
+ }
+ }
+ }
+ }
+
+ // convert elements NOT assigned to sub-meshes
+ int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
+ if ( nbCheckedElems < totalNbElems ) // not all elements are in sub-meshes
+ {
+ aHelper.SetElementsOnShape(false);
+ SMESHDS_SubMesh *smDS = 0;
+
+ // convert edges
+ SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator();
+ while( aEdgeItr->more() )
+ {
+ const SMDS_MeshEdge* edge = aEdgeItr->next();
+ if ( !edge->IsQuadratic() )
+ {
+ int id = edge->GetID();
+ const SMDS_MeshNode* n1 = edge->GetNode(0);
+ const SMDS_MeshNode* n2 = edge->GetNode(1);
+
+ meshDS->RemoveFreeElement(edge, smDS, /*fromGroups=*/false);
+
+ const SMDS_MeshEdge* NewEdge = aHelper.AddEdge(n1, n2, id, theForce3d);
+ ReplaceElemInGroups( edge, NewEdge, GetMeshDS());
+ }
+ else
+ {
+ aHelper.AddTLinks( static_cast< const SMDS_MeshEdge* >( edge ));
+ }
+ }
+
+ // convert faces
+ SMDS_FaceIteratorPtr aFaceItr = meshDS->facesIterator();
+ while( aFaceItr->more() )
+ {
+ const SMDS_MeshFace* face = aFaceItr->next();
+ if ( !face ) continue;
+
+ const SMDSAbs_EntityType type = face->GetEntityType();
+ bool alreadyOK;
+ switch( type )
+ {
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ alreadyOK = !theToBiQuad;
+ aHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( face ));
+ break;
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ alreadyOK = theToBiQuad;
+ aHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( face ));
+ break;
+ default: alreadyOK = false;
+ }
+ if ( alreadyOK )
+ continue;
+
+ const int id = face->GetID();
+ vector<const SMDS_MeshNode *> nodes ( face->begin_nodes(), face->end_nodes());
+
+ meshDS->RemoveFreeElement(face, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshFace * NewFace = 0;
+ switch( type )
+ {
+ case SMDSEntity_Triangle:
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_BiQuad_Triangle:
+ NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
+ if ( nodes.size() == 7 && nodes[6]->NbInverseElements() == 0 ) // rm a central node
+ GetMeshDS()->RemoveFreeNode( nodes[6], /*sm=*/0, /*fromGroups=*/true );
+ break;
+
+ case SMDSEntity_Quadrangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ if ( nodes.size() == 9 && nodes[8]->NbInverseElements() == 0 ) // rm a central node
+ GetMeshDS()->RemoveFreeNode( nodes[8], /*sm=*/0, /*fromGroups=*/true );
+ break;
+
+ default:;
+ NewFace = aHelper.AddPolygonalFace(nodes, id, theForce3d);
+ }
+ ReplaceElemInGroups( face, NewFace, GetMeshDS());
+ }
+
+ // convert volumes
+ vector<int> nbNodeInFaces;
+ SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator();
+ while(aVolumeItr->more())
+ {
+ const SMDS_MeshVolume* volume = aVolumeItr->next();
+ if ( !volume ) continue;
+
+ const SMDSAbs_EntityType type = volume->GetEntityType();
+ if ( volume->IsQuadratic() )
+ {
+ bool alreadyOK;
+ switch ( type )
+ {
+ case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break;
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break;
+ case SMDSEntity_Quad_Penta: alreadyOK = !theToBiQuad; break;
+ case SMDSEntity_BiQuad_Penta: alreadyOK = theToBiQuad; break;
+ default: alreadyOK = true;
+ }
+ if ( alreadyOK )
+ {
+ aHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( volume ));
+ 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_MeshVolume* >(volume)->GetQuantities();
+ else if ( type == SMDSEntity_Hexagonal_Prism )
+ volumeToPolyhedron( volume, nodes, nbNodeInFaces );
+
+ meshDS->RemoveFreeElement(volume, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshVolume * NewVolume = 0;
+ switch ( type )
+ {
+ case SMDSEntity_Tetra:
+ 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);
+ for ( size_t i = 20; i < nodes.size(); ++i ) // rm central nodes
+ if ( nodes[i]->NbInverseElements() == 0 )
+ GetMeshDS()->RemoveFreeNode( nodes[i], /*sm=*/0, /*fromGroups=*/true );
+ break;
+ case SMDSEntity_Pyramid:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+ nodes[3], nodes[4], id, theForce3d);
+ break;
+ case SMDSEntity_Penta:
+ case SMDSEntity_Quad_Penta:
+ case SMDSEntity_BiQuad_Penta:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+ nodes[3], nodes[4], nodes[5], id, theForce3d);
+ for ( size_t i = 15; i < nodes.size(); ++i ) // rm central nodes
+ if ( nodes[i]->NbInverseElements() == 0 )
+ GetMeshDS()->RemoveFreeNode( nodes[i], /*sm=*/0, /*fromGroups=*/true );
+ break;
+ case SMDSEntity_Hexagonal_Prism:
+ default:
+ NewVolume = aHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
+ }
+ ReplaceElemInGroups(volume, NewVolume, meshDS);
+ }
+ }
+
+ 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(myError);
+ SMESH_MesherHelper( *myMesh ).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
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d,
+ TIDSortedElemSet& theElements,
+ const bool theToBiQuad)
+{
+ if ( theElements.empty() ) return;
+
+ // we believe that all theElements are of the same type
+ const SMDSAbs_ElementType elemType = (*theElements.begin())->GetType();
+
+ // get all nodes shared by theElements
+ TIDSortedNodeSet allNodes;
+ TIDSortedElemSet::iterator eIt = theElements.begin();
+ for ( ; eIt != theElements.end(); ++eIt )
+ allNodes.insert( (*eIt)->begin_nodes(), (*eIt)->end_nodes() );
+
+ // complete theElements with elements of lower dim whose all nodes are in allNodes
+
+ TIDSortedElemSet quadAdjacentElems [ SMDSAbs_NbElementTypes ]; // quadratic adjacent elements
+ TIDSortedElemSet checkedAdjacentElems [ SMDSAbs_NbElementTypes ];
+ TIDSortedNodeSet::iterator nIt = allNodes.begin();
+ for ( ; nIt != allNodes.end(); ++nIt )
+ {
+ const SMDS_MeshNode* n = *nIt;
+ SMDS_ElemIteratorPtr invIt = n->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* e = invIt->next();
+ const SMDSAbs_ElementType type = e->GetType();
+ if ( e->IsQuadratic() )
+ {
+ quadAdjacentElems[ type ].insert( e );
+
+ bool alreadyOK;
+ switch ( e->GetEntityType() ) {
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break;
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break;
+ default: alreadyOK = true;
+ }
+ if ( alreadyOK )
+ continue;
+ }
+ if ( type >= elemType )
+ continue; // same type or more complex linear element
+
+ if ( !checkedAdjacentElems[ type ].insert( e ).second )
+ continue; // e is already checked
+
+ // check nodes
+ bool allIn = true;
+ SMDS_NodeIteratorPtr nodeIt = e->nodeIterator();
+ while ( nodeIt->more() && allIn )
+ allIn = allNodes.count( nodeIt->next() );
+ if ( allIn )
+ theElements.insert(e );
+ }
+ }
+
+ SMESH_MesherHelper helper(*myMesh);
+ helper.SetIsQuadratic( true );
+ helper.SetIsBiQuadratic( theToBiQuad );
+
+ // add links of quadratic adjacent elements to the helper
+
+ if ( !quadAdjacentElems[SMDSAbs_Edge].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Edge].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Edge].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshEdge*> (*eIt) );
+ }
+ if ( !quadAdjacentElems[SMDSAbs_Face].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Face].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Face].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshFace*> (*eIt) );
+ }
+ if ( !quadAdjacentElems[SMDSAbs_Volume].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Volume].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Volume].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshVolume*> (*eIt) );
+ }
+
+ // make quadratic (or bi-tri-quadratic) elements instead of linear ones
+
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+ SMESHDS_SubMesh* smDS = 0;
+ for ( eIt = theElements.begin(); eIt != theElements.end(); ++eIt )
+ {
+ const SMDS_MeshElement* elem = *eIt;
+
+ bool alreadyOK;
+ int nbCentralNodes = 0;
+ switch ( elem->GetEntityType() ) {
+ // linear convertible
+ case SMDSEntity_Edge:
+ case SMDSEntity_Triangle:
+ case SMDSEntity_Quadrangle:
+ case SMDSEntity_Tetra:
+ case SMDSEntity_Pyramid:
+ case SMDSEntity_Hexa:
+ case SMDSEntity_Penta: alreadyOK = false; nbCentralNodes = 0; break;
+ // quadratic that can become bi-quadratic
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa: alreadyOK =!theToBiQuad; nbCentralNodes = 0; break;
+ // bi-quadratic
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle: alreadyOK = theToBiQuad; nbCentralNodes = 1; break;
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; nbCentralNodes = 7; break;
+ // the rest
+ 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());
+
+ helper.SetSubShape( elem->getshapeId() );
+
+ if ( !smDS || !smDS->Contains( elem ))
+ smDS = meshDS->MeshElements( elem->getshapeId() );
+ meshDS->RemoveFreeElement(elem, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshElement * newElem = 0;
+ switch( nbNodes )
+ {
+ 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
+ newElem = helper.AddFace (nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ case 8:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+ break;
+ case 3:
+ newElem = helper.AddFace (nodes[0], nodes[1], nodes[2], id, theForce3d);
+ break;
+ case 2:
+ newElem = helper.AddEdge(nodes[0], nodes[1], id, theForce3d);
+ break;
+ case 5:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], id, theForce3d);
+ break;
+ case 6:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], id, theForce3d);
+ break;
+ default:;
+ }
+ ReplaceElemInGroups( elem, newElem, meshDS);
+ if( newElem && smDS )
+ smDS->AddElement( newElem );
+
+ // remove central nodes
+ for ( size_t i = nodes.size() - nbCentralNodes; i < nodes.size(); ++i )
+ if ( nodes[i]->NbInverseElements() == 0 )
+ meshDS->RemoveFreeNode( nodes[i], smDS, /*fromGroups=*/true );
+
+ } // loop on theElements
+
+ if ( !theForce3d )
+ { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
+ // helper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
+ // helper.FixQuadraticElements( myError );
+ SMESH_MesherHelper( *myMesh ).FixQuadraticElements(myError);
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Convert quadratic elements to linear ones and remove quadratic nodes
+ * \return int - nb of checked elements
+ */
+//=======================================================================
+
+int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh * theSm,
+ SMDS_ElemIteratorPtr theItr,
+ const int theShapeID)
+{
+ int nbElem = 0;
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+ ElemFeatures elemType;
+ vector<const SMDS_MeshNode *> nodes;
+
+ while( theItr->more() )
+ {
+ const SMDS_MeshElement* elem = theItr->next();
+ nbElem++;
+ if( elem && elem->IsQuadratic())
+ {
+ // get elem data
+ int nbCornerNodes = elem->NbCornerNodes();
+ nodes.assign( elem->begin_nodes(), elem->end_nodes() );
+
+ elemType.Init( elem, /*basicOnly=*/false ).SetID( elem->GetID() ).SetQuad( false );
+
+ //remove a quadratic element
+ if ( !theSm || !theSm->Contains( elem ))
+ theSm = meshDS->MeshElements( elem->getshapeId() );
+ meshDS->RemoveFreeElement( elem, theSm, /*fromGroups=*/false );
+
+ // remove medium nodes
+ for ( size_t i = nbCornerNodes; i < nodes.size(); ++i )
+ if ( nodes[i]->NbInverseElements() == 0 )
+ meshDS->RemoveFreeNode( nodes[i], theSm );
+
+ // add a linear element
+ nodes.resize( nbCornerNodes );
+ SMDS_MeshElement * newElem = AddElement( nodes, elemType );
+ ReplaceElemInGroups(elem, newElem, meshDS);
+ if( theSm && newElem )
+ theSm->AddElement( newElem );
+ }
+ }
+ return nbElem;
+}
+
+//=======================================================================
+//function : ConvertFromQuadratic
+//purpose :
+//=======================================================================
+
+bool SMESH_MeshEditor::ConvertFromQuadratic()
+{
+ int nbCheckedElems = 0;
+ if ( myMesh->HasShapeToMesh() )
+ {
+ if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
+ {
+ SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
+ while ( smIt->more() ) {
+ SMESH_subMesh* sm = smIt->next();
+ if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() )
+ nbCheckedElems += removeQuadElem( smDS, smDS->GetElements(), sm->GetId() );
+ }
+ }
+ }
+
+ int totalNbElems =
+ GetMeshDS()->NbEdges() + GetMeshDS()->NbFaces() + GetMeshDS()->NbVolumes();
+ if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
+ {
+ SMESHDS_SubMesh *aSM = 0;
+ removeQuadElem( aSM, GetMeshDS()->elementsIterator(), 0 );
+ }
+
+ return true;
+}
+
+namespace
+{
+ //================================================================================
+ /*!
+ * \brief Return true if all medium nodes of the element are in the node set
+ */
+ //================================================================================
+
+ bool allMediumNodesIn(const SMDS_MeshElement* elem, TIDSortedNodeSet& nodeSet )
+ {
+ for ( int i = elem->NbCornerNodes(); i < elem->NbNodes(); ++i )
+ if ( !nodeSet.count( elem->GetNode(i) ))
+ return false;
+ return true;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Makes given elements linear
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
+{
+ if ( theElements.empty() ) return;
+
+ // collect IDs of medium nodes of theElements; some of these nodes will be removed
+ set<int> mediumNodeIDs;
+ TIDSortedElemSet::iterator eIt = theElements.begin();
+ for ( ; eIt != theElements.end(); ++eIt )
+ {
+ const SMDS_MeshElement* e = *eIt;
+ for ( int i = e->NbCornerNodes(); i < e->NbNodes(); ++i )
+ mediumNodeIDs.insert( e->GetNode(i)->GetID() );
+ }
+
+ // replace given elements by linear ones
+ SMDS_ElemIteratorPtr elemIt = SMESHUtils::elemSetIterator( theElements );
+ removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
+
+ // we need to convert remaining elements whose all medium nodes are in mediumNodeIDs
+ // except those elements sharing medium nodes of quadratic element whose medium nodes
+ // are not all in mediumNodeIDs
+
+ // get remaining medium nodes
+ TIDSortedNodeSet mediumNodes;
+ set<int>::iterator nIdsIt = mediumNodeIDs.begin();
+ for ( ; nIdsIt != mediumNodeIDs.end(); ++nIdsIt )
+ if ( const SMDS_MeshNode* n = GetMeshDS()->FindNode( *nIdsIt ))
+ mediumNodes.insert( mediumNodes.end(), n );
+
+ // find more quadratic elements to convert
+ TIDSortedElemSet moreElemsToConvert;
+ TIDSortedNodeSet::iterator nIt = mediumNodes.begin();
+ for ( ; nIt != mediumNodes.end(); ++nIt )
+ {
+ SMDS_ElemIteratorPtr invIt = (*nIt)->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* e = invIt->next();
+ if ( e->IsQuadratic() && allMediumNodesIn( e, mediumNodes ))
+ {
+ // find a more complex element including e and
+ // whose medium nodes are not in mediumNodes
+ bool complexFound = false;
+ for ( int type = e->GetType() + 1; type < SMDSAbs_0DElement; ++type )
+ {
+ SMDS_ElemIteratorPtr invIt2 =
+ (*nIt)->GetInverseElementIterator( SMDSAbs_ElementType( type ));
+ while ( invIt2->more() )
+ {
+ const SMDS_MeshElement* eComplex = invIt2->next();
+ if ( eComplex->IsQuadratic() && !allMediumNodesIn( eComplex, mediumNodes))
+ {
+ int nbCommonNodes = SMESH_MeshAlgos::NbCommonNodes( e, eComplex );
+ if ( nbCommonNodes == e->NbNodes())
+ {
+ complexFound = true;
+ type = SMDSAbs_NbElementTypes; // to quit from the outer loop
+ break;
+ }
+ }
+ }
+ }
+ if ( !complexFound )
+ moreElemsToConvert.insert( e );
+ }
+ }
+ }
+ elemIt = SMESHUtils::elemSetIterator( moreElemsToConvert );
+ removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
+}
+
+//=======================================================================
+//function : SewSideElements
+//purpose :
+//=======================================================================
+
+SMESH_MeshEditor::Sew_Error
+SMESH_MeshEditor::SewSideElements (TIDSortedElemSet& theSide1,
+ TIDSortedElemSet& theSide2,
+ const SMDS_MeshNode* theFirstNode1,
+ const SMDS_MeshNode* theFirstNode2,
+ const SMDS_MeshNode* theSecondNode1,
+ const SMDS_MeshNode* theSecondNode2)
+{
+ ClearLastCreated();
+
+ if ( theSide1.size() != theSide2.size() )
+ return SEW_DIFF_NB_OF_ELEMENTS;
+
+ Sew_Error aResult = SEW_OK;
+ // Algo:
+ // 1. Build set of faces representing each side
+ // 2. Find which nodes of the side 1 to merge with ones on the side 2
+ // 3. Replace nodes in elements of the side 1 and remove replaced nodes
+
+ // =======================================================================
+ // 1. Build set of faces representing each side:
+ // =======================================================================
+ // a. build set of nodes belonging to faces
+ // b. complete set of faces: find missing faces whose nodes are in set of nodes
+ // c. create temporary faces representing side of volumes if correspondent
+ // face does not exist
+
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+ // TODO algorithm not OK with vtkUnstructuredGrid: 2 meshes can't share nodes
+ //SMDS_Mesh aTmpFacesMesh; // try to use the same mesh
+ TIDSortedElemSet faceSet1, faceSet2;
+ set<const SMDS_MeshElement*> volSet1, volSet2;
+ set<const SMDS_MeshNode*> nodeSet1, nodeSet2;
+ TIDSortedElemSet * faceSetPtr[] = { &faceSet1, &faceSet2 };
+ set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
+ set<const SMDS_MeshNode*> * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
+ 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 ];
+ TIDSortedElemSet * faceSet = faceSetPtr[ iSide ];
+ set<const SMDS_MeshElement*> * volSet = volSetPtr [ iSide ];
+ set<const SMDS_MeshElement*>::iterator vIt;
+ TIDSortedElemSet::iterator eIt;
+ set<const SMDS_MeshNode*>::iterator nIt;
+
+ // check that given nodes belong to given elements
+ const SMDS_MeshNode* n1 = ( iSide == 0 ) ? theFirstNode1 : theFirstNode2;
+ const SMDS_MeshNode* n2 = ( iSide == 0 ) ? theSecondNode1 : theSecondNode2;
+ int firstIndex = -1, secondIndex = -1;
+ for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
+ const SMDS_MeshElement* elem = *eIt;
+ if ( firstIndex < 0 ) firstIndex = elem->GetNodeIndex( n1 );
+ if ( secondIndex < 0 ) secondIndex = elem->GetNodeIndex( n2 );
+ if ( firstIndex > -1 && secondIndex > -1 ) break;
+ }
+ if ( firstIndex < 0 || secondIndex < 0 ) {
+ // we can simply return until temporary faces created
+ return (iSide == 0 ) ? SEW_BAD_SIDE1_NODES : SEW_BAD_SIDE2_NODES;
+ }
+
+ // -----------------------------------------------------------
+ // 1a. Collect nodes of existing faces
+ // and build set of face nodes in order to detect missing
+ // faces corresponding to sides of volumes
+ // -----------------------------------------------------------
+
+ set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
+
+ // loop on the given element of a side
+ for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
+ //const SMDS_MeshElement* elem = *eIt;
+ const SMDS_MeshElement* elem = *eIt;
+ if ( elem->GetType() == SMDSAbs_Face ) {
+ faceSet->insert( elem );
+ set <const SMDS_MeshNode*> faceNodeSet;
+ SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+ while ( nodeIt->more() ) {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ nodeSet->insert( n );
+ faceNodeSet.insert( n );
+ }
+ setOfFaceNodeSet.insert( faceNodeSet );
+ }
+ else if ( elem->GetType() == SMDSAbs_Volume )
+ volSet->insert( elem );
+ }
+ // ------------------------------------------------------------------------------
+ // 1b. Complete set of faces: find missing faces whose nodes are in set of nodes
+ // ------------------------------------------------------------------------------
+
+ for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
+ SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() ) { // loop on faces sharing a node
+ const SMDS_MeshElement* f = fIt->next();
+ if ( faceSet->find( f ) == faceSet->end() ) {
+ // check if all nodes are in nodeSet and
+ // complete setOfFaceNodeSet if they are
+ set <const SMDS_MeshNode*> faceNodeSet;
+ SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
+ bool allInSet = true;
+ while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ if ( nodeSet->find( n ) == nodeSet->end() )
+ allInSet = false;
+ else
+ faceNodeSet.insert( n );
+ }
+ if ( allInSet ) {
+ faceSet->insert( f );
+ setOfFaceNodeSet.insert( faceNodeSet );
+ }
+ }
+ }
+ }
+
+ // -------------------------------------------------------------------------
+ // 1c. Create temporary faces representing sides of volumes if correspondent
+ // face does not exist
+ // -------------------------------------------------------------------------
+
+ if ( !volSet->empty() ) {
+ //int nodeSetSize = nodeSet->size();
+
+ // loop on given volumes
+ for ( vIt = volSet->begin(); vIt != volSet->end(); vIt++ ) {
+ SMDS_VolumeTool vol (*vIt);
+ // loop on volume faces: find free faces
+ // --------------------------------------
+ list<const SMDS_MeshElement* > freeFaceList;
+ for ( iFace = 0; iFace < vol.NbFaces(); iFace++ ) {
+ if ( !vol.IsFreeFace( iFace ))
+ continue;
+ // check if there is already a face with same nodes in a face set
+ const SMDS_MeshElement* aFreeFace = 0;
+ const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iFace );
+ int nbNodes = vol.NbFaceNodes( iFace );
+ set <const SMDS_MeshNode*> faceNodeSet;
+ vol.GetFaceNodes( iFace, faceNodeSet );
+ bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
+ if ( isNewFace ) {
+ // no such a face is given but it still can exist, check it
+ vector<const SMDS_MeshNode *> nodes ( fNodes, fNodes + nbNodes);
+ aFreeFace = aMesh->FindElement( nodes, SMDSAbs_Face, /*noMedium=*/false );
+ }
+ if ( !aFreeFace ) {
+ // create a temporary face
+ if ( nbNodes == 3 ) {
+ //aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
+ aFreeFace = aMesh->AddFace( fNodes[0],fNodes[1],fNodes[2] );
+ }
+ else if ( nbNodes == 4 ) {
+ //aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
+ aFreeFace = aMesh->AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
+ }
+ else {
+ vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
+ //aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
+ aFreeFace = aMesh->AddPolygonalFace(poly_nodes);
+ }
+ if ( aFreeFace )
+ tempFaceList.push_back( aFreeFace );
+ }
+
+ if ( aFreeFace )
+ freeFaceList.push_back( aFreeFace );
+
+ } // loop on faces of a volume
+
+ // 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;
+ list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
+ while ( fIt != freeFaceList.end() ) { // loop on free faces
+ int nbSharedNodes = 0;
+ SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
+ while ( nodeIt->more() ) { // loop on free face nodes
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
+ while ( invElemIt->more() ) {
+ const SMDS_MeshElement* e = invElemIt->next();
+ nbSharedNodes += faceSet->count( e );
+ nbSharedNodes += elemSet->count( e );
+ }
+ }
+ if ( nbSharedNodes > maxNbNodes ) {
+ maxNbNodes = nbSharedNodes;
+ freeFaceList.erase( freeFaceList.begin(), fIt++ );
+ }
+ else if ( nbSharedNodes == maxNbNodes ) {
+ fIt++;
+ }
+ else {
+ freeFaceList.erase( fIt++ ); // here fIt++ occurs before erase
+ }
+ }
+ if ( freeFaceList.size() > 1 )
+ {
+ // could not choose one face, use another way
+ // choose a face most close to the bary center of the opposite side
+ gp_XYZ aBC( 0., 0., 0. );
+ set <const SMDS_MeshNode*> addedNodes;
+ TIDSortedElemSet * elemSet2 = elemSetPtr[ 1 - iSide ];
+ eIt = elemSet2->begin();
+ for ( eIt = elemSet2->begin(); eIt != elemSet2->end(); eIt++ ) {
+ SMDS_ElemIteratorPtr nodeIt = (*eIt)->nodesIterator();
+ while ( nodeIt->more() ) { // loop on free face nodes
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ if ( addedNodes.insert( n ).second )
+ aBC += gp_XYZ( n->X(),n->Y(),n->Z() );
+ }
+ }
+ aBC /= addedNodes.size();
+ double minDist = DBL_MAX;
+ fIt = freeFaceList.begin();
+ while ( fIt != freeFaceList.end() ) { // loop on free faces
+ double dist = 0;
+ SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
+ while ( nodeIt->more() ) { // loop on free face nodes
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ gp_XYZ p( n->X(),n->Y(),n->Z() );
+ dist += ( aBC - p ).SquareModulus();
+ }
+ if ( dist < minDist ) {
+ minDist = dist;
+ freeFaceList.erase( freeFaceList.begin(), fIt++ );
+ }
+ else
+ fIt = freeFaceList.erase( fIt++ );
+ }
+ }
+ } // choose one of several free faces of a volume
+
+ if ( freeFaceList.size() == 1 ) {
+ const SMDS_MeshElement* aFreeFace = freeFaceList.front();
+ faceSet->insert( aFreeFace );
+ // complete a node set with nodes of a found free face
+ // for ( iNode = 0; iNode < ; iNode++ )
+ // nodeSet->insert( fNodes[ iNode ] );
+ }
+
+ } // loop on volumes of a side
+
+ // // complete a set of faces if new nodes in a nodeSet appeared
+ // // ----------------------------------------------------------
+ // if ( nodeSetSize != nodeSet->size() ) {
+ // for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
+ // SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
+ // while ( fIt->more() ) { // loop on faces sharing a node
+ // const SMDS_MeshElement* f = fIt->next();
+ // if ( faceSet->find( f ) == faceSet->end() ) {
+ // // check if all nodes are in nodeSet and
+ // // complete setOfFaceNodeSet if they are
+ // set <const SMDS_MeshNode*> faceNodeSet;
+ // SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
+ // bool allInSet = true;
+ // while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
+ // const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ // if ( nodeSet->find( n ) == nodeSet->end() )
+ // allInSet = false;
+ // else
+ // faceNodeSet.insert( n );
+ // }
+ // if ( allInSet ) {
+ // faceSet->insert( f );
+ // setOfFaceNodeSet.insert( faceNodeSet );
+ // }
+ // }
+ // }
+ // }
+ // }
+ } // Create temporary faces, if there are volumes given
+ } // loop on sides
+
+ if ( faceSet1.size() != faceSet2.size() ) {
+ // delete temporary faces: they are in reverseElements of actual nodes
+ // 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);
+ MESSAGE("Diff nb of faces");
+ return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ }
+
+ // ============================================================
+ // 2. Find nodes to merge:
+ // bind a node to remove to a node to put instead
+ // ============================================================
+
+ TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
+ if ( theFirstNode1 != theFirstNode2 )
+ nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
+ if ( theSecondNode1 != theSecondNode2 )
+ nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
+
+ LinkID_Gen aLinkID_Gen( GetMeshDS() );
+ set< long > linkIdSet; // links to process
+ linkIdSet.insert( aLinkID_Gen.GetLinkID( theFirstNode1, theSecondNode1 ));
+
+ typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > NLink;
+ list< NLink > linkList[2];
+ linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
+ linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
+ // 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]++ )
+ {
+ NLink link[] = { *linkIt[0], *linkIt[1] };
+ long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
+ if ( !linkIdSet.count( linkID ) )
+ continue;
+
+ // by links, find faces in the face sets,
+ // and find indices of link nodes in the found faces;
+ // in a face set, there is only one or no face sharing a link
+ // ---------------------------------------------------------------
+
+ const SMDS_MeshElement* face[] = { 0, 0 };
+ 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;
+ //cout << "Side " << iSide << " ";
+ //cout << "L( " << n1->GetID() << ", " << n2->GetID() << " ) " << endl;
+ // find a face by two link nodes
+ face[ iSide ] = SMESH_MeshAlgos::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
+ SMDS_MeshElement::iterator nIt( face[ iSide ]->interlacedNodesIterator() ), nEnd;
+ fnodes[ iSide ].assign( nIt, nEnd );
+ fnodes[ iSide ].push_back( fnodes[ iSide ].front());
+ }
+ }
+
+ // check similarity of elements of the sides
+ if (aResult == SEW_OK && (( face[0] && !face[1] ) || ( !face[0] && face[1] ))) {
+ MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
+ if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
+ aResult = ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
+ }
+ else {
+ aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ }
+ break; // do not return because it's necessary to remove tmp faces
+ }
+
+ // set nodes to merge
+ // -------------------
+
+ if ( face[0] && face[1] ) {
+ 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 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;
+ }
+ 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
+ // -----------------------------------------
+
+ for ( iNode = 0; iNode < nbNodes; iNode++ ) {
+ 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 = 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() ||
+ !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
+ MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
+ " " << (faceSetPtr[1]->empty()));
+ aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ }
+
+ // ====================================================================
+ // 3. Replace nodes in elements of the side 1 and remove replaced 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);
+
+ if ( aResult != SEW_OK)
+ return aResult;
+
+ list< int > nodeIDsToRemove;
+ vector< const SMDS_MeshNode*> nodes;
+ ElemFeatures elemType;
+
+ // loop on nodes replacement map
+ TNodeNodeMap::iterator nReplaceMapIt = nReplaceMap.begin(), nnIt;
+ for ( ; nReplaceMapIt != nReplaceMap.end(); nReplaceMapIt++ )
+ if ( (*nReplaceMapIt).first != (*nReplaceMapIt).second )
+ {
+ const SMDS_MeshNode* nToRemove = (*nReplaceMapIt).first;
+ nodeIDsToRemove.push_back( nToRemove->GetID() );
+ // loop on elements sharing nToRemove
+ SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
+ while ( invElemIt->more() ) {
+ const SMDS_MeshElement* e = invElemIt->next();
+ // get a new suite of nodes: make replacement
+ int nbReplaced = 0, i = 0, nbNodes = e->NbNodes();
+ nodes.resize( nbNodes );
+ SMDS_ElemIteratorPtr nIt = e->nodesIterator();
+ while ( nIt->more() ) {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
+ nnIt = nReplaceMap.find( n );
+ if ( nnIt != nReplaceMap.end() ) {
+ nbReplaced++;
+ n = (*nnIt).second;
+ }
+ nodes[ i++ ] = n;
+ }
+ // if ( nbReplaced == nbNodes && e->GetType() == SMDSAbs_Face )
+ // elemIDsToRemove.push_back( e->GetID() );
+ // else
+ if ( nbReplaced )
+ {
+ elemType.Init( e, /*basicOnly=*/false ).SetID( e->GetID() );
+ aMesh->RemoveElement( e );
+
+ if ( SMDS_MeshElement* newElem = this->AddElement( nodes, elemType ))
+ {
+ AddToSameGroups( newElem, e, aMesh );
+ if ( int aShapeId = e->getshapeId() )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ }
+ }
+ }
+ }
+
+ Remove( nodeIDsToRemove, true );
+
+ return aResult;
+}
+
+//================================================================================
+/*!
+ * \brief Find corresponding nodes in two sets of faces
+ * \param theSide1 - first face set
+ * \param theSide2 - second first face
+ * \param theFirstNode1 - a boundary node of set 1
+ * \param theFirstNode2 - a node of set 2 corresponding to theFirstNode1
+ * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
+ * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
+ * \param nReplaceMap - output map of corresponding nodes
+ * \return bool - is a success or not
+ */
+//================================================================================
+
+#ifdef _DEBUG_
+//#define DEBUG_MATCHING_NODES
+#endif
+
+SMESH_MeshEditor::Sew_Error
+SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
+ set<const SMDS_MeshElement*>& theSide2,
+ const SMDS_MeshNode* theFirstNode1,
+ const SMDS_MeshNode* theFirstNode2,
+ const SMDS_MeshNode* theSecondNode1,
+ const SMDS_MeshNode* theSecondNode2,
+ TNodeNodeMap & nReplaceMap)
+{
+ set<const SMDS_MeshElement*> * faceSetPtr[] = { &theSide1, &theSide2 };
+
+ nReplaceMap.clear();
+ //if ( theFirstNode1 != theFirstNode2 )
+ nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
+ //if ( theSecondNode1 != theSecondNode2 )
+ nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
+
+ set< SMESH_TLink > linkSet; // set of nodes where order of nodes is ignored
+ linkSet.insert( SMESH_TLink( theFirstNode1, theSecondNode1 ));
+
+ list< NLink > linkList[2];
+ linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
+ linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
+
+ // 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]++ ) {
+ NLink link[] = { *linkIt[0], *linkIt[1] };
+ if ( linkSet.find( link[0] ) == linkSet.end() )
+ continue;
+
+ // by links, find faces in the face sets,
+ // and find indices of link nodes in the found faces;
+ // in a face set, there is only one or no face sharing a link
+ // ---------------------------------------------------------------
+
+ const SMDS_MeshElement* face[] = { 0, 0 };
+ list<const SMDS_MeshNode*> notLinkNodes[2];
+ //bool reverse[] = { false, false }; // order of notLinkNodes
+ int nbNodes[2];
+ for ( int 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* > facesOfNode1;
+ for ( int iNode = 0; iNode < 2; iNode++ ) // loop on 2 nodes of a link
+ {
+ // during a loop of the first node, we find all faces around n1,
+ // during a loop of the second node, we find one face sharing both n1 and n2
+ const SMDS_MeshNode* n = iNode ? 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
+ ! facesOfNode1.insert( f ).second ) // f encounters twice
+ {
+ if ( face[ iSide ] ) {
+ MESSAGE( "2 faces per link " );
+ return ( iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
+ }
+ face[ iSide ] = f;
+ faceSet->erase( f );
+
+ // get not link nodes
+ int nbN = f->NbNodes();
+ if ( f->IsQuadratic() )
+ nbN /= 2;
+ nbNodes[ iSide ] = nbN;
+ list< const SMDS_MeshNode* > & nodes = notLinkNodes[ iSide ];
+ int i1 = f->GetNodeIndex( n1 );
+ int i2 = f->GetNodeIndex( n2 );
+ int iEnd = nbN, iBeg = -1, iDelta = 1;
+ bool reverse = ( Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1 );
+ if ( reverse ) {
+ std::swap( iEnd, iBeg ); iDelta = -1;
+ }
+ int i = i2;
+ while ( true ) {
+ i += iDelta;
+ if ( i == iEnd ) i = iBeg + iDelta;
+ if ( i == i1 ) break;
+ nodes.push_back ( f->GetNode( i ) );
+ }
+ }
+ }
+ }
+ }
+ // check similarity of elements of the sides
+ if (( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
+ MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
+ if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
+ return ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
+ }
+ else {
+ return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ }
+ }
+
+ // set nodes to merge
+ // -------------------
+
+ if ( face[0] && face[1] ) {
+ if ( nbNodes[0] != nbNodes[1] ) {
+ MESSAGE("Diff nb of face nodes");
+ return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ }
+#ifdef DEBUG_MATCHING_NODES
+ MESSAGE ( " Link 1: " << link[0].first->GetID() <<" "<< link[0].second->GetID()
+ << " F 1: " << face[0] << "| Link 2: " << link[1].first->GetID() <<" "
+ << link[1].second->GetID() << " F 2: " << face[1] << " | Bind: " ) ;
+#endif
+ int nbN = nbNodes[0];
+ {
+ list<const SMDS_MeshNode*>::iterator n1 = notLinkNodes[0].begin();
+ list<const SMDS_MeshNode*>::iterator n2 = notLinkNodes[1].begin();
+ for ( int i = 0 ; i < nbN - 2; ++i ) {
+#ifdef DEBUG_MATCHING_NODES
+ MESSAGE ( (*n1)->GetID() << " to " << (*n2)->GetID() );
+#endif
+ nReplaceMap.insert( make_pair( *(n1++), *(n2++) ));
+ }
+ }
+
+ // add other links of the face 1 to linkList
+ // -----------------------------------------
+
+ const SMDS_MeshElement* f0 = face[0];
+ const SMDS_MeshNode* n1 = f0->GetNode( nbN - 1 );
+ for ( int i = 0; i < nbN; i++ )
+ {
+ const SMDS_MeshNode* n2 = f0->GetNode( i );
+ pair< set< SMESH_TLink >::iterator, bool > iter_isnew =
+ linkSet.insert( SMESH_TLink( n1, n2 ));
+ if ( !iter_isnew.second ) { // already in a set: no need to process
+ linkSet.erase( iter_isnew.first );
+ }
+ else // new in set == encountered for the first time: add
+ {
+#ifdef DEBUG_MATCHING_NODES
+ MESSAGE ( "Add link 1: " << n1->GetID() << " " << n2->GetID() << " "
+ << " | link 2: " << nReplaceMap[n1]->GetID() << " " << nReplaceMap[n2]->GetID() << " " );
+#endif
+ linkList[0].push_back ( NLink( n1, n2 ));
+ linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
+ }
+ n1 = n2;
+ }
+ } // 2 faces found
+ } // loop on link lists
+
+ return SEW_OK;
+}
+
+namespace // automatically find theAffectedElems for DoubleNodes()
+{
+ bool isOut( const SMDS_MeshNode* n, const gp_XYZ& norm, const SMDS_MeshElement* elem );
+
+ //--------------------------------------------------------------------------------
+ // Nodes shared by adjacent FissureBorder's.
+ // 1 node if FissureBorder separates faces
+ // 2 nodes if FissureBorder separates volumes
+ struct SubBorder
+ {
+ const SMDS_MeshNode* _nodes[2];
+ int _nbNodes;
+
+ SubBorder( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2 = 0 )
+ {
+ _nodes[0] = n1;
+ _nodes[1] = n2;
+ _nbNodes = bool( n1 ) + bool( n2 );
+ if ( _nbNodes == 2 && n1 > n2 )
+ std::swap( _nodes[0], _nodes[1] );
+ }
+ bool operator<( const SubBorder& other ) const
+ {
+ for ( int i = 0; i < _nbNodes; ++i )
+ {
+ if ( _nodes[i] < other._nodes[i] ) return true;
+ if ( _nodes[i] > other._nodes[i] ) return false;
+ }
+ return false;
+ }
+ };
+
+ //--------------------------------------------------------------------------------
+ // Map a SubBorder to all FissureBorder it bounds
+ struct FissureBorder;
+ typedef std::map< SubBorder, std::vector< FissureBorder* > > TBorderLinks;
+ typedef TBorderLinks::iterator TMappedSub;
+
+ //--------------------------------------------------------------------------------
+ /*!
+ * \brief Element border (volume facet or face edge) at a fissure
+ */
+ struct FissureBorder
+ {
+ std::vector< const SMDS_MeshNode* > _nodes; // border nodes
+ const SMDS_MeshElement* _elems[2]; // volume or face adjacent to fissure
+
+ std::vector< TMappedSub > _mappedSubs; // Sub() in TBorderLinks map
+ std::vector< const SMDS_MeshNode* > _sortedNodes; // to compare FissureBorder's
+
+ FissureBorder( FissureBorder && from ) // move constructor
+ {
+ std::swap( _nodes, from._nodes );
+ std::swap( _sortedNodes, from._sortedNodes );
+ _elems[0] = from._elems[0];
+ _elems[1] = from._elems[1];
+ }
+
+ FissureBorder( const SMDS_MeshElement* elemToDuplicate,
+ std::vector< const SMDS_MeshElement* > & adjElems)
+ : _nodes( elemToDuplicate->NbCornerNodes() )
+ {
+ for ( size_t i = 0; i < _nodes.size(); ++i )
+ _nodes[i] = elemToDuplicate->GetNode( i );
+
+ SMDSAbs_ElementType type = SMDSAbs_ElementType( elemToDuplicate->GetType() + 1 );
+ findAdjacent( type, adjElems );
+ }
+
+ FissureBorder( const SMDS_MeshNode** nodes,
+ const size_t nbNodes,
+ const SMDSAbs_ElementType adjElemsType,
+ std::vector< const SMDS_MeshElement* > & adjElems)
+ : _nodes( nodes, nodes + nbNodes )
+ {
+ findAdjacent( adjElemsType, adjElems );
+ }
+
+ void findAdjacent( const SMDSAbs_ElementType adjElemsType,
+ std::vector< const SMDS_MeshElement* > & adjElems)
+ {
+ _elems[0] = _elems[1] = 0;
+ adjElems.clear();
+ if ( SMDS_Mesh::GetElementsByNodes( _nodes, adjElems, adjElemsType ))
+ for ( size_t i = 0; i < adjElems.size() && i < 2; ++i )
+ _elems[i] = adjElems[i];
+ }
+
+ bool operator<( const FissureBorder& other ) const
+ {
+ return GetSortedNodes() < other.GetSortedNodes();
+ }
+
+ const std::vector< const SMDS_MeshNode* >& GetSortedNodes() const
+ {
+ if ( _sortedNodes.empty() && !_nodes.empty() )
+ {
+ FissureBorder* me = const_cast<FissureBorder*>( this );
+ me->_sortedNodes = me->_nodes;
+ std::sort( me->_sortedNodes.begin(), me->_sortedNodes.end() );
+ }
+ return _sortedNodes;
+ }
+
+ size_t NbSub() const
+ {
+ return _nodes.size();
+ }
+
+ SubBorder Sub(size_t i) const
+ {
+ return SubBorder( _nodes[i], NbSub() > 2 ? _nodes[ (i+1)%NbSub() ] : 0 );
+ }
+
+ void AddSelfTo( TBorderLinks& borderLinks )
+ {
+ _mappedSubs.resize( NbSub() );
+ for ( size_t i = 0; i < NbSub(); ++i )
+ {
+ TBorderLinks::iterator s2b =
+ borderLinks.insert( std::make_pair( Sub(i), TBorderLinks::mapped_type() )).first;
+ s2b->second.push_back( this );
+ _mappedSubs[ i ] = s2b;
+ }
+ }
+
+ void Clear()
+ {
+ _nodes.clear();
+ }
+
+ const SMDS_MeshElement* GetMarkedElem() const
+ {
+ if ( _nodes.empty() ) return 0; // cleared
+ if ( _elems[0] && _elems[0]->isMarked() ) return _elems[0];
+ if ( _elems[1] && _elems[1]->isMarked() ) return _elems[1];
+ return 0;
+ }
+
+ gp_XYZ GetNorm() const // normal to the border
+ {
+ gp_XYZ norm;
+ if ( _nodes.size() == 2 )
+ {
+ gp_XYZ avgNorm( 0,0,0 ); // sum of normals of adjacent faces
+ if ( SMESH_MeshAlgos::FaceNormal( _elems[0], norm ))
+ avgNorm += norm;
+ if ( SMESH_MeshAlgos::FaceNormal( _elems[1], norm ))
+ avgNorm += norm;
+
+ gp_XYZ bordDir( SMESH_NodeXYZ( _nodes[0] ) - SMESH_NodeXYZ( _nodes[1] ));
+ norm = bordDir ^ avgNorm;
+ }
+ else
+ {
+ SMESH_NodeXYZ p0( _nodes[0] );
+ SMESH_NodeXYZ p1( _nodes[1] );
+ SMESH_NodeXYZ p2( _nodes[2] );
+ norm = ( p0 - p1 ) ^ ( p2 - p1 );
+ }
+ if ( isOut( _nodes[0], norm, GetMarkedElem() ))
+ norm.Reverse();
+
+ return norm;
+ }
+
+ void ChooseSide() // mark an _elem located at positive side of fissure
+ {
+ _elems[0]->setIsMarked( true );
+ gp_XYZ norm = GetNorm();
+ double maxX = norm.Coord(1);
+ if ( Abs( maxX ) < Abs( norm.Coord(2)) ) maxX = norm.Coord(2);
+ if ( Abs( maxX ) < Abs( norm.Coord(3)) ) maxX = norm.Coord(3);
+ if ( maxX < 0 )
+ {
+ _elems[0]->setIsMarked( false );
+ _elems[1]->setIsMarked( true );
+ }
+ }
+
+ }; // struct FissureBorder