+ }
+ }
+ list< list<SMESH_MeshEditor_PathPoint> >::iterator itLLPP = LLPPs.begin();
+ list<SMESH_MeshEditor_PathPoint> firstList = *itLLPP;
+ list<SMESH_MeshEditor_PathPoint>::iterator itPP = firstList.begin();
+ for(; itPP!=firstList.end(); itPP++) {
+ fullList.push_back( *itPP );
+ }
+ SMESH_MeshEditor_PathPoint PP1 = fullList.back();
+ fullList.pop_back();
+ itLLPP++;
+ for(; itLLPP!=LLPPs.end(); itLLPP++) {
+ list<SMESH_MeshEditor_PathPoint> currList = *itLLPP;
+ itPP = currList.begin();
+ SMESH_MeshEditor_PathPoint PP2 = currList.front();
+ gp_Dir D1 = PP1.Tangent();
+ gp_Dir D2 = PP2.Tangent();
+ gp_Dir Dnew( gp_Vec( (D1.X()+D2.X())/2, (D1.Y()+D2.Y())/2,
+ (D1.Z()+D2.Z())/2 ) );
+ PP1.SetTangent(Dnew);
+ fullList.push_back(PP1);
+ itPP++;
+ for(; itPP!=firstList.end(); itPP++) {
+ fullList.push_back( *itPP );
+ }
+ PP1 = fullList.back();
+ fullList.pop_back();
+ }
+ // if wire not closed
+ fullList.push_back(PP1);
+ // else ???
+ }
+ else {
+ return EXTR_BAD_PATH_SHAPE;
+ }
+
+ return MakeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
+ theHasRefPoint, theRefPoint, theMakeGroups);
+}
+
+
+//=======================================================================
+//function : ExtrusionAlongTrack
+//purpose :
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet theElements[2],
+ SMESH_Mesh* theTrack,
+ const SMDS_MeshNode* theN1,
+ const bool theHasAngles,
+ list<double>& theAngles,
+ const bool theLinearVariation,
+ const bool theHasRefPoint,
+ const gp_Pnt& theRefPoint,
+ const bool theMakeGroups)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ int aNbE;
+ std::list<double> aPrms;
+ TIDSortedElemSet::iterator itElem;
+
+ gp_XYZ aGC;
+ TopoDS_Edge aTrackEdge;
+ TopoDS_Vertex aV1, aV2;
+
+ SMDS_ElemIteratorPtr aItE;
+ SMDS_NodeIteratorPtr aItN;
+ SMDSAbs_ElementType aTypeE;
+
+ TNodeOfNodeListMap mapNewNodes;
+
+ // 1. Check data
+ aNbE = theElements[0].size() + theElements[1].size();
+ // nothing to do
+ if ( !aNbE )
+ return EXTR_NO_ELEMENTS;
+
+ // 1.1 Track Pattern
+ ASSERT( theTrack );
+
+ SMESHDS_Mesh* pMeshDS = theTrack->GetMeshDS();
+
+ aItE = pMeshDS->elementsIterator();
+ while ( aItE->more() ) {
+ const SMDS_MeshElement* pE = aItE->next();
+ aTypeE = pE->GetType();
+ // Pattern must contain links only
+ if ( aTypeE != SMDSAbs_Edge )
+ return EXTR_PATH_NOT_EDGE;
+ }
+
+ list<SMESH_MeshEditor_PathPoint> fullList;
+
+ const TopoDS_Shape& aS = theTrack->GetShapeToMesh();
+
+ if ( !theTrack->HasShapeToMesh() ) {
+ //Mesh without shape
+ const SMDS_MeshNode* currentNode = NULL;
+ const SMDS_MeshNode* prevNode = theN1;
+ std::vector<const SMDS_MeshNode*> aNodesList;
+ aNodesList.push_back(theN1);
+ int nbEdges = 0, conn=0;
+ const SMDS_MeshElement* prevElem = NULL;
+ const SMDS_MeshElement* currentElem = NULL;
+ int totalNbEdges = theTrack->NbEdges();
+ SMDS_ElemIteratorPtr nIt;
+
+ //check start node
+ if( !theTrack->GetMeshDS()->Contains(theN1) ) {
+ return EXTR_BAD_STARTING_NODE;
+ }
+
+ conn = nbEdgeConnectivity(theN1);
+ if( conn != 1 )
+ return EXTR_PATH_NOT_EDGE;
+
+ aItE = theN1->GetInverseElementIterator();
+ prevElem = aItE->next();
+ currentElem = prevElem;
+ //Get all nodes
+ if(totalNbEdges == 1 ) {
+ nIt = currentElem->nodesIterator();
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ if(currentNode == prevNode)
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ aNodesList.push_back(currentNode);
+ } else {
+ nIt = currentElem->nodesIterator();
+ while( nIt->more() ) {
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ if(currentNode == prevNode)
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ aNodesList.push_back(currentNode);
+
+ //case of the closed mesh
+ if(currentNode == theN1) {
+ nbEdges++;
+ break;
+ }
+
+ conn = nbEdgeConnectivity(currentNode);
+ if(conn > 2) {
+ return EXTR_PATH_NOT_EDGE;
+ }else if( conn == 1 && nbEdges > 0 ) {
+ //End of the path
+ nbEdges++;
+ break;
+ }else {
+ prevNode = currentNode;
+ aItE = currentNode->GetInverseElementIterator();
+ currentElem = aItE->next();
+ if( currentElem == prevElem)
+ currentElem = aItE->next();
+ nIt = currentElem->nodesIterator();
+ prevElem = currentElem;
+ nbEdges++;
+ }
+ }
+ }
+
+ if(nbEdges != totalNbEdges)
+ return EXTR_PATH_NOT_EDGE;
+
+ TopTools_SequenceOfShape Edges;
+ double x1,x2,y1,y2,z1,z2;
+ list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
+ int startNid = theN1->GetID();
+ for(int i = 1; i < aNodesList.size(); i++) {
+ x1 = aNodesList[i-1]->X();x2 = aNodesList[i]->X();
+ y1 = aNodesList[i-1]->Y();y2 = aNodesList[i]->Y();
+ z1 = aNodesList[i-1]->Z();z2 = aNodesList[i]->Z();
+ TopoDS_Edge e = BRepBuilderAPI_MakeEdge(gp_Pnt(x1,y1,z1),gp_Pnt(x2,y2,z2));
+ list<SMESH_MeshEditor_PathPoint> LPP;
+ aPrms.clear();
+ MakeEdgePathPoints(aPrms, e, (aNodesList[i-1]->GetID()==startNid), LPP);
+ LLPPs.push_back(LPP);
+ if( aNodesList[i-1]->GetID() == startNid ) startNid = aNodesList[i]->GetID();
+ else startNid = aNodesList[i-1]->GetID();
+
+ }
+
+ list< list<SMESH_MeshEditor_PathPoint> >::iterator itLLPP = LLPPs.begin();
+ list<SMESH_MeshEditor_PathPoint> firstList = *itLLPP;
+ list<SMESH_MeshEditor_PathPoint>::iterator itPP = firstList.begin();
+ for(; itPP!=firstList.end(); itPP++) {
+ fullList.push_back( *itPP );
+ }
+
+ SMESH_MeshEditor_PathPoint PP1 = fullList.back();
+ SMESH_MeshEditor_PathPoint PP2;
+ fullList.pop_back();
+ itLLPP++;
+ for(; itLLPP!=LLPPs.end(); itLLPP++) {
+ list<SMESH_MeshEditor_PathPoint> currList = *itLLPP;
+ itPP = currList.begin();
+ PP2 = currList.front();
+ gp_Dir D1 = PP1.Tangent();
+ gp_Dir D2 = PP2.Tangent();
+ gp_Dir Dnew( gp_Vec( (D1.X()+D2.X())/2, (D1.Y()+D2.Y())/2,
+ (D1.Z()+D2.Z())/2 ) );
+ PP1.SetTangent(Dnew);
+ fullList.push_back(PP1);
+ itPP++;
+ for(; itPP!=currList.end(); itPP++) {
+ fullList.push_back( *itPP );
+ }
+ PP1 = fullList.back();
+ fullList.pop_back();
+ }
+ fullList.push_back(PP1);
+
+ } // Sub-shape for the Pattern must be an Edge or Wire
+ else if( aS.ShapeType() == TopAbs_EDGE ) {
+ aTrackEdge = TopoDS::Edge( aS );
+ // the Edge must not be degenerated
+ if ( SMESH_Algo::isDegenerated( aTrackEdge ) )
+ return EXTR_BAD_PATH_SHAPE;
+ TopExp::Vertices( aTrackEdge, aV1, aV2 );
+ const SMDS_MeshNode* aN1 = SMESH_Algo::VertexNode( aV1, pMeshDS );
+ const SMDS_MeshNode* aN2 = SMESH_Algo::VertexNode( aV2, pMeshDS );
+ // starting node must be aN1 or aN2
+ if ( !( aN1 == theN1 || aN2 == theN1 ) )
+ return EXTR_BAD_STARTING_NODE;
+ aItN = pMeshDS->nodesIterator();
+ while ( aItN->more() ) {
+ const SMDS_MeshNode* pNode = aItN->next();
+ if( pNode==aN1 || pNode==aN2 ) continue;
+ const SMDS_EdgePosition* pEPos =
+ static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
+ double aT = pEPos->GetUParameter();
+ aPrms.push_back( aT );
+ }
+ //Extrusion_Error err =
+ MakeEdgePathPoints(aPrms, aTrackEdge, (aN1==theN1), fullList);
+ }
+ else if( aS.ShapeType() == TopAbs_WIRE ) {
+ list< SMESH_subMesh* > LSM;
+ TopTools_SequenceOfShape Edges;
+ TopExp_Explorer eExp(aS, TopAbs_EDGE);
+ for(; eExp.More(); eExp.Next()) {
+ TopoDS_Edge E = TopoDS::Edge( eExp.Current() );
+ if( SMESH_Algo::isDegenerated(E) ) continue;
+ SMESH_subMesh* SM = theTrack->GetSubMesh(E);
+ if(SM) {
+ LSM.push_back(SM);
+ Edges.Append(E);
+ }
+ }
+ list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
+ TopoDS_Vertex aVprev;
+ TColStd_MapOfInteger UsedNums;
+ int NbEdges = Edges.Length();
+ int i = 1;
+ for(; i<=NbEdges; i++) {
+ int k = 0;
+ list< SMESH_subMesh* >::iterator itLSM = LSM.begin();
+ for(; itLSM!=LSM.end(); itLSM++) {
+ k++;
+ if(UsedNums.Contains(k)) continue;
+ aTrackEdge = TopoDS::Edge( Edges.Value(k) );
+ SMESH_subMesh* locTrack = *itLSM;
+ SMESHDS_SubMesh* locMeshDS = locTrack->GetSubMeshDS();
+ TopExp::Vertices( aTrackEdge, aV1, aV2 );
+ bool aN1isOK = false, aN2isOK = false;
+ if ( aVprev.IsNull() ) {
+ // if previous vertex is not yet defined, it means that we in the beginning of wire
+ // and we have to find initial vertex corresponding to starting node theN1
+ const SMDS_MeshNode* aN1 = SMESH_Algo::VertexNode( aV1, pMeshDS );
+ const SMDS_MeshNode* aN2 = SMESH_Algo::VertexNode( aV2, pMeshDS );
+ // starting node must be aN1 or aN2
+ aN1isOK = ( aN1 && aN1 == theN1 );
+ aN2isOK = ( aN2 && aN2 == theN1 );
+ }
+ else {
+ // we have specified ending vertex of the previous edge on the previous iteration
+ // and we have just to check that it corresponds to any vertex in current segment
+ aN1isOK = aVprev.IsSame( aV1 );
+ aN2isOK = aVprev.IsSame( aV2 );
+ }
+ if ( !aN1isOK && !aN2isOK ) continue;
+ // 2. Collect parameters on the track edge
+ aPrms.clear();
+ aItN = locMeshDS->GetNodes();
+ while ( aItN->more() ) {
+ const SMDS_MeshNode* pNode = aItN->next();
+ const SMDS_EdgePosition* pEPos =
+ static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
+ double aT = pEPos->GetUParameter();
+ aPrms.push_back( aT );
+ }
+ list<SMESH_MeshEditor_PathPoint> LPP;
+ //Extrusion_Error err =
+ MakeEdgePathPoints(aPrms, aTrackEdge, aN1isOK, LPP);
+ LLPPs.push_back(LPP);
+ UsedNums.Add(k);
+ // update startN for search following egde
+ if ( aN1isOK ) aVprev = aV2;
+ else aVprev = aV1;
+ break;
+ }
+ }
+ list< list<SMESH_MeshEditor_PathPoint> >::iterator itLLPP = LLPPs.begin();
+ list<SMESH_MeshEditor_PathPoint>& firstList = *itLLPP;
+ fullList.splice( fullList.end(), firstList );
+
+ SMESH_MeshEditor_PathPoint PP1 = fullList.back();
+ fullList.pop_back();
+ itLLPP++;
+ for(; itLLPP!=LLPPs.end(); itLLPP++) {
+ list<SMESH_MeshEditor_PathPoint>& currList = *itLLPP;
+ SMESH_MeshEditor_PathPoint PP2 = currList.front();
+ gp_Dir D1 = PP1.Tangent();
+ gp_Dir D2 = PP2.Tangent();
+ gp_Dir Dnew( ( D1.XYZ() + D2.XYZ() ) / 2 );
+ PP1.SetTangent(Dnew);
+ fullList.push_back(PP1);
+ fullList.splice( fullList.end(), currList, ++currList.begin(), currList.end() );
+ PP1 = fullList.back();
+ fullList.pop_back();
+ }
+ // if wire not closed
+ fullList.push_back(PP1);
+ // else ???
+ }
+ else {
+ return EXTR_BAD_PATH_SHAPE;
+ }
+
+ return MakeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
+ theHasRefPoint, theRefPoint, theMakeGroups);
+}
+
+
+//=======================================================================
+//function : MakeEdgePathPoints
+//purpose : auxilary for ExtrusionAlongTrack
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::MakeEdgePathPoints(std::list<double>& aPrms,
+ const TopoDS_Edge& aTrackEdge,
+ bool FirstIsStart,
+ list<SMESH_MeshEditor_PathPoint>& LPP)
+{
+ Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2;
+ aTolVec=1.e-7;
+ aTolVec2=aTolVec*aTolVec;
+ double aT1, aT2;
+ TopoDS_Vertex aV1, aV2;
+ TopExp::Vertices( aTrackEdge, aV1, aV2 );
+ aT1=BRep_Tool::Parameter( aV1, aTrackEdge );
+ aT2=BRep_Tool::Parameter( aV2, aTrackEdge );
+ // 2. Collect parameters on the track edge
+ aPrms.push_front( aT1 );
+ aPrms.push_back( aT2 );
+ // sort parameters
+ aPrms.sort();
+ if( FirstIsStart ) {
+ if ( aT1 > aT2 ) {
+ aPrms.reverse();
+ }
+ }
+ else {
+ if ( aT2 > aT1 ) {
+ aPrms.reverse();
+ }
+ }
+ // 3. Path Points
+ SMESH_MeshEditor_PathPoint aPP;
+ Handle(Geom_Curve) aC3D = BRep_Tool::Curve( aTrackEdge, aTx1, aTx2 );
+ std::list<double>::iterator aItD = aPrms.begin();
+ for(; aItD != aPrms.end(); ++aItD) {
+ double aT = *aItD;
+ gp_Pnt aP3D;
+ gp_Vec aVec;
+ aC3D->D1( aT, aP3D, aVec );
+ aL2 = aVec.SquareMagnitude();
+ if ( aL2 < aTolVec2 )
+ return EXTR_CANT_GET_TANGENT;
+ gp_Dir aTgt( aVec );
+ aPP.SetPnt( aP3D );
+ aPP.SetTangent( aTgt );
+ aPP.SetParameter( aT );
+ LPP.push_back(aPP);
+ }
+ return EXTR_OK;
+}
+
+
+//=======================================================================
+//function : MakeExtrElements
+//purpose : auxilary for ExtrusionAlongTrack
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet theElemSets[2],
+ list<SMESH_MeshEditor_PathPoint>& fullList,
+ const bool theHasAngles,
+ list<double>& theAngles,
+ const bool theLinearVariation,
+ const bool theHasRefPoint,
+ const gp_Pnt& theRefPoint,
+ const bool theMakeGroups)
+{
+ const int aNbTP = fullList.size();
+ // Angles
+ if( theHasAngles && !theAngles.empty() && theLinearVariation )
+ LinearAngleVariation(aNbTP-1, theAngles);
+ // fill vector of path points with angles
+ vector<SMESH_MeshEditor_PathPoint> aPPs;
+ list<SMESH_MeshEditor_PathPoint>::iterator itPP = fullList.begin();
+ list<double>::iterator itAngles = theAngles.begin();
+ aPPs.push_back( *itPP++ );
+ for( ; itPP != fullList.end(); itPP++) {
+ aPPs.push_back( *itPP );
+ if ( theHasAngles && itAngles != theAngles.end() )
+ aPPs.back().SetAngle( *itAngles++ );
+ }
+
+ TNodeOfNodeListMap mapNewNodes;
+ TElemOfVecOfNnlmiMap mapElemNewNodes;
+ TTElemOfElemListMap newElemsMap;
+ TIDSortedElemSet::iterator itElem;
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+
+ // 3. Center of rotation aV0
+ gp_Pnt aV0 = theRefPoint;
+ if ( !theHasRefPoint )
+ {
+ gp_XYZ aGC( 0.,0.,0. );
+ TIDSortedElemSet newNodes;
+
+ for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+ {
+ TIDSortedElemSet& theElements = theElemSets[ is2ndSet ];
+ itElem = theElements.begin();
+ for ( ; itElem != theElements.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 )
+ aGC += SMESH_TNodeXYZ( node );
+ }
+ }
+ }
+ aGC /= newNodes.size();
+ aV0.SetXYZ( aGC );
+ } // if (!theHasRefPoint) {
+
+ // 4. Processing the elements
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ setElemsFirst( theElemSets );
+ for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+ {
+ TIDSortedElemSet& theElements = theElemSets[ is2ndSet ];
+ for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) {
+ // check element type
+ const SMDS_MeshElement* elem = *itElem;
+ if ( !elem )
+ continue;
+ // SMDSAbs_ElementType aTypeE = elem->GetType();
+ // if ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge )
+ // continue;
+
+ vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
+ newNodesItVec.reserve( elem->NbNodes() );
+
+ // loop on elem nodes
+ int nodeIndex = -1;
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() )
+ {
+ ++nodeIndex;
+ // check if a node has been already processed
+ const SMDS_MeshNode* node =
+ static_cast<const SMDS_MeshNode*>( itN->next() );
+ TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
+ if ( nIt == mapNewNodes.end() ) {
+ nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
+ list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
+
+ // make new nodes
+ Standard_Real aAngle1x, aAngleT1T0, aTolAng;
+ gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x;
+ gp_Ax1 anAx1, anAxT1T0;
+ gp_Dir aDT1x, aDT0x, aDT1T0;
+
+ aTolAng=1.e-4;
+
+ aV0x = aV0;
+ aPN0 = SMESH_TNodeXYZ( node );
+
+ const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
+ aP0x = aPP0.Pnt();
+ aDT0x= aPP0.Tangent();
+ //cout<<"j = 0 PP: Pnt("<<aP0x.X()<<","<<aP0x.Y()<<","<<aP0x.Z()<<")"<<endl;
+
+ for ( int j = 1; j < aNbTP; ++j ) {
+ const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j];
+ aP1x = aPP1.Pnt();
+ aDT1x = aPP1.Tangent();
+ aAngle1x = aPP1.Angle();
+
+ gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
+ // Translation
+ gp_Vec aV01x( aP0x, aP1x );
+ aTrsf.SetTranslation( aV01x );
+
+ // traslated point
+ aV1x = aV0x.Transformed( aTrsf );
+ aPN1 = aPN0.Transformed( aTrsf );
+
+ // rotation 1 [ T1,T0 ]
+ aAngleT1T0=-aDT1x.Angle( aDT0x );
+ if (fabs(aAngleT1T0) > aTolAng) {
+ aDT1T0=aDT1x^aDT0x;
+ anAxT1T0.SetLocation( aV1x );
+ anAxT1T0.SetDirection( aDT1T0 );
+ aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 );
+
+ aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
+ }
+
+ // rotation 2
+ if ( theHasAngles ) {
+ anAx1.SetLocation( aV1x );
+ anAx1.SetDirection( aDT1x );
+ aTrsfRot.SetRotation( anAx1, aAngle1x );
+
+ aPN1 = aPN1.Transformed( aTrsfRot );
+ }
+
+ // make new node
+ //MESSAGE("elem->IsQuadratic " << elem->IsQuadratic() << " " << elem->IsMediumNode(node));
+ if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
+ // create additional node
+ double x = ( aPN1.X() + aPN0.X() )/2.;
+ double y = ( aPN1.Y() + aPN0.Y() )/2.;
+ double z = ( aPN1.Z() + aPN0.Z() )/2.;
+ const SMDS_MeshNode* newNode = aMesh->AddNode(x,y,z);
+ myLastCreatedNodes.Append(newNode);
+ srcNodes.Append( node );
+ listNewNodes.push_back( newNode );
+ }
+ const SMDS_MeshNode* newNode = aMesh->AddNode( aPN1.X(), aPN1.Y(), aPN1.Z() );
+ myLastCreatedNodes.Append(newNode);
+ srcNodes.Append( node );
+ listNewNodes.push_back( newNode );
+
+ aPN0 = aPN1;
+ aP0x = aP1x;
+ aV0x = aV1x;
+ aDT0x = aDT1x;
+ }
+ }
+
+ else {
+ // if current elem is quadratic and current node is not medium
+ // we have to check - may be it is needed to insert additional nodes
+ if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
+ list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
+ if(listNewNodes.size()==aNbTP-1) {
+ vector<const SMDS_MeshNode*> aNodes(2*(aNbTP-1));
+ gp_XYZ P(node->X(), node->Y(), node->Z());
+ list< const SMDS_MeshNode* >::iterator it = listNewNodes.begin();
+ int i;
+ for(i=0; i<aNbTP-1; i++) {
+ const SMDS_MeshNode* N = *it;
+ double x = ( N->X() + P.X() )/2.;
+ double y = ( N->Y() + P.Y() )/2.;
+ double z = ( N->Z() + P.Z() )/2.;
+ const SMDS_MeshNode* newN = aMesh->AddNode(x,y,z);
+ srcNodes.Append( node );
+ myLastCreatedNodes.Append(newN);
+ aNodes[2*i] = newN;
+ aNodes[2*i+1] = N;
+ P = gp_XYZ(N->X(),N->Y(),N->Z());
+ }
+ listNewNodes.clear();
+ for(i=0; i<2*(aNbTP-1); i++) {
+ listNewNodes.push_back(aNodes[i]);
+ }
+ }
+ }
+ }
+
+ newNodesItVec.push_back( nIt );
+ }
+ // make new elements
+ //sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem],
+ // newNodesItVec[0]->second.size(), myLastCreatedElems );
+ sweepElement( elem, newNodesItVec, newElemsMap[elem], aNbTP-1, srcElems );
+ }
+ }
+
+ makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElemSets[0], aNbTP-1, srcElems );
+
+ if ( theMakeGroups )
+ generateGroups( srcNodes, srcElems, "extruded");
+
+ return EXTR_OK;
+}
+
+
+//=======================================================================
+//function : LinearAngleVariation
+//purpose : auxilary for ExtrusionAlongTrack
+//=======================================================================
+void SMESH_MeshEditor::LinearAngleVariation(const int nbSteps,
+ list<double>& Angles)
+{
+ int nbAngles = Angles.size();
+ if( nbSteps > nbAngles ) {
+ vector<double> theAngles(nbAngles);
+ list<double>::iterator it = Angles.begin();
+ int i = -1;
+ for(; it!=Angles.end(); it++) {
+ i++;
+ theAngles[i] = (*it);
+ }
+ 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.clear();
+ it = res.begin();
+ for(; it!=res.end(); it++)
+ Angles.push_back( *it );
+ }
+}
+
+
+//================================================================================
+/*!
+ * \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)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ bool needReverse = false;
+ string groupPostfix;
+ switch ( theTrsf.Form() ) {
+ case gp_PntMirror:
+ MESSAGE("gp_PntMirror");
+ needReverse = true;
+ groupPostfix = "mirrored";
+ break;
+ case gp_Ax1Mirror:
+ MESSAGE("gp_Ax1Mirror");
+ groupPostfix = "mirrored";
+ break;
+ case gp_Ax2Mirror:
+ MESSAGE("gp_Ax2Mirror");
+ needReverse = true;
+ groupPostfix = "mirrored";
+ break;
+ case gp_Rotation:
+ MESSAGE("gp_Rotation");
+ groupPostfix = "rotated";
+ break;
+ case gp_Translation:
+ MESSAGE("gp_Translation");
+ groupPostfix = "translated";
+ break;
+ case gp_Scale:
+ MESSAGE("gp_Scale");
+ groupPostfix = "scaled";
+ break;
+ case gp_CompoundTrsf: // different scale by axis
+ MESSAGE("gp_CompoundTrsf");
+ groupPostfix = "scaled";
+ break;
+ default:
+ MESSAGE("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.Append(newNode);
+ srcNodes.Append( node );
+ }
+ else if ( theCopy ) {
+ const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+ n2n_isnew.first->second = newNode;
+ myLastCreatedNodes.Append(newNode);
+ srcNodes.Append( 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();
+ int nbNodes = elem->NbNodes();
+ if ( geomType == SMDSGeom_NONE ) continue; // node
+
+ nodes.resize( nbNodes );
+
+ if ( geomType == SMDSGeom_POLYHEDRA ) // ------------------ polyhedral volume
+ {
+ const SMDS_VtkVolume* aPolyedre = dynamic_cast<const SMDS_VtkVolume*>( 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
+ int iNode = 0;
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() ) {
+ const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
+ TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
+ if ( nodeMapIt == nodeMap.end() )
+ break; // not all nodes transformed
+ nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
+ }
+ if ( iNode != nbNodes )
+ continue; // not all nodes transformed
+ }
+
+ if ( editor ) {
+ // copy in this or a new mesh
+ if ( editor->AddElement( nodes, elemType.Init( elem, /*basicOnly=*/false )))
+ srcElems.Append( 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 = editor->myLastCreatedElems;
+
+ PGroupIDs newGroupIDs;
+
+ if ( ( theMakeGroups && theCopy ) ||
+ ( theMakeGroups && theTargetMesh ) )
+ newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh, false );
+
+ 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 append to names of new groups
+ * \param targetMesh - mesh to create groups in
+ * \param topPresent - is there "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.Length() != elems.Length() )
+ throw SALOME_Exception("SMESH_MeshEditor::generateGroups(): invalid args");
+
+ // loop on created elements
+ for (int iElem = 1; iElem <= elems.Length(); ++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 < gens.Length() && 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 < gens.Length() && 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)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ 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(/*idInceasingOrder=*/true);
+ 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();
+
+ if (nbNodes < 3)
+ return 0;
+
+ set<const SMDS_MeshNode*> nodeSet;
+
+ // get simple seq of nodes
+ //const SMDS_MeshNode* simpleNodes[ nbNodes ];
+ vector<const SMDS_MeshNode*> simpleNodes( nbNodes );
+ int iSimple = 0, nbUnique = 0;
+
+ simpleNodes[iSimple++] = faceNodes[0];
+ nbUnique++;
+ for (int iCur = 1; iCur < nbNodes; iCur++) {
+ if (faceNodes[iCur] != simpleNodes[iSimple - 1]) {
+ simpleNodes[iSimple++] = faceNodes[iCur];
+ if (nodeSet.insert( faceNodes[iCur] ).second)
+ nbUnique++;
+ }
+ }
+ int nbSimple = iSimple;
+ if (simpleNodes[nbSimple - 1] == simpleNodes[0]) {
+ nbSimple--;
+ iSimple--;
+ }
+
+ if (nbUnique < 3)
+ return 0;
+
+ // separate loops
+ int nbNew = 0;
+ bool foundLoop = (nbSimple > nbUnique);
+ while (foundLoop) {
+ foundLoop = false;
+ set<const SMDS_MeshNode*> loopSet;
+ for (iSimple = 0; iSimple < nbSimple && !foundLoop; iSimple++) {
+ const SMDS_MeshNode* n = simpleNodes[iSimple];
+ if (!loopSet.insert( n ).second) {
+ foundLoop = true;
+
+ // separate loop
+ int iC = 0, curLast = iSimple;
+ for (; iC < curLast; iC++) {
+ if (simpleNodes[iC] == n) break;
+ }
+ int loopLen = curLast - iC;
+ if (loopLen > 2) {
+ // create sub-element
+ nbNew++;
+ quantities.push_back(loopLen);
+ for (; iC < curLast; iC++) {
+ poly_nodes.push_back(simpleNodes[iC]);
+ }
+ }
+ // shift the rest nodes (place from the first loop position)
+ for (iC = curLast + 1; iC < nbSimple; iC++) {
+ simpleNodes[iC - loopLen] = simpleNodes[iC];
+ }
+ nbSimple -= loopLen;
+ iSimple -= loopLen;
+ }
+ } // for (iSimple = 0; iSimple < nbSimple; iSimple++)
+ } // while (foundLoop)
+
+ if (iSimple > 2) {
+ nbNew++;
+ quantities.push_back(iSimple);
+ for (int i = 0; i < iSimple; i++)
+ poly_nodes.push_back(simpleNodes[i]);
+ }
+
+ return nbNew;
+}
+
+//=======================================================================
+//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)
+{
+ MESSAGE("MergeNodes");
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ TNodeNodeMap nodeNodeMap; // node to replace - new node
+ set<const SMDS_MeshElement*> elems; // all elements with changed nodes
+ list< int > rmElemIds, rmNodeIds;
+
+ // 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 ));
+ if ( nToRemove != nToKeep )
+ {
+ rmNodeIds.push_back( nToRemove->GetID() );
+ AddToSameGroups( nToKeep, nToRemove, aMesh );
+ // set _alwaysComputed to a sub-mesh of VERTEX to enable mesh computing
+ // after MergeNodes() w/o creating node in place of merged ones.
+ const SMDS_PositionPtr& pos = nToRemove->GetPosition();
+ if ( pos && pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
+ if ( SMESH_subMesh* sm = myMesh->GetSubMeshContaining( nToRemove->getshapeId() ))
+ sm->SetIsAlwaysComputed( true );
+ }
+ SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
+ while ( invElemIt->more() ) {
+ const SMDS_MeshElement* elem = invElemIt->next();
+ elems.insert(elem);
+ }
+ }
+ }
+ // Change element nodes or remove an element
+
+ set<const SMDS_MeshNode*> nodeSet;
+ vector< const SMDS_MeshNode*> curNodes, uniqueNodes;
+ vector<int> iRepl;
+ ElemFeatures elemType;
+
+ set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
+ for ( ; eIt != elems.end(); eIt++ )
+ {
+ const SMDS_MeshElement* elem = *eIt;
+ int nbNodes = elem->NbNodes();
+ int aShapeId = FindShape( elem );
+
+ nodeSet.clear();
+ 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 sticks
+ n = (*nnIt).second;
+ { ////////// BUG 0020185: begin
+ bool stopRecur = false;
+ set<const SMDS_MeshNode*> nodesRecur;
+ nodesRecur.insert(n);
+ while (!stopRecur) {
+ TNodeNodeMap::iterator nnIt_i = nodeNodeMap.find( n );
+ if ( nnIt_i != nodeNodeMap.end() ) { // n sticks
+ n = (*nnIt_i).second;
+ if (!nodesRecur.insert(n).second) {
+ // error: recursive dependancy
+ stopRecur = true;
+ }
+ }
+ else
+ stopRecur = true;
+ }
+ } ////////// BUG 0020185: end
+ }
+ curNodes[ iCur ] = n;
+ bool isUnique = nodeSet.insert( n ).second;
+ if ( isUnique )
+ uniqueNodes[ iUnique++ ] = n;
+ else
+ iRepl[ nbRepl++ ] = iCur;
+ iCur++;
+ }
+
+ // Analyse element topology after replacement
+
+ bool isOk = true;
+ int nbUniqueNodes = nodeSet.size();
+ if ( nbNodes != nbUniqueNodes ) // some nodes stick
+ {
+ if (elem->IsPoly()) // Polygons and Polyhedral volumes
+ {
+ if (elem->GetType() == SMDSAbs_Face) // Polygon
+ {
+ elemType.Init( elem );
+ 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;
+ int nbNew = SimplifyFace( curNodes, polygons_nodes, quantities );
+ if (nbNew > 0)
+ {
+ vector<const SMDS_MeshNode *> face_nodes;
+ int inode = 0;
+ for (int iface = 0; iface < nbNew; iface++)
+ {
+ 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 );
+ }
+ SMDS_MeshElement* newElem = AddElement( face_nodes, elemType );
+ if ( aShapeId )
+ aMesh->SetMeshElementOnShape(newElem, aShapeId);
+ }
+ }
+ rmElemIds.push_back(elem->GetID());
+
+ } // Polygon
+
+ else if (elem->GetType() == SMDSAbs_Volume) // Polyhedral volume
+ {
+ if (nbUniqueNodes < 4) {
+ rmElemIds.push_back(elem->GetID());
+ }
+ else {
+ // each face has to be analyzed in order to check volume validity
+ const SMDS_VtkVolume* aPolyedre =
+ dynamic_cast<const SMDS_VtkVolume*>( elem );
+ if (aPolyedre) {
+ int nbFaces = aPolyedre->NbFaces();
+
+ vector<const SMDS_MeshNode *> poly_nodes;
+ vector<int> quantities;
+
+ for (int iface = 1; iface <= nbFaces; iface++) {
+ int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
+ vector<const SMDS_MeshNode *> faceNodes (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) {
+ // to be done: remove coincident faces
+ }
+
+ if (quantities.size() > 3)
+ {
+ //aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
+ const SMDS_MeshElement* newElem =
+ aMesh->AddPolyhedralVolume(poly_nodes, quantities);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ rmElemIds.push_back(elem->GetID());
+ }
+ }
+ else {
+ rmElemIds.push_back(elem->GetID());
+ }
+ }
+ }
+ else {
+ }
+
+ continue;
+ } // poly element
+
+ // Regular elements
+ // TODO not all the possible cases are solved. Find something more generic?
+ switch ( nbNodes ) {
+ case 2: ///////////////////////////////////// EDGE
+ isOk = false; break;
+ case 3: ///////////////////////////////////// TRIANGLE
+ isOk = false; break;
+ case 4:
+ if ( elem->GetType() == SMDSAbs_Volume ) // TETRAHEDRON
+ isOk = false;
+ else { //////////////////////////////////// QUADRANGLE
+ if ( nbUniqueNodes < 3 )
+ isOk = false;
+ else if ( nbRepl == 2 && iRepl[ 1 ] - iRepl[ 0 ] == 2 )
+ isOk = false; // opposite nodes stick
+ //MESSAGE("isOk " << isOk);
+ }
+ break;
+ case 6: ///////////////////////////////////// PENTAHEDRON
+ if ( nbUniqueNodes == 4 ) {
+ // ---------------------------------> tetrahedron
+ if (nbRepl == 3 &&
+ iRepl[ 0 ] > 2 && iRepl[ 1 ] > 2 && iRepl[ 2 ] > 2 ) {
+ // all top nodes stick: reverse a bottom
+ uniqueNodes[ 0 ] = curNodes [ 1 ];
+ uniqueNodes[ 1 ] = curNodes [ 0 ];
+ }
+ else if (nbRepl == 3 &&
+ iRepl[ 0 ] < 3 && iRepl[ 1 ] < 3 && iRepl[ 2 ] < 3 ) {
+ // all bottom nodes stick: set a top before
+ uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
+ uniqueNodes[ 0 ] = curNodes [ 3 ];
+ uniqueNodes[ 1 ] = curNodes [ 4 ];
+ uniqueNodes[ 2 ] = curNodes [ 5 ];
+ }
+ else if (nbRepl == 4 &&
+ iRepl[ 2 ] - iRepl [ 0 ] == 3 && iRepl[ 3 ] - iRepl [ 1 ] == 3 ) {
+ // a lateral face turns into a line: reverse a bottom
+ uniqueNodes[ 0 ] = curNodes [ 1 ];
+ uniqueNodes[ 1 ] = curNodes [ 0 ];
+ }
+ else
+ isOk = false;
+ }
+ else if ( nbUniqueNodes == 5 ) {
+ // PENTAHEDRON --------------------> 2 tetrahedrons
+ if ( nbRepl == 2 && iRepl[ 1 ] - iRepl [ 0 ] == 3 ) {
+ // a bottom node sticks with a linked top one
+ // 1.
+ SMDS_MeshElement* newElem =
+ aMesh->AddVolume(curNodes[ 3 ],
+ curNodes[ 4 ],
+ curNodes[ 5 ],
+ curNodes[ iRepl[ 0 ] == 2 ? 1 : 2 ]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ // 2. : reverse a bottom
+ uniqueNodes[ 0 ] = curNodes [ 1 ];
+ uniqueNodes[ 1 ] = curNodes [ 0 ];
+ nbUniqueNodes = 4;
+ }
+ else
+ isOk = false;
+ }
+ else
+ isOk = false;
+ break;
+ case 8: {
+ if(elem->IsQuadratic()) { // Quadratic quadrangle
+ // 1 5 2
+ // +---+---+
+ // | |
+ // | |
+ // 4+ +6
+ // | |
+ // | |
+ // +---+---+
+ // 0 7 3
+ isOk = false;
+ if(nbRepl==2) {
+ MESSAGE("nbRepl=2: " << iRepl[0] << " " << iRepl[1]);
+ }
+ if(nbRepl==3) {
+ MESSAGE("nbRepl=3: " << iRepl[0] << " " << iRepl[1] << " " << iRepl[2]);
+ nbUniqueNodes = 6;
+ if( iRepl[0]==0 && iRepl[1]==1 && iRepl[2]==4 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[2];
+ uniqueNodes[2] = curNodes[3];
+ uniqueNodes[3] = curNodes[5];
+ uniqueNodes[4] = curNodes[6];
+ uniqueNodes[5] = curNodes[7];
+ isOk = true;
+ }
+ if( iRepl[0]==0 && iRepl[1]==3 && iRepl[2]==7 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[1];
+ uniqueNodes[2] = curNodes[2];
+ uniqueNodes[3] = curNodes[4];
+ uniqueNodes[4] = curNodes[5];
+ uniqueNodes[5] = curNodes[6];
+ isOk = true;
+ }
+ if( iRepl[0]==0 && iRepl[1]==4 && iRepl[2]==7 ) {
+ uniqueNodes[0] = curNodes[1];
+ uniqueNodes[1] = curNodes[2];
+ uniqueNodes[2] = curNodes[3];
+ uniqueNodes[3] = curNodes[5];
+ uniqueNodes[4] = curNodes[6];
+ uniqueNodes[5] = curNodes[0];
+ isOk = true;
+ }
+ if( iRepl[0]==1 && iRepl[1]==2 && iRepl[2]==5 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[1];
+ uniqueNodes[2] = curNodes[3];
+ uniqueNodes[3] = curNodes[4];
+ uniqueNodes[4] = curNodes[6];
+ uniqueNodes[5] = curNodes[7];
+ isOk = true;
+ }
+ if( iRepl[0]==1 && iRepl[1]==4 && iRepl[2]==5 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[2];
+ uniqueNodes[2] = curNodes[3];
+ uniqueNodes[3] = curNodes[1];
+ uniqueNodes[4] = curNodes[6];
+ uniqueNodes[5] = curNodes[7];
+ isOk = true;
+ }
+ if( iRepl[0]==2 && iRepl[1]==3 && iRepl[2]==6 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[1];
+ uniqueNodes[2] = curNodes[2];
+ uniqueNodes[3] = curNodes[4];
+ uniqueNodes[4] = curNodes[5];
+ uniqueNodes[5] = curNodes[7];
+ isOk = true;
+ }
+ if( iRepl[0]==2 && iRepl[1]==5 && iRepl[2]==6 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[1];
+ uniqueNodes[2] = curNodes[3];
+ uniqueNodes[3] = curNodes[4];
+ uniqueNodes[4] = curNodes[2];
+ uniqueNodes[5] = curNodes[7];
+ isOk = true;
+ }
+ if( iRepl[0]==3 && iRepl[1]==6 && iRepl[2]==7 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[1];
+ uniqueNodes[2] = curNodes[2];
+ uniqueNodes[3] = curNodes[4];
+ uniqueNodes[4] = curNodes[5];
+ uniqueNodes[5] = curNodes[3];
+ isOk = true;
+ }
+ }
+ if(nbRepl==4) {
+ MESSAGE("nbRepl=4: " << iRepl[0] << " " << iRepl[1] << " " << iRepl[2] << " " << iRepl[3]);
+ }
+ if(nbRepl==5) {
+ MESSAGE("nbRepl=5: " << iRepl[0] << " " << iRepl[1] << " " << iRepl[2] << " " << iRepl[3] << " " << iRepl[4]);
+ }
+ break;
+ }
+ //////////////////////////////////// HEXAHEDRON
+ isOk = false;
+ SMDS_VolumeTool hexa (elem);
+ hexa.SetExternalNormal();
+ if ( nbUniqueNodes == 4 && nbRepl == 4 ) {
+ //////////////////////// HEX ---> 1 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 iOppFace = hexa.GetOppFaceIndex( iFace );
+ ind = hexa.GetFaceNodesIndices( iOppFace );
+ int nbStick = 0;
+ for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
+ if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
+ nbStick++;
+ }
+ if ( nbStick == 1 ) {
+ // ... and the opposite one - into a triangle.
+ // set a top node
+ ind = hexa.GetFaceNodesIndices( iFace );
+ uniqueNodes[ 3 ] = curNodes[ind[ 0 ]];
+ isOk = true;
+ }
+ break;
+ }
+ }
+ }
+ else if ( nbUniqueNodes == 6 && nbRepl == 2 ) {
+ //////////////////////// HEX ---> 1 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] ))
+ {
+ isOk = true;
+ // 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;
+ }
+ }
+ break;
+ }
+ else if (nbUniqueNodes == 5 && nbRepl == 4 ) {
+ //////////////////// HEXAHEDRON ---> 2 tetrahedrons
+ 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 );
+ int nbStick = 0;
+ iUnique = 2; // reverse a tetrahedron 1 bottom
+ for ( iCur = 0; iCur < 4 && nbStick == 0; iCur++ ) {
+ if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
+ nbStick++;
+ else if ( iUnique >= 0 )
+ uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
+ }
+ if ( nbStick == 0 ) {
+ // ... and the opposite one is a quadrangle
+ // set a top node
+ const int* indTop = hexa.GetFaceNodesIndices( iFace );
+ uniqueNodes[ 3 ] = curNodes[indTop[ 0 ]];
+ nbUniqueNodes = 4;
+ // tetrahedron 2
+ SMDS_MeshElement* newElem =
+ aMesh->AddVolume(curNodes[ind[ 0 ]],
+ curNodes[ind[ 3 ]],
+ curNodes[ind[ 2 ]],
+ curNodes[indTop[ 0 ]]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ isOk = true;
+ }
+ break;
+ }
+ }
+ }
+ else if ( nbUniqueNodes == 6 && nbRepl == 4 ) {
+ ////////////////// HEXAHEDRON ---> 2 tetrahedrons or 1 prism
+ // find indices of quad and tri faces
+ int iQuadFace[ 6 ], iTriFace[ 6 ], nbQuad = 0, nbTri = 0, iFace;
+ for ( iFace = 0; iFace < 6; iFace++ ) {
+ const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
+ nodeSet.clear();
+ for ( iCur = 0; iCur < 4; iCur++ )
+ nodeSet.insert( curNodes[ind[ iCur ]] );
+ nbUniqueNodes = nodeSet.size();
+ if ( nbUniqueNodes == 3 )
+ iTriFace[ nbTri++ ] = iFace;
+ else if ( nbUniqueNodes == 4 )
+ iQuadFace[ nbQuad++ ] = iFace;
+ }
+ if (nbQuad == 2 && nbTri == 4 &&
+ hexa.GetOppFaceIndex( iQuadFace[ 0 ] ) == iQuadFace[ 1 ]) {
+ // 2 opposite quadrangles stuck with a diagonal;
+ // sample groups of merged indices: (0-4)(2-6)
+ // --------------------------------------------> 2 tetrahedrons
+ const int *ind1 = hexa.GetFaceNodesIndices( iQuadFace[ 0 ]); // indices of quad1 nodes
+ const int *ind2 = hexa.GetFaceNodesIndices( iQuadFace[ 1 ]);
+ int i0, i1d, i2, i3d, i0t, i2t; // d-daigonal, t-top
+ if (curNodes[ind1[ 0 ]] == curNodes[ind2[ 0 ]] &&
+ curNodes[ind1[ 2 ]] == curNodes[ind2[ 2 ]]) {
+ // stuck with 0-2 diagonal
+ i0 = ind1[ 3 ];
+ i1d = ind1[ 0 ];
+ i2 = ind1[ 1 ];
+ i3d = ind1[ 2 ];
+ i0t = ind2[ 1 ];
+ i2t = ind2[ 3 ];
+ }
+ else if (curNodes[ind1[ 1 ]] == curNodes[ind2[ 3 ]] &&
+ curNodes[ind1[ 3 ]] == curNodes[ind2[ 1 ]]) {
+ // stuck with 1-3 diagonal
+ i0 = ind1[ 0 ];
+ i1d = ind1[ 1 ];
+ i2 = ind1[ 2 ];
+ i3d = ind1[ 3 ];
+ i0t = ind2[ 0 ];
+ i2t = ind2[ 1 ];
+ }
+ else {
+ ASSERT(0);
+ }
+ // tetrahedron 1
+ uniqueNodes[ 0 ] = curNodes [ i0 ];
+ uniqueNodes[ 1 ] = curNodes [ i1d ];
+ uniqueNodes[ 2 ] = curNodes [ i3d ];
+ uniqueNodes[ 3 ] = curNodes [ i0t ];
+ nbUniqueNodes = 4;
+ // tetrahedron 2
+ SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ i1d ],
+ curNodes[ i2 ],
+ curNodes[ i3d ],
+ curNodes[ i2t ]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ isOk = true;
+ }
+ else if (( nbTri == 2 && nbQuad == 3 ) || // merged (0-4)(1-5)
+ ( nbTri == 4 && nbQuad == 2 )) { // merged (7-4)(1-5)
+ // --------------------------------------------> prism
+ // find 2 opposite triangles
+ nbUniqueNodes = 6;
+ for ( iFace = 0; iFace + 1 < nbTri; iFace++ ) {
+ if ( hexa.GetOppFaceIndex( iTriFace[ iFace ] ) == iTriFace[ iFace + 1 ]) {
+ // find indices of kept and replaced nodes
+ // and fill unique nodes of 2 opposite triangles
+ const int *ind1 = hexa.GetFaceNodesIndices( iTriFace[ iFace ]);
+ const int *ind2 = hexa.GetFaceNodesIndices( iTriFace[ iFace + 1 ]);
+ const SMDS_MeshNode** hexanodes = hexa.GetNodes();
+ // fill unique nodes
+ iUnique = 0;
+ isOk = true;
+ for ( iCur = 0; iCur < 4 && isOk; iCur++ ) {
+ const SMDS_MeshNode* n = curNodes[ind1[ iCur ]];
+ const SMDS_MeshNode* nInit = hexanodes[ind1[ iCur ]];
+ if ( n == nInit ) {
+ // iCur of a linked node of the opposite face (make normals co-directed):
+ int iCurOpp = ( iCur == 1 || iCur == 3 ) ? 4 - iCur : iCur;
+ // check that correspondent corners of triangles are linked
+ if ( !hexa.IsLinked( ind1[ iCur ], ind2[ iCurOpp ] ))
+ isOk = false;
+ else {
+ uniqueNodes[ iUnique ] = n;
+ uniqueNodes[ iUnique + 3 ] = curNodes[ind2[ iCurOpp ]];
+ iUnique++;
+ }
+ }
+ }
+ break;
+ }
+ }
+ }
+ } // if ( nbUniqueNodes == 6 && nbRepl == 4 )
+ else
+ {
+ MESSAGE("MergeNodes() removes hexahedron "<< elem);
+ }
+ break;
+ } // HEXAHEDRON
+
+ default:
+ isOk = false;
+ } // switch ( nbNodes )
+
+ } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
+
+ if ( isOk ) // the non-poly elem remains valid after sticking nodes
+ {
+ elemType.Init( elem ).SetID( elem->GetID() );
+
+ SMESHDS_SubMesh * sm = aShapeId > 0 ? aMesh->MeshElements(aShapeId) : 0;
+ aMesh->RemoveFreeElement(elem, sm, /*fromGroups=*/false);
+
+ uniqueNodes.resize(nbUniqueNodes);
+ SMDS_MeshElement* newElem = this->AddElement( uniqueNodes, elemType );
+ if ( sm && newElem )
+ sm->AddElement( newElem );
+ if ( elem != newElem )
+ ReplaceElemInGroups( elem, newElem, aMesh );
+ }
+ else {
+ // Remove invalid regular element or invalid polygon
+ rmElemIds.push_back( elem->GetID() );
+ }
+
+ } // loop on elements
+
+ // Remove bad elements, then equal nodes (order important)
+
+ Remove( rmElemIds, false );
+ Remove( rmNodeIds, true );
+
+ return;
+}
+
+
+// ========================================================
+// class : SortableElement
+// purpose : allow sorting elements basing on their nodes
+// ========================================================
+class SortableElement : public set <const SMDS_MeshElement*>
+{
+public:
+
+ SortableElement( const SMDS_MeshElement* theElem )
+ {
+ myElem = theElem;
+ SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
+ while ( nodeIt->more() )
+ this->insert( nodeIt->next() );
+ }
+
+ const SMDS_MeshElement* Get() const
+ { return myElem; }
+
+ void Set(const SMDS_MeshElement* e) const
+ { myElem = e; }
+
+
+private:
+ mutable const SMDS_MeshElement* myElem;
+};
+
+//=======================================================================
+//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)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ typedef map< SortableElement, int > TMapOfNodeSet;
+ typedef list<int> TGroupOfElems;
+
+ if ( theElements.empty() )
+ { // get all elements in the mesh
+ SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
+ while ( eIt->more() )
+ theElements.insert( theElements.end(), eIt->next() );
+ }
+
+ vector< TGroupOfElems > arrayOfGroups;
+ TGroupOfElems groupOfElems;
+ TMapOfNodeSet mapOfNodeSet;
+
+ TIDSortedElemSet::iterator elemIt = theElements.begin();
+ for ( int i = 0; elemIt != theElements.end(); ++elemIt )
+ {
+ const SMDS_MeshElement* curElem = *elemIt;
+ SortableElement SE(curElem);
+ // check uniqueness
+ pair< TMapOfNodeSet::iterator, bool> pp = mapOfNodeSet.insert(make_pair(SE, i));
+ if ( !pp.second ) { // one more coincident elem
+ TMapOfNodeSet::iterator& itSE = pp.first;
+ int ind = (*itSE).second;
+ arrayOfGroups[ind].push_back( curElem->GetID() );
+ }
+ else {
+ arrayOfGroups.push_back( groupOfElems );
+ arrayOfGroups.back().push_back( curElem->GetID() );
+ i++;
+ }
+ }
+
+ groupOfElems.clear();
+ vector< TGroupOfElems >::iterator groupIt = arrayOfGroups.begin();
+ for ( ; groupIt != arrayOfGroups.end(); ++groupIt )
+ {
+ if ( groupIt->size() > 1 ) {
+ //groupOfElems.sort(); -- theElements is sorted already
+ theGroupsOfElementsID.push_back( groupOfElems );
+ theGroupsOfElementsID.back().splice( theGroupsOfElementsID.back().end(), *groupIt );
+ }
+ }
+}
+
+//=======================================================================
+//function : MergeElements
+//purpose : In each given group, substitute all elements by the first one.
+//=======================================================================
+
+void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElementsID)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ 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 : 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 );
+
+ //vector<const SMDS_MeshNode*> nodes;
+ const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
+ TIDSortedElemSet foundElems;
+ bool needTheLast = ( theLastNode != 0 );
+
+ while ( nStart != theLastNode ) {
+ if ( nStart == theFirstNode )
+ return !needTheLast;
+
+ // find all free border faces sharing form 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 ) {
+ // get nodes
+ int iNode = 0, nbNodes = e->NbNodes();
+ //const SMDS_MeshNode* nodes[nbNodes+1];
+ vector<const SMDS_MeshNode*> nodes(nbNodes+1);
+
+ if(e->IsQuadratic()) {
+ const SMDS_VtkFace* F =
+ dynamic_cast<const SMDS_VtkFace*>(e);
+ if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
+ // use special nodes iterator
+ SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
+ while( anIter->more() ) {
+ nodes[ iNode++ ] = cast2Node(anIter->next());
+ }
+ }
+ else {
+ SMDS_ElemIteratorPtr nIt = e->nodesIterator();
+ while ( nIt->more() )
+ nodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
+ }
+ nodes[ iNode ] = nodes[ 0 ];
+ // check 2 links
+ 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 ? 1 : 0 )]);
+ 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;
+
+ // append 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
+ list< const SMDS_MeshNode* >::iterator nIt = cNL->begin();
+ list< const SMDS_MeshElement* >::iterator fIt = cFL->begin();
+ for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt );
+ for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt );
+ 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 :
+//=======================================================================
+
+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)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ MESSAGE("::SewFreeBorder()");
+ 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 wont 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;
+ 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, 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
+ if(elem->IsQuadratic()) {
+ const SMDS_VtkFace* F =
+ dynamic_cast<const SMDS_VtkFace*>(elem);
+ if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
+ // use special nodes iterator
+ SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
+ while( anIter->more() ) {
+ nodes[ iNode ] = cast2Node(anIter->next());
+ if ( nodes[ iNode++ ] == prevSideNode )
+ iPrevNode = iNode - 1;
+ }
+ }
+ else {
+ SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
+ 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(" Cant 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[] = { nSide[0].size(), nSide[1].size() };
+ int maxNbNodes = Max( nbNodes[0], nbNodes[1] );
+
+ TListOfListOfNodes nodeGroupsToMerge;
+ if ( nbNodes[0] == nbNodes[1] ||
+ ( 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
+ //double param[ 2 ][ maxNbNodes ];
+ double* param[ 2 ];
+ param[0] = new double [ maxNbNodes ];
+ param[1] = new double [ 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();
+
+ 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
+ gp_XYZ p1 (n1->X(), n1->Y(), n1->Z());
+ gp_XYZ p2 (n2->X(), n2->Y(), n2->Z());
+ 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() );
+ }
+ 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 (true) {
+ const SMDS_MeshElement* adjElem = findAdjacentFace( n12, n22, elem );
+ if ( adjElem )
+ InsertNodesIntoLink( adjElem, n12, n22, nodeList, toCreatePolygons );
+ else
+ break;
+ }
+ 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( elem );
+ elem = findAdjacentFace( n1, n2, 0 );
+ }
+ if ( notFound || otherLink ) {
+ // add element and nodes of the side into the insertMap
+ insertMapIt = insertMap.insert
+ ( TElemOfNodeListMap::value_type( 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;
+ const SMDS_MeshNode* n1 = nodeList.front(); nodeList.pop_front();
+ const SMDS_MeshNode* n2 = nodeList.front(); nodeList.pop_front();
+
+ InsertNodesIntoLink( elem, n1, n2, nodeList, toCreatePolygons );
+
+ if ( !theSideIsFreeBorder ) {
+ // look for and insert nodes into the faces adjacent to elem
+ while (true) {
+ const SMDS_MeshElement* adjElem = findAdjacentFace( n1, n2, elem );
+ if ( adjElem )
+ InsertNodesIntoLink( adjElem, n1, n2, nodeList, toCreatePolygons );
+ else
+ break;
+ }
+ }
+ if (toCreatePolyedrs) {
+ // perform insertion into the links of adjacent volumes
+ UpdateVolumes(n1, n2, nodeList);
+ }
+ }
+
+ delete param[0];
+ delete param[1];
+ } // end: insert new nodes
+
+ MergeNodes ( nodeGroupsToMerge );
+
+ return aResult;
+}
+
+//=======================================================================
+//function : InsertNodesIntoLink
+//purpose : insert theNodesToInsert into theFace between theBetweenNode1
+// and theBetweenNode2 and split theElement
+//=======================================================================
+
+void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theFace,
+ const SMDS_MeshNode* theBetweenNode1,
+ const SMDS_MeshNode* theBetweenNode2,
+ list<const SMDS_MeshNode*>& theNodesToInsert,
+ const bool toCreatePoly)
+{
+ 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;
+ //const SMDS_MeshNode* nodes[ theFace->NbNodes() ];
+ vector<const SMDS_MeshNode*> nodes( theFace->NbNodes() );
+
+ if(theFace->IsQuadratic()) {
+ const SMDS_VtkFace* F =
+ dynamic_cast<const SMDS_VtkFace*>(theFace);
+ if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
+ // use special nodes iterator
+ SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
+ while( anIter->more() ) {
+ const SMDS_MeshNode* n = cast2Node(anIter->next());
+ if ( n == theBetweenNode1 )
+ il1 = iNode;
+ else if ( n == theBetweenNode2 )
+ il2 = iNode;
+ else if ( i3 < 0 )
+ i3 = iNode;
+ else
+ i4 = iNode;
+ nodes[ iNode++ ] = n;
+ }
+ }
+ else {
+ SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
+ while ( nodeIt->more() ) {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( 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 ;