+ else if ( nbn == 6 && isQuadratic ) /////// quadratic triangle
+ {
+ f = aMesh->FindFace( nodes[0], nodes[2], nodes[4], nodes[1], nodes[3], nodes[5] );
+ if ( !f ||
+ nodes[2] != f->GetNodeWrap( f->GetNodeIndex( nodes[0] ) + 2*nextShift ))
+ {
+ const SMDS_MeshNode* newOrder[6] = { nodes[2 - 2*nextShift],
+ nodes[2],
+ nodes[2 + 2*nextShift],
+ nodes[3 - 2*nextShift],
+ nodes[3],
+ nodes[3 + 2*nextShift]};
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace( newOrder[ 0 ],
+ newOrder[ 1 ],
+ newOrder[ 2 ],
+ newOrder[ 3 ],
+ newOrder[ 4 ],
+ newOrder[ 5 ] ));
+ }
+ }
+ else if ( nbn == 8 && isQuadratic ) /////// quadratic quadrangle
+ {
+ f = aMesh->FindFace( nodes[0], nodes[2], nodes[4], nodes[6],
+ nodes[1], nodes[3], nodes[5], nodes[7] );
+ if ( !f ||
+ nodes[ 2 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 2*nextShift ))
+ {
+ const SMDS_MeshNode* newOrder[8] = { nodes[0],
+ nodes[4 - 2*nextShift],
+ nodes[4],
+ nodes[4 + 2*nextShift],
+ nodes[1],
+ nodes[5 - 2*nextShift],
+ nodes[5],
+ nodes[5 + 2*nextShift] };
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace(newOrder[ 0 ], newOrder[ 1 ],
+ newOrder[ 2 ], newOrder[ 3 ],
+ newOrder[ 4 ], newOrder[ 5 ],
+ newOrder[ 6 ], newOrder[ 7 ]));
+ }
+ }
+ else if ( nbn == 9 && isQuadratic ) /////// bi-quadratic quadrangle
+ {
+ f = aMesh->FindElement( vector<const SMDS_MeshNode*>( nodes, nodes+nbn ),
+ SMDSAbs_Face, /*noMedium=*/false);
+ if ( !f ||
+ nodes[ 2 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 2*nextShift ))
+ {
+ const SMDS_MeshNode* newOrder[9] = { nodes[0],
+ nodes[4 - 2*nextShift],
+ nodes[4],
+ nodes[4 + 2*nextShift],
+ nodes[1],
+ nodes[5 - 2*nextShift],
+ nodes[5],
+ nodes[5 + 2*nextShift],
+ nodes[8] };
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace(newOrder[ 0 ], newOrder[ 1 ],
+ newOrder[ 2 ], newOrder[ 3 ],
+ newOrder[ 4 ], newOrder[ 5 ],
+ newOrder[ 6 ], newOrder[ 7 ],
+ newOrder[ 8 ]));
+ }
+ }
+ else //////// polygon
+ {
+ vector<const SMDS_MeshNode*> polygon_nodes ( nodes, nodes+nbn );
+ const SMDS_MeshFace * f = aMesh->FindFace( polygon_nodes );
+ if ( !f ||
+ nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + nextShift ))
+ {
+ if ( !vTool.IsForward() )
+ std::reverse( polygon_nodes.begin(), polygon_nodes.end());
+ if ( f )
+ aMesh->ChangeElementNodes( f, &polygon_nodes[0], nbn );
+ else
+ AddElement(polygon_nodes, SMDSAbs_Face, polygon_nodes.size()>4);
+ }
+ }
+
+ while ( srcElements.Length() < myLastCreatedElems.Length() )
+ srcElements.Append( *srcEdge );
+
+ } // loop on free faces
+
+ // go to the next volume
+ iVol = 0;
+ while ( iVol++ < nbVolumesByStep ) v++;
+
+ } // loop on steps
+ } // loop on volumes of one step
+ } // sweep free links into faces
+
+ // Make a ceiling face with a normal external to a volume
+
+ SMDS_VolumeTool lastVol( itElem->second.back(), /*ignoreCentralNodes=*/false );
+
+ int iF = lastVol.GetFaceIndex( aFaceLastNodes );
+ if ( iF >= 0 ) {
+ lastVol.SetExternalNormal();
+ const SMDS_MeshNode** nodes = lastVol.GetFaceNodes( iF );
+ int nbn = lastVol.NbFaceNodes( iF );
+ if ( nbn == 3 ) {
+ if (!hasFreeLinks ||
+ !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]))
+ myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
+ }
+ else if ( nbn == 4 )
+ {
+ if (!hasFreeLinks ||
+ !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]))
+ myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]));
+ }
+ else if ( nbn == 6 && isQuadratic )
+ {
+ if (!hasFreeLinks ||
+ !aMesh->FindFace(nodes[0], nodes[2], nodes[4], nodes[1], nodes[3], nodes[5]) )
+ myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
+ nodes[1], nodes[3], nodes[5]));
+ }
+ else if ( nbn == 8 && isQuadratic )
+ {
+ if (!hasFreeLinks ||
+ !aMesh->FindFace(nodes[0], nodes[2], nodes[4], nodes[6],
+ nodes[1], nodes[3], nodes[5], nodes[7]) )
+ myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
+ nodes[1], nodes[3], nodes[5], nodes[7]));
+ }
+ else if ( nbn == 9 && isQuadratic )
+ {
+ if (!hasFreeLinks ||
+ !aMesh->FindElement(vector<const SMDS_MeshNode*>( nodes, nodes+nbn ),
+ SMDSAbs_Face, /*noMedium=*/false) )
+ myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
+ nodes[1], nodes[3], nodes[5], nodes[7],
+ nodes[8]));
+ }
+ else {
+ vector<const SMDS_MeshNode*> polygon_nodes ( nodes, nodes + nbn );
+ if (!hasFreeLinks || !aMesh->FindFace(polygon_nodes))
+ myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
+ }
+
+ while ( srcElements.Length() < myLastCreatedElems.Length() )
+ srcElements.Append( myLastCreatedElems.Last() );
+ }
+ } // loop on swept elements
+}
+
+//=======================================================================
+//function : RotationSweep
+//purpose :
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
+ const gp_Ax1& theAxis,
+ const double theAngle,
+ const int theNbSteps,
+ const double theTol,
+ const bool theMakeGroups,
+ const bool theMakeWalls)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+
+ MESSAGE( "RotationSweep()");
+ 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;
+ TElemOfElemListMap newElemsMap;
+
+ const bool isQuadraticMesh = bool( myMesh->NbEdges(ORDER_QUADRATIC) +
+ myMesh->NbFaces(ORDER_QUADRATIC) +
+ myMesh->NbVolumes(ORDER_QUADRATIC) );
+ // loop on theElems
+ TIDSortedElemSet::iterator itElem;
+ 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() )
+ {
+ // check if a node has been already sweeped
+ 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 );
+
+ 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.Append(newNode);
+ srcNodes.Append( 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.Append(newNode);
+ srcNodes.Append( 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, theElems, theNbSteps, srcElems );
+
+ PGroupIDs newGroupIDs;
+ if ( theMakeGroups )
+ newGroupIDs = generateGroups( srcNodes, srcElems, "rotated");
+
+ return newGroupIDs;
+}
+
+
+//=======================================================================
+//function : CreateNode
+//purpose :
+//=======================================================================
+const SMDS_MeshNode* SMESH_MeshEditor::CreateNode(const double x,
+ const double y,
+ const double z,
+ const double tolnode,
+ SMESH_SequenceOfNode& aNodes)
+{
+ // myLastCreatedElems.Clear();
+ // myLastCreatedNodes.Clear();
+
+ gp_Pnt P1(x,y,z);
+ SMESHDS_Mesh * aMesh = myMesh->GetMeshDS();
+
+ // try to search in sequence of existing nodes
+ // if aNodes.Length()>0 we 'nave to use given sequence
+ // else - use all nodes of mesh
+ if(aNodes.Length()>0) {
+ int i;
+ for(i=1; i<=aNodes.Length(); i++) {
+ gp_Pnt P2(aNodes.Value(i)->X(),aNodes.Value(i)->Y(),aNodes.Value(i)->Z());
+ if(P1.Distance(P2)<tolnode)
+ return aNodes.Value(i);
+ }
+ }
+ else {
+ SMDS_NodeIteratorPtr itn = aMesh->nodesIterator();
+ while(itn->more()) {
+ const SMDS_MeshNode* aN = static_cast<const SMDS_MeshNode*> (itn->next());
+ gp_Pnt P2(aN->X(),aN->Y(),aN->Z());
+ if(P1.Distance(P2)<tolnode)
+ return aN;
+ }
+ }
+
+ // create new node and return it
+ const SMDS_MeshNode* NewNode = aMesh->AddNode(x,y,z);
+ //myLastCreatedNodes.Append(NewNode);
+ return NewNode;
+}
+
+
+//=======================================================================
+//function : ExtrusionSweep
+//purpose :
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems,
+ const gp_Vec& theStep,
+ const int theNbSteps,
+ TElemOfElemListMap& newElemsMap,
+ const bool theMakeGroups,
+ const int theFlags,
+ const double theTolerance)
+{
+ ExtrusParam aParams;
+ aParams.myDir = gp_Dir(theStep);
+ aParams.myNodes.Clear();
+ aParams.mySteps = new TColStd_HSequenceOfReal;
+ int i;
+ for(i=1; i<=theNbSteps; i++)
+ aParams.mySteps->Append(theStep.Magnitude());
+
+ return
+ ExtrusionSweep(theElems,aParams,newElemsMap,theMakeGroups,theFlags,theTolerance);
+}
+
+
+//=======================================================================
+//function : ExtrusionSweep
+//purpose :
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems,
+ ExtrusParam& theParams,
+ TElemOfElemListMap& newElemsMap,
+ const bool theMakeGroups,
+ const int theFlags,
+ const double theTolerance)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ int nbsteps = theParams.mySteps->Length();
+
+ TNodeOfNodeListMap mapNewNodes;
+ //TNodeOfNodeVecMap mapNewNodes;
+ TElemOfVecOfNnlmiMap mapElemNewNodes;
+ //TElemOfVecOfMapNodesMap mapElemNewNodes;
+
+ const bool isQuadraticMesh = bool( myMesh->NbEdges(ORDER_QUADRATIC) +
+ myMesh->NbFaces(ORDER_QUADRATIC) +
+ myMesh->NbVolumes(ORDER_QUADRATIC) );
+ // loop on theElems
+ TIDSortedElemSet::iterator itElem;
+ for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
+ // check element type
+ 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() )
+ {
+ // check if a node has been already sweeped
+ const SMDS_MeshNode* node = cast2Node( 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;
+ }
+ }
+
+ double coord[] = { node->X(), node->Y(), node->Z() };
+ for ( int i = 0; i < nbsteps; i++ )
+ {
+ if ( needMediumNodes ) // create a medium node
+ {
+ double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1)/2.;
+ double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1)/2.;
+ double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1)/2.;
+ if( theFlags & EXTRUSION_FLAG_SEW ) {
+ const SMDS_MeshNode * newNode = CreateNode(x, y, z,
+ theTolerance, theParams.myNodes);
+ listNewNodes.push_back( newNode );
+ }
+ else {
+ const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
+ myLastCreatedNodes.Append(newNode);
+ srcNodes.Append( node );
+ listNewNodes.push_back( newNode );
+ }
+ }
+ // create a corner node
+ coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
+ coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
+ coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
+ if( theFlags & EXTRUSION_FLAG_SEW ) {
+ const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
+ theTolerance, theParams.myNodes);
+ listNewNodes.push_back( newNode );
+ }
+ else {
+ const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+ myLastCreatedNodes.Append(newNode);
+ srcNodes.Append( node );
+ listNewNodes.push_back( newNode );
+ }
+ }
+ }
+ newNodesItVec.push_back( nIt );
+ }
+ // make new elements
+ sweepElement( elem, newNodesItVec, newElemsMap[elem], nbsteps, srcElems );
+ }
+
+ if( theFlags & EXTRUSION_FLAG_BOUNDARY ) {
+ makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, nbsteps, srcElems );
+ }
+ PGroupIDs newGroupIDs;
+ if ( theMakeGroups )
+ newGroupIDs = generateGroups( srcNodes, srcElems, "extruded");
+
+ return newGroupIDs;
+}
+
+//=======================================================================
+//function : ExtrusionAlongTrack
+//purpose :
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements,
+ SMESH_subMesh* theTrack,
+ const SMDS_MeshNode* theN1,
+ const bool theHasAngles,
+ list<double>& theAngles,
+ const bool theLinearVariation,
+ const bool theHasRefPoint,
+ const gp_Pnt& theRefPoint,
+ const bool theMakeGroups)
+{
+ MESSAGE("ExtrusionAlongTrack");
+ 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.size();
+ // nothing to do
+ if ( !aNbE )
+ return EXTR_NO_ELEMENTS;
+
+ // 1.1 Track Pattern
+ ASSERT( theTrack );
+
+ SMESHDS_SubMesh* pSubMeshDS = theTrack->GetSubMeshDS();
+
+ aItE = pSubMeshDS->GetElements();
+ 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->GetSubShape();
+ // Sub-shape for the Pattern must be an Edge or Wire
+ if( aS.ShapeType() == TopAbs_EDGE ) {
+ aTrackEdge = TopoDS::Edge( aS );
+ // the Edge must not be degenerated
+ if ( BRep_Tool::Degenerated( aTrackEdge ) )
+ return EXTR_BAD_PATH_SHAPE;
+ TopExp::Vertices( aTrackEdge, aV1, aV2 );
+ aItN = theTrack->GetFather()->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN1 = aItN->next();
+ aItN = theTrack->GetFather()->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN2 = aItN->next();
+ // starting node must be aN1 or aN2
+ if ( !( aN1 == theN1 || aN2 == theN1 ) )
+ return EXTR_BAD_STARTING_NODE;
+ aItN = pSubMeshDS->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 );
+ }
+ //Extrusion_Error err =
+ MakeEdgePathPoints(aPrms, aTrackEdge, (aN1==theN1), fullList);
+ } else if( aS.ShapeType() == TopAbs_WIRE ) {
+ list< SMESH_subMesh* > LSM;
+ TopTools_SequenceOfShape Edges;
+ SMESH_subMeshIteratorPtr itSM = theTrack->getDependsOnIterator(false,true);
+ while(itSM->more()) {
+ SMESH_subMesh* SM = itSM->next();
+ LSM.push_back(SM);
+ const TopoDS_Shape& aS = SM->GetSubShape();
+ Edges.Append(aS);
+ }
+ list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
+ int startNid = theN1->GetID();
+ 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 );
+ aItN = locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN1 = aItN->next();
+ aItN = locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN2 = aItN->next();
+ // starting node must be aN1 or aN2
+ if ( !( aN1->GetID() == startNid || aN2->GetID() == startNid ) ) continue;
+ // 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,(aN1->GetID()==startNid), LPP);
+ LLPPs.push_back(LPP);
+ UsedNums.Add(k);
+ // update startN for search following egde
+ if( aN1->GetID() == startNid ) startNid = aN2->GetID();
+ else startNid = aN1->GetID();
+ break;
+ }
+ }
+ 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,
+ 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.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( aS == SMESH_Mesh::PseudoShape() ) {
+ //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 > 2)
+ 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 ( BRep_Tool::Degenerated( aTrackEdge ) )
+ return EXTR_BAD_PATH_SHAPE;
+ TopExp::Vertices( aTrackEdge, aV1, aV2 );
+ aItN = theTrack->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN1 = aItN->next();
+ aItN = theTrack->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN2 = aItN->next();
+ // 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( BRep_Tool::Degenerated(E) ) continue;
+ SMESH_subMesh* SM = theTrack->GetSubMesh(E);
+ if(SM) {
+ LSM.push_back(SM);
+ Edges.Append(E);
+ }
+ }
+ list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
+ int startNid = theN1->GetID();
+ 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 );
+ aItN = locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN1 = aItN->next();
+ aItN = locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN2 = aItN->next();
+ // starting node must be aN1 or aN2
+ if ( !( aN1->GetID() == startNid || aN2->GetID() == startNid ) ) continue;
+ // 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,(aN1->GetID()==startNid), LPP);
+ LLPPs.push_back(LPP);
+ UsedNums.Add(k);
+ // update startN for search following egde
+ if( aN1->GetID() == startNid ) startNid = aN2->GetID();
+ else startNid = aN1->GetID();
+ break;
+ }
+ }
+ 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!=currList.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 : 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& theElements,
+ list<SMESH_MeshEditor_PathPoint>& fullList,
+ const bool theHasAngles,
+ list<double>& theAngles,
+ const bool theLinearVariation,
+ const bool theHasRefPoint,
+ const gp_Pnt& theRefPoint,
+ const bool theMakeGroups)
+{
+ MESSAGE("MakeExtrElements");
+ //cout<<"MakeExtrElements fullList.size() = "<<fullList.size()<<endl;
+ int aNbTP = fullList.size();
+ vector<SMESH_MeshEditor_PathPoint> aPPs(aNbTP);
+ // Angles
+ if( theHasAngles && theAngles.size()>0 && theLinearVariation ) {
+ LinearAngleVariation(aNbTP-1, theAngles);
+ }
+ vector<double> aAngles( aNbTP );
+ int j = 0;
+ for(; j<aNbTP; ++j) {
+ aAngles[j] = 0.;
+ }
+ if ( theHasAngles ) {
+ double anAngle;;
+ std::list<double>::iterator aItD = theAngles.begin();
+ for ( j=1; (aItD != theAngles.end()) && (j<aNbTP); ++aItD, ++j ) {
+ anAngle = *aItD;
+ aAngles[j] = anAngle;
+ }
+ }
+ // fill vector of path points with angles
+ //aPPs.resize(fullList.size());
+ j = -1;
+ list<SMESH_MeshEditor_PathPoint>::iterator itPP = fullList.begin();
+ for(; itPP!=fullList.end(); itPP++) {
+ j++;
+ SMESH_MeshEditor_PathPoint PP = *itPP;
+ PP.SetAngle(aAngles[j]);
+ aPPs[j] = PP;
+ }
+
+ TNodeOfNodeListMap mapNewNodes;
+ TElemOfVecOfNnlmiMap mapElemNewNodes;
+ TElemOfElemListMap newElemsMap;
+ TIDSortedElemSet::iterator itElem;
+ double aX, aY, aZ;
+ int aNb;
+ SMDSAbs_ElementType aTypeE;
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+
+ // 3. Center of rotation aV0
+ gp_Pnt aV0 = theRefPoint;
+ gp_XYZ aGC;
+ if ( !theHasRefPoint ) {
+ aNb = 0;
+ aGC.SetCoord( 0.,0.,0. );
+
+ itElem = theElements.begin();
+ for ( ; itElem != theElements.end(); itElem++ ) {
+ const SMDS_MeshElement* elem = *itElem;
+
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() ) {
+ const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
+ aX = node->X();
+ aY = node->Y();
+ aZ = node->Z();
+
+ if ( mapNewNodes.find( node ) == mapNewNodes.end() ) {
+ list<const SMDS_MeshNode*> aLNx;
+ mapNewNodes[node] = aLNx;
+ //
+ gp_XYZ aXYZ( aX, aY, aZ );
+ aGC += aXYZ;
+ ++aNb;
+ }
+ }
+ }
+ aGC /= aNb;
+ aV0.SetXYZ( aGC );
+ } // if (!theHasRefPoint) {
+ mapNewNodes.clear();
+
+ // 4. Processing the elements
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) {
+ // check element type
+ const SMDS_MeshElement* elem = *itElem;
+ aTypeE = elem->GetType();
+ if ( !elem || ( 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
+ aX = node->X(); aY = node->Y(); aZ = node->Z();
+
+ 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.SetCoord(aX, aY, aZ);
+
+ 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 ( 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 );
+ }
+ aX = aPN1.X();
+ aY = aPN1.Y();
+ aZ = aPN1.Z();
+ const SMDS_MeshNode* newNode = aMesh->AddNode( aX, aY, aZ );
+ 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, theElements, 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";
+ }
+
+ SMESH_MeshEditor targetMeshEditor( theTargetMesh );
+ SMESHDS_Mesh* aTgtMesh = theTargetMesh ? theTargetMesh->GetMeshDS() : 0;
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+
+ // 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
+ 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;
+
+ double coord[3];
+ coord[0] = node->X();
+ coord[1] = node->Y();
+ coord[2] = node->Z();
+ 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 );
+ }
+ }
+ }
+ }
+
+ // either create new elements or reverse mirrored ones
+ if ( !theCopy && !needReverse && !theTargetMesh )
+ return PGroupIDs();
+
+ TIDSortedElemSet::iterator invElemIt = inverseElemSet.begin();
+ for ( ; invElemIt != inverseElemSet.end(); invElemIt++ )
+ theElems.insert( *invElemIt );
+
+ // Replicate or reverse elements
+
+ std::vector<int> iForw;
+ 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
+
+ switch ( geomType ) {
+
+ case SMDSGeom_POLYGON: // ---------------------- polygon
+ {
+ vector<const SMDS_MeshNode*> poly_nodes (nbNodes);
+ int iNode = 0;
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while (itN->more()) {
+ const SMDS_MeshNode* node =
+ static_cast<const SMDS_MeshNode*>(itN->next());
+ TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
+ if (nodeMapIt == nodeMap.end())
+ break; // not all nodes transformed
+ if (needReverse) {
+ // reverse mirrored faces and volumes
+ poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second;
+ } else {
+ poly_nodes[iNode] = (*nodeMapIt).second;
+ }
+ iNode++;
+ }
+ if ( iNode != nbNodes )
+ continue; // not all nodes transformed
+
+ if ( theTargetMesh ) {
+ myLastCreatedElems.Append(aTgtMesh->AddPolygonalFace(poly_nodes));
+ srcElems.Append( elem );
+ }
+ else if ( theCopy ) {
+ myLastCreatedElems.Append(aMesh->AddPolygonalFace(poly_nodes));
+ srcElems.Append( elem );
+ }
+ else {
+ aMesh->ChangePolygonNodes(elem, poly_nodes);
+ }
+ }
+ break;
+
+ case SMDSGeom_POLYHEDRA: // ------------------ polyhedral volume
+ {
+ const SMDS_VtkVolume* aPolyedre =
+ dynamic_cast<const SMDS_VtkVolume*>( elem );
+ if (!aPolyedre) {
+ MESSAGE("Warning: bad volumic element");
+ continue;
+ }
+
+ vector<const SMDS_MeshNode*> poly_nodes; poly_nodes.reserve( nbNodes );
+ vector<int> quantities; quantities.reserve( nbNodes );
+
+ bool allTransformed = true;
+ int nbFaces = aPolyedre->NbFaces();
+ for (int iface = 1; iface <= nbFaces && allTransformed; iface++) {
+ int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
+ for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) {
+ const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode);
+ TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
+ if (nodeMapIt == nodeMap.end()) {
+ allTransformed = false; // not all nodes transformed
+ } else {
+ poly_nodes.push_back((*nodeMapIt).second);
+ }
+ if ( needReverse && allTransformed )
+ std::reverse( poly_nodes.end() - nbFaceNodes, poly_nodes.end() );
+ }
+ quantities.push_back(nbFaceNodes);
+ }
+ if ( !allTransformed )
+ continue; // not all nodes transformed
+
+ if ( theTargetMesh ) {
+ myLastCreatedElems.Append(aTgtMesh->AddPolyhedralVolume(poly_nodes, quantities));
+ srcElems.Append( elem );
+ }
+ else if ( theCopy ) {
+ myLastCreatedElems.Append(aMesh->AddPolyhedralVolume(poly_nodes, quantities));
+ srcElems.Append( elem );
+ }
+ else {
+ aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
+ }
+ }
+ break;
+
+ case SMDSGeom_BALL: // -------------------- Ball
+ {
+ if ( !theCopy && !theTargetMesh ) continue;
+
+ TNodeNodeMap::iterator nodeMapIt = nodeMap.find( elem->GetNode(0) );
+ if (nodeMapIt == nodeMap.end())
+ continue; // not all nodes transformed
+
+ double diameter = static_cast<const SMDS_BallElement*>(elem)->GetDiameter();
+ if ( theTargetMesh ) {
+ myLastCreatedElems.Append(aTgtMesh->AddBall( nodeMapIt->second, diameter ));
+ srcElems.Append( elem );
+ }
+ else {
+ myLastCreatedElems.Append(aMesh->AddBall( nodeMapIt->second, diameter ));
+ srcElems.Append( elem );
+ }
+ }
+ break;
+
+ default: // ----------------------- Regular elements
+
+ while ( iForw.size() < nbNodes ) iForw.push_back( iForw.size() );
+ const std::vector<int>& iRev = SMDS_MeshCell::reverseSmdsOrder( elem->GetEntityType() );
+ const std::vector<int>& i = needReverse ? iRev : iForw;
+
+ // find transformed nodes
+ vector<const SMDS_MeshNode*> nodes(nbNodes);
+ int iNode = 0;
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() ) {
+ const SMDS_MeshNode* node =
+ static_cast<const SMDS_MeshNode*>( itN->next() );
+ TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
+ if ( nodeMapIt == nodeMap.end() )
+ break; // not all nodes transformed
+ nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
+ }
+ if ( iNode != nbNodes )
+ continue; // not all nodes transformed
+
+ if ( theTargetMesh ) {
+ if ( SMDS_MeshElement* copy =
+ targetMeshEditor.AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
+ myLastCreatedElems.Append( copy );
+ srcElems.Append( elem );
+ }
+ }
+ else if ( theCopy ) {
+ if ( AddElement( nodes, elem->GetType(), elem->IsPoly() ))
+ srcElems.Append( elem );
+ }
+ else {
+ // reverse element as it was reversed by transformation
+ if ( nbNodes > 2 )
+ aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes );
+ }
+ } // switch ( geomType )
+
+ } // loop on elements
+
+ PGroupIDs newGroupIDs;
+
+ if ( ( theMakeGroups && theCopy ) ||
+ ( theMakeGroups && theTargetMesh ) )
+ newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh );
+
+ 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
+ */
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
+ const SMESH_SequenceOfElemPtr& elemGens,
+ const std::string& postfix,
+ SMESH_Mesh* targetMesh)
+{
+ PGroupIDs newGroupIDs( new list<int> );
+ SMESH_Mesh* mesh = targetMesh ? targetMesh : GetMesh();
+
+ // Sort existing groups by types and collect their names
+
+ // to store an old group and a generated new one
+ typedef pair< SMESHDS_GroupBase*, 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() );
+ SMESHDS_Group* newGroup = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(),
+ groupDS->GetType() );
+ groupsByType[ groupDS->GetType() ].push_back( make_pair( groupDS, newGroup ));
+ orderedOldNewGroups.push_back( & groupsByType[ groupDS->GetType() ].back() );
+ }
+
+ // Loop on nodes and elements to add them in new groups
+
+ for ( int isNodes = 0; isNodes < 2; ++isNodes )
+ {
+ const SMESH_SequenceOfElemPtr& gens = isNodes ? nodeGens : elemGens;
+ const SMESH_SequenceOfElemPtr& elems = isNodes ? myLastCreatedNodes : myLastCreatedElems;
+ if ( gens.Length() != elems.Length() )
+ throw SALOME_Exception(LOCALIZED("invalid args"));
+
+ // 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 sourceElem
+ list< const SMDS_MeshElement* > resultElems;
+ 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 );
+
+ // add resultElems to groups made by ones the sourceElem belongs to
+ list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end();
+ for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew )
+ {
+ SMESHDS_GroupBase* oldGroup = gOldNew->first;
+ if ( oldGroup->Contains( sourceElem )) // sourceElem is in oldGroup
+ {
+ // fill in a new group
+ SMDS_MeshGroup & newGroup = gOldNew->second->SMDSGroup();
+ list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
+ for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt )
+ newGroup.Add( *resElemIt );
+ }
+ }
+ } // loop on created elements
+ }// loop on nodes and elements
+
+ // Create new SMESH_Groups from SMESHDS_Groups and remove empty SMESHDS_Groups
+
+ for ( size_t i = 0; i < orderedOldNewGroups.size(); ++i )
+ {
+ SMESHDS_GroupBase* oldGroupDS = orderedOldNewGroups[i]->first;
+ SMESHDS_Group* newGroupDS = orderedOldNewGroups[i]->second;
+ if ( newGroupDS->IsEmpty() )
+ {
+ mesh->GetMeshDS()->RemoveGroup( newGroupDS );
+ }
+ else
+ {
+ // make a name
+ string name = oldGroupDS->GetStoreName();
+ if ( !targetMesh ) {
+ name += "_";
+ name += postfix;
+ int nb = 1;
+ while ( !groupNames.insert( name ).second ) // name exists
+ name = SMESH_Comment( oldGroupDS->GetStoreName() ) << "_" << postfix << "_" << nb++;
+ }
+ newGroupDS->SetStoreName( name.c_str() );
+
+ // make a SMESH_Groups
+ mesh->AddGroup( newGroupDS );
+ newGroupIDs->push_back( newGroupDS->GetID() );
+
+ // set group type
+ newGroupDS->SetType( newGroupDS->GetElements()->next()->GetType() );
+ }
+ }
+
+ 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
+ */
+//================================================================================
+
+void SMESH_MeshEditor::FindCoincidentNodes (TIDSortedNodeSet & theNodes,
+ const double theTolerance,
+ TListOfListOfNodes & theGroupsOfNodes)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ if ( theNodes.empty() )
+ { // get all nodes in the mesh
+ SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator(/*idInceasingOrder=*/true);
+ while ( nIt->more() )
+ theNodes.insert( theNodes.end(),nIt->next());
+ }
+
+ SMESH_OctreeNode::FindCoincidentNodes ( theNodes, &theGroupsOfNodes, theTolerance);
+}
+
+
+//=======================================================================
+/*!
+ * \brief Implementation of search for the node closest to point
+ */
+//=======================================================================
+
+struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
+{
+ //---------------------------------------------------------------------
+ /*!
+ * \brief Constructor
+ */
+ SMESH_NodeSearcherImpl( const SMESHDS_Mesh* theMesh )
+ {
+ myMesh = ( SMESHDS_Mesh* ) theMesh;
+
+ TIDSortedNodeSet nodes;
+ if ( theMesh ) {
+ SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true);
+ while ( nIt->more() )
+ nodes.insert( nodes.end(), nIt->next() );
+ }
+ myOctreeNode = new SMESH_OctreeNode(nodes) ;
+
+ // get max size of a leaf box
+ SMESH_OctreeNode* tree = myOctreeNode;
+ while ( !tree->isLeaf() )
+ {
+ SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
+ if ( cIt->more() )
+ tree = cIt->next();
+ }
+ myHalfLeafSize = tree->maxSize() / 2.;
+ }
+
+ //---------------------------------------------------------------------
+ /*!
+ * \brief Move node and update myOctreeNode accordingly
+ */
+ void MoveNode( const SMDS_MeshNode* node, const gp_Pnt& toPnt )
+ {
+ myOctreeNode->UpdateByMoveNode( node, toPnt );
+ myMesh->MoveNode( node, toPnt.X(), toPnt.Y(), toPnt.Z() );
+ }
+
+ //---------------------------------------------------------------------
+ /*!
+ * \brief Do it's job
+ */
+ const SMDS_MeshNode* FindClosestTo( const gp_Pnt& thePnt )
+ {
+ map<double, const SMDS_MeshNode*> dist2Nodes;
+ myOctreeNode->NodesAround( thePnt.Coord(), dist2Nodes, myHalfLeafSize );
+ if ( !dist2Nodes.empty() )
+ return dist2Nodes.begin()->second;
+ list<const SMDS_MeshNode*> nodes;
+ //myOctreeNode->NodesAround( &tgtNode, &nodes, myHalfLeafSize );
+
+ double minSqDist = DBL_MAX;
+ if ( nodes.empty() ) // get all nodes of OctreeNode's closest to thePnt
+ {
+ // sort leafs by their distance from thePnt
+ typedef map< double, SMESH_OctreeNode* > TDistTreeMap;
+ TDistTreeMap treeMap;
+ list< SMESH_OctreeNode* > treeList;
+ list< SMESH_OctreeNode* >::iterator trIt;
+ treeList.push_back( myOctreeNode );
+
+ gp_XYZ pointNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
+ bool pointInside = myOctreeNode->isInside( pointNode, myHalfLeafSize );
+ for ( trIt = treeList.begin(); trIt != treeList.end(); ++trIt)
+ {
+ SMESH_OctreeNode* tree = *trIt;
+ if ( !tree->isLeaf() ) // put children to the queue
+ {
+ if ( pointInside && !tree->isInside( pointNode, myHalfLeafSize )) continue;
+ SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
+ while ( cIt->more() )
+ treeList.push_back( cIt->next() );
+ }
+ else if ( tree->NbNodes() ) // put a tree to the treeMap
+ {
+ const Bnd_B3d& box = tree->getBox();
+ double sqDist = thePnt.SquareDistance( 0.5 * ( box.CornerMin() + box.CornerMax() ));
+ pair<TDistTreeMap::iterator,bool> it_in = treeMap.insert( make_pair( sqDist, tree ));
+ if ( !it_in.second ) // not unique distance to box center
+ treeMap.insert( it_in.first, make_pair( sqDist + 1e-13*treeMap.size(), tree ));
+ }
+ }
+ // find distance after which there is no sense to check tree's
+ double sqLimit = DBL_MAX;
+ TDistTreeMap::iterator sqDist_tree = treeMap.begin();
+ if ( treeMap.size() > 5 ) {
+ SMESH_OctreeNode* closestTree = sqDist_tree->second;
+ const Bnd_B3d& box = closestTree->getBox();
+ double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() );
+ sqLimit = limit * limit;
+ }
+ // get all nodes from trees
+ for ( ; sqDist_tree != treeMap.end(); ++sqDist_tree) {
+ if ( sqDist_tree->first > sqLimit )
+ break;
+ SMESH_OctreeNode* tree = sqDist_tree->second;
+ tree->NodesAround( tree->GetNodeIterator()->next(), &nodes );
+ }
+ }
+ // find closest among nodes
+ minSqDist = DBL_MAX;
+ const SMDS_MeshNode* closestNode = 0;
+ list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
+ for ( ; nIt != nodes.end(); ++nIt ) {
+ double sqDist = thePnt.SquareDistance( SMESH_TNodeXYZ( *nIt ) );
+ if ( minSqDist > sqDist ) {
+ closestNode = *nIt;
+ minSqDist = sqDist;
+ }
+ }
+ return closestNode;
+ }
+
+ //---------------------------------------------------------------------
+ /*!
+ * \brief Destructor
+ */
+ ~SMESH_NodeSearcherImpl() { delete myOctreeNode; }
+
+ //---------------------------------------------------------------------
+ /*!
+ * \brief Return the node tree
+ */
+ const SMESH_OctreeNode* getTree() const { return myOctreeNode; }
+
+private:
+ SMESH_OctreeNode* myOctreeNode;
+ SMESHDS_Mesh* myMesh;
+ double myHalfLeafSize; // max size of a leaf box
+};
+
+//=======================================================================
+/*!
+ * \brief Return SMESH_NodeSearcher
+ */
+//=======================================================================
+
+SMESH_NodeSearcher* SMESH_MeshEditor::GetNodeSearcher()
+{
+ return new SMESH_NodeSearcherImpl( GetMeshDS() );
+}
+
+// ========================================================================
+namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
+{
+ const int MaxNbElemsInLeaf = 10; // maximal number of elements in a leaf of tree
+ const int MaxLevel = 7; // maximal tree height -> nb terminal boxes: 8^7 = 2097152
+ const double NodeRadius = 1e-9; // to enlarge bnd box of element
+
+ //=======================================================================
+ /*!
+ * \brief Octal tree of bounding boxes of elements
+ */
+ //=======================================================================
+
+ class ElementBndBoxTree : public SMESH_Octree
+ {
+ public:
+
+ ElementBndBoxTree(const SMDS_Mesh& mesh,
+ SMDSAbs_ElementType elemType,
+ SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(),
+ double tolerance = NodeRadius );
+ void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems );
+ void getElementsNearLine ( const gp_Ax1& line, TIDSortedElemSet& foundElems);
+ void getElementsInSphere ( const gp_XYZ& center,
+ const double radius, TIDSortedElemSet& foundElems);
+ size_t getSize() { return std::max( _size, _elements.size() ); }
+ ~ElementBndBoxTree();
+
+ protected:
+ ElementBndBoxTree():_size(0) {}
+ SMESH_Octree* allocateOctreeChild() const { return new ElementBndBoxTree; }
+ void buildChildrenData();
+ Bnd_B3d* buildRootBox();
+ private:
+ //!< Bounding box of element
+ struct ElementBox : public Bnd_B3d
+ {
+ const SMDS_MeshElement* _element;
+ int _refCount; // an ElementBox can be included in several tree branches
+ ElementBox(const SMDS_MeshElement* elem, double tolerance);
+ };
+ vector< ElementBox* > _elements;
+ size_t _size;
+ };
+
+ //================================================================================
+ /*!
+ * \brief ElementBndBoxTree creation
+ */
+ //================================================================================
+
+ ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt, double tolerance)
+ :SMESH_Octree( new SMESH_Octree::Limit( MaxLevel, /*minSize=*/0. ))
+ {
+ int nbElems = mesh.GetMeshInfo().NbElements( elemType );
+ _elements.reserve( nbElems );
+
+ SMDS_ElemIteratorPtr elemIt = theElemIt ? theElemIt : mesh.elementsIterator( elemType );
+ while ( elemIt->more() )
+ _elements.push_back( new ElementBox( elemIt->next(),tolerance ));
+
+ compute();
+ }
+
+ //================================================================================
+ /*!
+ * \brief Destructor
+ */
+ //================================================================================
+
+ ElementBndBoxTree::~ElementBndBoxTree()
+ {
+ for ( int i = 0; i < _elements.size(); ++i )
+ if ( --_elements[i]->_refCount <= 0 )
+ delete _elements[i];
+ }
+
+ //================================================================================
+ /*!
+ * \brief Return the maximal box
+ */
+ //================================================================================
+
+ Bnd_B3d* ElementBndBoxTree::buildRootBox()
+ {
+ Bnd_B3d* box = new Bnd_B3d;
+ for ( int i = 0; i < _elements.size(); ++i )
+ box->Add( *_elements[i] );
+ return box;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Redistrubute element boxes among children
+ */
+ //================================================================================
+
+ void ElementBndBoxTree::buildChildrenData()
+ {
+ for ( int i = 0; i < _elements.size(); ++i )
+ {
+ for (int j = 0; j < 8; j++)
+ {
+ if ( !_elements[i]->IsOut( myChildren[j]->getBox() ))
+ {
+ _elements[i]->_refCount++;
+ ((ElementBndBoxTree*)myChildren[j])->_elements.push_back( _elements[i]);
+ }
+ }
+ _elements[i]->_refCount--;
+ }
+ _size = _elements.size();
+ SMESHUtils::FreeVector( _elements ); // = _elements.clear() + free memory
+
+ for (int j = 0; j < 8; j++)
+ {
+ ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j]);
+ if ( child->_elements.size() <= MaxNbElemsInLeaf )
+ child->myIsLeaf = true;
+
+ if ( child->_elements.capacity() - child->_elements.size() > 1000 )
+ SMESHUtils::CompactVector( child->_elements );
+ }
+ }
+
+ //================================================================================
+ /*!
+ * \brief Return elements which can include the point
+ */
+ //================================================================================
+
+ void ElementBndBoxTree::getElementsNearPoint( const gp_Pnt& point,
+ TIDSortedElemSet& foundElems)
+ {
+ if ( getBox().IsOut( point.XYZ() ))
+ return;
+
+ if ( isLeaf() )
+ {
+ for ( int i = 0; i < _elements.size(); ++i )
+ if ( !_elements[i]->IsOut( point.XYZ() ))
+ foundElems.insert( _elements[i]->_element );
+ }
+ else
+ {
+ for (int i = 0; i < 8; i++)
+ ((ElementBndBoxTree*) myChildren[i])->getElementsNearPoint( point, foundElems );
+ }
+ }
+
+ //================================================================================
+ /*!
+ * \brief Return elements which can be intersected by the line
+ */
+ //================================================================================
+
+ void ElementBndBoxTree::getElementsNearLine( const gp_Ax1& line,
+ TIDSortedElemSet& foundElems)
+ {
+ if ( getBox().IsOut( line ))
+ return;
+
+ if ( isLeaf() )
+ {
+ for ( int i = 0; i < _elements.size(); ++i )
+ if ( !_elements[i]->IsOut( line ))
+ foundElems.insert( _elements[i]->_element );
+ }
+ else
+ {
+ for (int i = 0; i < 8; i++)
+ ((ElementBndBoxTree*) myChildren[i])->getElementsNearLine( line, foundElems );
+ }
+ }
+
+ //================================================================================
+ /*!
+ * \brief Return elements from leaves intersecting the sphere
+ */
+ //================================================================================
+
+ void ElementBndBoxTree::getElementsInSphere ( const gp_XYZ& center,
+ const double radius,
+ TIDSortedElemSet& foundElems)
+ {
+ if ( getBox().IsOut( center, radius ))
+ return;
+
+ if ( isLeaf() )
+ {
+ for ( int i = 0; i < _elements.size(); ++i )
+ if ( !_elements[i]->IsOut( center, radius ))
+ foundElems.insert( _elements[i]->_element );
+ }
+ else
+ {
+ for (int i = 0; i < 8; i++)
+ ((ElementBndBoxTree*) myChildren[i])->getElementsInSphere( center, radius, foundElems );
+ }
+ }
+
+ //================================================================================
+ /*!
+ * \brief Construct the element box
+ */
+ //================================================================================
+
+ ElementBndBoxTree::ElementBox::ElementBox(const SMDS_MeshElement* elem, double tolerance)
+ {
+ _element = elem;
+ _refCount = 1;
+ SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
+ while ( nIt->more() )
+ Add( SMESH_TNodeXYZ( nIt->next() ));
+ Enlarge( tolerance );
+ }
+
+} // namespace
+
+//=======================================================================
+/*!
+ * \brief Implementation of search for the elements by point and
+ * of classification of point in 2D mesh
+ */
+//=======================================================================
+
+struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
+{
+ SMESHDS_Mesh* _mesh;
+ SMDS_ElemIteratorPtr _meshPartIt;
+ ElementBndBoxTree* _ebbTree;
+ SMESH_NodeSearcherImpl* _nodeSearcher;
+ SMDSAbs_ElementType _elementType;
+ double _tolerance;
+ bool _outerFacesFound;
+ set<const SMDS_MeshElement*> _outerFaces; // empty means "no internal faces at all"
+
+ SMESH_ElementSearcherImpl( SMESHDS_Mesh& mesh, SMDS_ElemIteratorPtr elemIt=SMDS_ElemIteratorPtr())
+ : _mesh(&mesh),_meshPartIt(elemIt),_ebbTree(0),_nodeSearcher(0),_tolerance(-1),_outerFacesFound(false) {}
+ ~SMESH_ElementSearcherImpl()
+ {
+ if ( _ebbTree ) delete _ebbTree; _ebbTree = 0;
+ if ( _nodeSearcher ) delete _nodeSearcher; _nodeSearcher = 0;
+ }
+ virtual int FindElementsByPoint(const gp_Pnt& point,
+ SMDSAbs_ElementType type,
+ vector< const SMDS_MeshElement* >& foundElements);
+ virtual TopAbs_State GetPointState(const gp_Pnt& point);
+ virtual const SMDS_MeshElement* FindClosestTo( const gp_Pnt& point,
+ SMDSAbs_ElementType type );
+
+ void GetElementsNearLine( const gp_Ax1& line,
+ SMDSAbs_ElementType type,
+ vector< const SMDS_MeshElement* >& foundElems);
+ double getTolerance();
+ bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face,
+ const double tolerance, double & param);
+ void findOuterBoundary(const SMDS_MeshElement* anyOuterFace);
+ bool isOuterBoundary(const SMDS_MeshElement* face) const
+ {
+ return _outerFaces.empty() || _outerFaces.count(face);
+ }
+ struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState())
+ {
+ const SMDS_MeshElement* _face;
+ gp_Vec _faceNorm;
+ bool _coincides; //!< the line lays in face plane
+ TInters(const SMDS_MeshElement* face, const gp_Vec& faceNorm, bool coinc=false)
+ : _face(face), _faceNorm( faceNorm ), _coincides( coinc ) {}
+ };
+ struct TFaceLink //!< link and faces sharing it (used in findOuterBoundary())
+ {
+ SMESH_TLink _link;
+ TIDSortedElemSet _faces;
+ TFaceLink( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshElement* face)
+ : _link( n1, n2 ), _faces( &face, &face + 1) {}
+ };
+};
+
+ostream& operator<< (ostream& out, const SMESH_ElementSearcherImpl::TInters& i)
+{
+ return out << "TInters(face=" << ( i._face ? i._face->GetID() : 0)
+ << ", _coincides="<<i._coincides << ")";
+}
+
+//=======================================================================
+/*!
+ * \brief define tolerance for search
+ */
+//=======================================================================
+
+double SMESH_ElementSearcherImpl::getTolerance()
+{
+ if ( _tolerance < 0 )
+ {
+ const SMDS_MeshInfo& meshInfo = _mesh->GetMeshInfo();
+
+ _tolerance = 0;
+ if ( _nodeSearcher && meshInfo.NbNodes() > 1 )
+ {
+ double boxSize = _nodeSearcher->getTree()->maxSize();
+ _tolerance = 1e-8 * boxSize/* / meshInfo.NbNodes()*/;
+ }
+ else if ( _ebbTree && meshInfo.NbElements() > 0 )
+ {
+ double boxSize = _ebbTree->maxSize();
+ _tolerance = 1e-8 * boxSize/* / meshInfo.NbElements()*/;
+ }
+ if ( _tolerance == 0 )
+ {
+ // define tolerance by size of a most complex element
+ int complexType = SMDSAbs_Volume;
+ while ( complexType > SMDSAbs_All &&
+ meshInfo.NbElements( SMDSAbs_ElementType( complexType )) < 1 )
+ --complexType;
+ if ( complexType == SMDSAbs_All ) return 0; // empty mesh
+ double elemSize;
+ if ( complexType == int( SMDSAbs_Node ))
+ {
+ SMDS_NodeIteratorPtr nodeIt = _mesh->nodesIterator();
+ elemSize = 1;
+ if ( meshInfo.NbNodes() > 2 )
+ elemSize = SMESH_TNodeXYZ( nodeIt->next() ).Distance( nodeIt->next() );
+ }
+ else
+ {
+ SMDS_ElemIteratorPtr elemIt =
+ _mesh->elementsIterator( SMDSAbs_ElementType( complexType ));
+ const SMDS_MeshElement* elem = elemIt->next();
+ SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+ SMESH_TNodeXYZ n1( cast2Node( nodeIt->next() ));
+ elemSize = 0;
+ while ( nodeIt->more() )
+ {
+ double dist = n1.Distance( cast2Node( nodeIt->next() ));
+ elemSize = max( dist, elemSize );
+ }
+ }
+ _tolerance = 1e-4 * elemSize;
+ }
+ }
+ return _tolerance;
+}
+
+//================================================================================
+/*!
+ * \brief Find intersection of the line and an edge of face and return parameter on line
+ */
+//================================================================================
+
+bool SMESH_ElementSearcherImpl::getIntersParamOnLine(const gp_Lin& line,
+ const SMDS_MeshElement* face,
+ const double tol,
+ double & param)
+{
+ int nbInts = 0;
+ param = 0;
+
+ GeomAPI_ExtremaCurveCurve anExtCC;
+ Handle(Geom_Curve) lineCurve = new Geom_Line( line );
+
+ int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes();
+ for ( int i = 0; i < nbNodes && nbInts < 2; ++i )
+ {
+ GC_MakeSegment edge( SMESH_TNodeXYZ( face->GetNode( i )),
+ SMESH_TNodeXYZ( face->GetNode( (i+1)%nbNodes) ));
+ anExtCC.Init( lineCurve, edge);
+ if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol)
+ {
+ Quantity_Parameter pl, pe;
+ anExtCC.LowerDistanceParameters( pl, pe );
+ param += pl;
+ if ( ++nbInts == 2 )
+ break;
+ }
+ }
+ if ( nbInts > 0 ) param /= nbInts;
+ return nbInts > 0;
+}
+//================================================================================
+/*!
+ * \brief Find all faces belonging to the outer boundary of mesh
+ */
+//================================================================================
+
+void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerFace)
+{
+ if ( _outerFacesFound ) return;
+
+ // Collect all outer faces by passing from one outer face to another via their links
+ // and BTW find out if there are internal faces at all.
+
+ // checked links and links where outer boundary meets internal one
+ set< SMESH_TLink > visitedLinks, seamLinks;
+
+ // links to treat with already visited faces sharing them
+ list < TFaceLink > startLinks;
+
+ // load startLinks with the first outerFace
+ startLinks.push_back( TFaceLink( outerFace->GetNode(0), outerFace->GetNode(1), outerFace));
+ _outerFaces.insert( outerFace );
+
+ TIDSortedElemSet emptySet;
+ while ( !startLinks.empty() )
+ {
+ const SMESH_TLink& link = startLinks.front()._link;
+ TIDSortedElemSet& faces = startLinks.front()._faces;
+
+ outerFace = *faces.begin();
+ // find other faces sharing the link
+ const SMDS_MeshElement* f;
+ while (( f = SMESH_MeshEditor::FindFaceInSet(link.node1(), link.node2(), emptySet, faces )))
+ faces.insert( f );
+
+ // select another outer face among the found
+ const SMDS_MeshElement* outerFace2 = 0;
+ if ( faces.size() == 2 )
+ {
+ outerFace2 = (outerFace == *faces.begin() ? *faces.rbegin() : *faces.begin());
+ }
+ else if ( faces.size() > 2 )
+ {
+ seamLinks.insert( link );
+
+ // link direction within the outerFace
+ gp_Vec n1n2( SMESH_TNodeXYZ( link.node1()),
+ SMESH_TNodeXYZ( link.node2()));
+ int i1 = outerFace->GetNodeIndex( link.node1() );
+ int i2 = outerFace->GetNodeIndex( link.node2() );
+ bool rev = ( abs(i2-i1) == 1 ? i1 > i2 : i2 > i1 );
+ if ( rev ) n1n2.Reverse();
+ // outerFace normal
+ gp_XYZ ofNorm, fNorm;
+ if ( SMESH_Algo::FaceNormal( outerFace, ofNorm, /*normalized=*/false ))
+ {
+ // direction from the link inside outerFace
+ gp_Vec dirInOF = gp_Vec( ofNorm ) ^ n1n2;
+ // sort all other faces by angle with the dirInOF
+ map< double, const SMDS_MeshElement* > angle2Face;
+ set< const SMDS_MeshElement*, TIDCompare >::const_iterator face = faces.begin();
+ for ( ; face != faces.end(); ++face )
+ {
+ if ( !SMESH_Algo::FaceNormal( *face, fNorm, /*normalized=*/false ))
+ continue;
+ gp_Vec dirInF = gp_Vec( fNorm ) ^ n1n2;
+ double angle = dirInOF.AngleWithRef( dirInF, n1n2 );
+ if ( angle < 0 ) angle += 2. * M_PI;
+ angle2Face.insert( make_pair( angle, *face ));
+ }
+ if ( !angle2Face.empty() )
+ outerFace2 = angle2Face.begin()->second;
+ }
+ }
+ // store the found outer face and add its links to continue seaching from
+ if ( outerFace2 )
+ {
+ _outerFaces.insert( outerFace );
+ int nbNodes = outerFace2->NbNodes()/( outerFace2->IsQuadratic() ? 2 : 1 );
+ for ( int i = 0; i < nbNodes; ++i )
+ {
+ SMESH_TLink link2( outerFace2->GetNode(i), outerFace2->GetNode((i+1)%nbNodes));
+ if ( visitedLinks.insert( link2 ).second )
+ startLinks.push_back( TFaceLink( link2.node1(), link2.node2(), outerFace2 ));
+ }
+ }
+ startLinks.pop_front();
+ }
+ _outerFacesFound = true;
+
+ if ( !seamLinks.empty() )
+ {
+ // There are internal boundaries touching the outher one,
+ // find all faces of internal boundaries in order to find
+ // faces of boundaries of holes, if any.
+
+ }
+ else
+ {
+ _outerFaces.clear();
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Find elements of given type where the given point is IN or ON.
+ * Returns nb of found elements and elements them-selves.
+ *
+ * 'ALL' type means elements of any type excluding nodes, balls and 0D elements
+ */
+//=======================================================================
+
+int SMESH_ElementSearcherImpl::
+FindElementsByPoint(const gp_Pnt& point,
+ SMDSAbs_ElementType type,
+ vector< const SMDS_MeshElement* >& foundElements)
+{
+ foundElements.clear();
+
+ double tolerance = getTolerance();
+
+ // =================================================================================
+ if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement || type == SMDSAbs_Ball)
+ {
+ if ( !_nodeSearcher )
+ _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
+
+ const SMDS_MeshNode* closeNode = _nodeSearcher->FindClosestTo( point );
+ if ( !closeNode ) return foundElements.size();
+
+ if ( point.Distance( SMESH_TNodeXYZ( closeNode )) > tolerance )
+ return foundElements.size(); // to far from any node
+
+ if ( type == SMDSAbs_Node )
+ {
+ foundElements.push_back( closeNode );
+ }
+ else
+ {
+ SMDS_ElemIteratorPtr elemIt = closeNode->GetInverseElementIterator( type );
+ while ( elemIt->more() )
+ foundElements.push_back( elemIt->next() );
+ }
+ }
+ // =================================================================================
+ else // elements more complex than 0D
+ {
+ if ( !_ebbTree || _elementType != type )
+ {
+ if ( _ebbTree ) delete _ebbTree;
+ _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt, tolerance );
+ }
+ TIDSortedElemSet suspectElems;
+ _ebbTree->getElementsNearPoint( point, suspectElems );
+ TIDSortedElemSet::iterator elem = suspectElems.begin();
+ for ( ; elem != suspectElems.end(); ++elem )
+ if ( !SMESH_MeshEditor::IsOut( *elem, point, tolerance ))
+ foundElements.push_back( *elem );
+ }
+ return foundElements.size();
+}
+
+//=======================================================================
+/*!
+ * \brief Find an element of given type most close to the given point
+ *
+ * WARNING: Only face search is implemeneted so far
+ */
+//=======================================================================
+
+const SMDS_MeshElement*
+SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt& point,
+ SMDSAbs_ElementType type )
+{
+ const SMDS_MeshElement* closestElem = 0;
+
+ if ( type == SMDSAbs_Face )
+ {
+ if ( !_ebbTree || _elementType != type )
+ {
+ if ( _ebbTree ) delete _ebbTree;
+ _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
+ }
+ TIDSortedElemSet suspectElems;
+ _ebbTree->getElementsNearPoint( point, suspectElems );
+
+ if ( suspectElems.empty() && _ebbTree->maxSize() > 0 )
+ {
+ gp_Pnt boxCenter = 0.5 * ( _ebbTree->getBox().CornerMin() +
+ _ebbTree->getBox().CornerMax() );
+ double radius;
+ if ( _ebbTree->getBox().IsOut( point.XYZ() ))
+ radius = point.Distance( boxCenter ) - 0.5 * _ebbTree->maxSize();
+ else
+ radius = _ebbTree->maxSize() / pow( 2., _ebbTree->getHeight()) / 2;
+ while ( suspectElems.empty() )
+ {
+ _ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems );
+ radius *= 1.1;
+ }
+ }
+ double minDist = std::numeric_limits<double>::max();
+ multimap< double, const SMDS_MeshElement* > dist2face;
+ TIDSortedElemSet::iterator elem = suspectElems.begin();
+ for ( ; elem != suspectElems.end(); ++elem )
+ {
+ double dist = SMESH_MeshEditor::GetDistance( dynamic_cast<const SMDS_MeshFace*>(*elem),
+ point );
+ if ( dist < minDist + 1e-10)
+ {
+ minDist = dist;
+ dist2face.insert( dist2face.begin(), make_pair( dist, *elem ));
+ }
+ }
+ if ( !dist2face.empty() )
+ {
+ multimap< double, const SMDS_MeshElement* >::iterator d2f = dist2face.begin();
+ closestElem = d2f->second;
+ // if there are several elements at the same distance, select one
+ // with GC closest to the point
+ typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
+ double minDistToGC = 0;
+ for ( ++d2f; d2f != dist2face.end() && fabs( d2f->first - minDist ) < 1e-10; ++d2f )
+ {
+ if ( minDistToGC == 0 )
+ {
+ gp_XYZ gc(0,0,0);
+ gc = accumulate( TXyzIterator(closestElem->nodesIterator()),
+ TXyzIterator(), gc ) / closestElem->NbNodes();
+ minDistToGC = point.SquareDistance( gc );
+ }
+ gp_XYZ gc(0,0,0);
+ gc = accumulate( TXyzIterator( d2f->second->nodesIterator()),
+ TXyzIterator(), gc ) / d2f->second->NbNodes();
+ double d = point.SquareDistance( gc );
+ if ( d < minDistToGC )
+ {
+ minDistToGC = d;
+ closestElem = d2f->second;
+ }
+ }
+ // cout << "FindClosestTo( " <<point.X()<<", "<<point.Y()<<", "<<point.Z()<<" ) FACE "
+ // <<closestElem->GetID() << " DIST " << minDist << endl;
+ }
+ }
+ else
+ {
+ // NOT IMPLEMENTED SO FAR
+ }
+ return closestElem;
+}
+
+
+//================================================================================
+/*!
+ * \brief Classify the given point in the closed 2D mesh
+ */
+//================================================================================
+
+TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
+{
+ double tolerance = getTolerance();
+ if ( !_ebbTree || _elementType != SMDSAbs_Face )
+ {
+ if ( _ebbTree ) delete _ebbTree;
+ _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = SMDSAbs_Face, _meshPartIt );
+ }
+ // Algo: analyse transition of a line starting at the point through mesh boundary;
+ // try three lines parallel to axis of the coordinate system and perform rough
+ // analysis. If solution is not clear perform thorough analysis.
+
+ const int nbAxes = 3;
+ gp_Dir axisDir[ nbAxes ] = { gp::DX(), gp::DY(), gp::DZ() };
+ map< double, TInters > paramOnLine2TInters[ nbAxes ];
+ list< TInters > tangentInters[ nbAxes ]; // of faces whose plane includes the line
+ multimap< int, int > nbInt2Axis; // to find the simplest case
+ for ( int axis = 0; axis < nbAxes; ++axis )
+ {
+ gp_Ax1 lineAxis( point, axisDir[axis]);
+ gp_Lin line ( lineAxis );
+
+ TIDSortedElemSet suspectFaces; // faces possibly intersecting the line
+ _ebbTree->getElementsNearLine( lineAxis, suspectFaces );
+
+ // Intersect faces with the line
+
+ map< double, TInters > & u2inters = paramOnLine2TInters[ axis ];
+ TIDSortedElemSet::iterator face = suspectFaces.begin();
+ for ( ; face != suspectFaces.end(); ++face )
+ {
+ // get face plane
+ gp_XYZ fNorm;
+ if ( !SMESH_Algo::FaceNormal( *face, fNorm, /*normalized=*/false)) continue;
+ gp_Pln facePlane( SMESH_TNodeXYZ( (*face)->GetNode(0)), fNorm );
+
+ // perform intersection
+ IntAna_IntConicQuad intersection( line, IntAna_Quadric( facePlane ));
+ if ( !intersection.IsDone() )
+ continue;
+ if ( intersection.IsInQuadric() )
+ {
+ tangentInters[ axis ].push_back( TInters( *face, fNorm, true ));
+ }
+ else if ( ! intersection.IsParallel() && intersection.NbPoints() > 0 )
+ {
+ gp_Pnt intersectionPoint = intersection.Point(1);
+ if ( !SMESH_MeshEditor::IsOut( *face, intersectionPoint, tolerance ))
+ u2inters.insert(make_pair( intersection.ParamOnConic(1), TInters( *face, fNorm )));
+ }
+ }
+ // Analyse intersections roughly
+
+ int nbInter = u2inters.size();
+ if ( nbInter == 0 )
+ return TopAbs_OUT;
+
+ double f = u2inters.begin()->first, l = u2inters.rbegin()->first;
+ if ( nbInter == 1 ) // not closed mesh
+ return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN;
+
+ if ( fabs( f ) < tolerance || fabs( l ) < tolerance )
+ return TopAbs_ON;
+
+ if ( (f<0) == (l<0) )
+ return TopAbs_OUT;
+
+ int nbIntBeforePoint = std::distance( u2inters.begin(), u2inters.lower_bound(0));
+ int nbIntAfterPoint = nbInter - nbIntBeforePoint;
+ if ( nbIntBeforePoint == 1 || nbIntAfterPoint == 1 )
+ return TopAbs_IN;
+
+ nbInt2Axis.insert( make_pair( min( nbIntBeforePoint, nbIntAfterPoint ), axis ));
+
+ if ( _outerFacesFound ) break; // pass to thorough analysis
+
+ } // three attempts - loop on CS axes
+
+ // Analyse intersections thoroughly.
+ // We make two loops maximum, on the first one we only exclude touching intersections,
+ // on the second, if situation is still unclear, we gather and use information on
+ // position of faces (internal or outer). If faces position is already gathered,
+ // we make the second loop right away.
+
+ for ( int hasPositionInfo = _outerFacesFound; hasPositionInfo < 2; ++hasPositionInfo )
+ {
+ multimap< int, int >::const_iterator nb_axis = nbInt2Axis.begin();
+ for ( ; nb_axis != nbInt2Axis.end(); ++nb_axis )
+ {
+ int axis = nb_axis->second;
+ map< double, TInters > & u2inters = paramOnLine2TInters[ axis ];
+
+ gp_Ax1 lineAxis( point, axisDir[axis]);
+ gp_Lin line ( lineAxis );
+
+ // add tangent intersections to u2inters
+ double param;
+ list< TInters >::const_iterator tgtInt = tangentInters[ axis ].begin();
+ for ( ; tgtInt != tangentInters[ axis ].end(); ++tgtInt )
+ if ( getIntersParamOnLine( line, tgtInt->_face, tolerance, param ))
+ u2inters.insert(make_pair( param, *tgtInt ));
+ tangentInters[ axis ].clear();
+
+ // Count intersections before and after the point excluding touching ones.
+ // If hasPositionInfo we count intersections of outer boundary only
+
+ int nbIntBeforePoint = 0, nbIntAfterPoint = 0;
+ double f = numeric_limits<double>::max(), l = -numeric_limits<double>::max();
+ map< double, TInters >::iterator u_int1 = u2inters.begin(), u_int2 = u_int1;
+ bool ok = ! u_int1->second._coincides;
+ while ( ok && u_int1 != u2inters.end() )
+ {
+ double u = u_int1->first;
+ bool touchingInt = false;
+ if ( ++u_int2 != u2inters.end() )
+ {
+ // skip intersections at the same point (if the line passes through edge or node)
+ int nbSamePnt = 0;
+ while ( u_int2 != u2inters.end() && fabs( u_int2->first - u ) < tolerance )
+ {
+ ++nbSamePnt;
+ ++u_int2;
+ }
+
+ // skip tangent intersections
+ int nbTgt = 0;
+ const SMDS_MeshElement* prevFace = u_int1->second._face;
+ while ( ok && u_int2->second._coincides )
+ {
+ if ( SMESH_Algo::GetCommonNodes(prevFace , u_int2->second._face).empty() )
+ ok = false;
+ else
+ {
+ nbTgt++;
+ u_int2++;
+ ok = ( u_int2 != u2inters.end() );
+ }
+ }
+ if ( !ok ) break;
+
+ // skip intersections at the same point after tangent intersections
+ if ( nbTgt > 0 )
+ {
+ double u2 = u_int2->first;
+ ++u_int2;
+ while ( u_int2 != u2inters.end() && fabs( u_int2->first - u2 ) < tolerance )
+ {
+ ++nbSamePnt;
+ ++u_int2;
+ }
+ }
+ // decide if we skipped a touching intersection
+ if ( nbSamePnt + nbTgt > 0 )
+ {
+ double minDot = numeric_limits<double>::max(), maxDot = -numeric_limits<double>::max();
+ map< double, TInters >::iterator u_int = u_int1;
+ for ( ; u_int != u_int2; ++u_int )
+ {
+ if ( u_int->second._coincides ) continue;
+ double dot = u_int->second._faceNorm * line.Direction();
+ if ( dot > maxDot ) maxDot = dot;
+ if ( dot < minDot ) minDot = dot;
+ }
+ touchingInt = ( minDot*maxDot < 0 );
+ }
+ }
+ if ( !touchingInt )
+ {
+ if ( !hasPositionInfo || isOuterBoundary( u_int1->second._face ))
+ {
+ if ( u < 0 )
+ ++nbIntBeforePoint;
+ else
+ ++nbIntAfterPoint;
+ }
+ if ( u < f ) f = u;
+ if ( u > l ) l = u;
+ }
+
+ u_int1 = u_int2; // to next intersection
+
+ } // loop on intersections with one line
+
+ if ( ok )
+ {
+ if ( fabs( f ) < tolerance || fabs( l ) < tolerance )
+ return TopAbs_ON;
+
+ if ( nbIntBeforePoint == 0 || nbIntAfterPoint == 0)
+ return TopAbs_OUT;
+
+ if ( nbIntBeforePoint + nbIntAfterPoint == 1 ) // not closed mesh
+ return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN;
+
+ if ( nbIntBeforePoint == 1 || nbIntAfterPoint == 1 )
+ return TopAbs_IN;
+
+ if ( (f<0) == (l<0) )
+ return TopAbs_OUT;
+
+ if ( hasPositionInfo )
+ return nbIntBeforePoint % 2 ? TopAbs_IN : TopAbs_OUT;
+ }
+ } // loop on intersections of the tree lines - thorough analysis
+
+ if ( !hasPositionInfo )
+ {
+ // gather info on faces position - is face in the outer boundary or not
+ map< double, TInters > & u2inters = paramOnLine2TInters[ 0 ];
+ findOuterBoundary( u2inters.begin()->second._face );
+ }
+
+ } // two attempts - with and w/o faces position info in the mesh
+
+ return TopAbs_UNKNOWN;
+}
+
+//=======================================================================
+/*!
+ * \brief Return elements possibly intersecting the line
+ */
+//=======================================================================
+
+void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1& line,
+ SMDSAbs_ElementType type,
+ vector< const SMDS_MeshElement* >& foundElems)
+{
+ if ( !_ebbTree || _elementType != type )
+ {
+ if ( _ebbTree ) delete _ebbTree;
+ _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
+ }
+ TIDSortedElemSet suspectFaces; // elements possibly intersecting the line
+ _ebbTree->getElementsNearLine( line, suspectFaces );
+ foundElems.assign( suspectFaces.begin(), suspectFaces.end());
+}
+
+//=======================================================================
+/*!
+ * \brief Return SMESH_ElementSearcher
+ */
+//=======================================================================
+
+SMESH_ElementSearcher* SMESH_MeshEditor::GetElementSearcher()
+{
+ return new SMESH_ElementSearcherImpl( *GetMeshDS() );
+}
+
+//=======================================================================
+/*!
+ * \brief Return SMESH_ElementSearcher acting on a sub-set of elements
+ */
+//=======================================================================
+
+SMESH_ElementSearcher* SMESH_MeshEditor::GetElementSearcher(SMDS_ElemIteratorPtr elemIt)
+{
+ return new SMESH_ElementSearcherImpl( *GetMeshDS(), elemIt );
+}
+
+//=======================================================================
+/*!
+ * \brief Return true if the point is IN or ON of the element
+ */
+//=======================================================================
+
+bool SMESH_MeshEditor::IsOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol )
+{
+ if ( element->GetType() == SMDSAbs_Volume)
+ {
+ return SMDS_VolumeTool( element ).IsOut( point.X(), point.Y(), point.Z(), tol );
+ }
+
+ // get ordered nodes
+
+ vector< gp_XYZ > xyz;
+ vector<const SMDS_MeshNode*> nodeList;
+
+ SMDS_ElemIteratorPtr nodeIt = element->nodesIterator();
+ if ( element->IsQuadratic() ) {
+ if (const SMDS_VtkFace* f=dynamic_cast<const SMDS_VtkFace*>(element))
+ nodeIt = f->interlacedNodesElemIterator();
+ else if (const SMDS_VtkEdge* e =dynamic_cast<const SMDS_VtkEdge*>(element))
+ nodeIt = e->interlacedNodesElemIterator();
+ }
+ while ( nodeIt->more() )
+ {
+ const SMDS_MeshNode* node = cast2Node( nodeIt->next() );
+ xyz.push_back( SMESH_TNodeXYZ(node) );
+ nodeList.push_back(node);
+ }
+
+ int i, nbNodes = element->NbNodes();
+
+ if ( element->GetType() == SMDSAbs_Face ) // --------------------------------------------------
+ {
+ // compute face normal
+ gp_Vec faceNorm(0,0,0);
+ xyz.push_back( xyz.front() );
+ nodeList.push_back( nodeList.front() );
+ for ( i = 0; i < nbNodes; ++i )
+ {
+ gp_Vec edge1( xyz[i+1], xyz[i]);
+ gp_Vec edge2( xyz[i+1], xyz[(i+2)%nbNodes] );
+ faceNorm += edge1 ^ edge2;
+ }
+ double normSize = faceNorm.Magnitude();
+ if ( normSize <= tol )
+ {
+ // degenerated face: point is out if it is out of all face edges
+ for ( i = 0; i < nbNodes; ++i )
+ {
+ SMDS_LinearEdge edge( nodeList[i], nodeList[i+1] );
+ if ( !IsOut( &edge, point, tol ))
+ return false;
+ }
+ return true;
+ }
+ faceNorm /= normSize;
+
+ // check if the point lays on face plane
+ gp_Vec n2p( xyz[0], point );
+ if ( fabs( n2p * faceNorm ) > tol )
+ return true; // not on face plane
+
+ // check if point is out of face boundary:
+ // define it by closest transition of a ray point->infinity through face boundary
+ // on the face plane.
+ // First, find normal of a plane perpendicular to face plane, to be used as a cutting tool
+ // to find intersections of the ray with the boundary.
+ gp_Vec ray = n2p;
+ gp_Vec plnNorm = ray ^ faceNorm;
+ normSize = plnNorm.Magnitude();
+ if ( normSize <= tol ) return false; // point coincides with the first node
+ plnNorm /= normSize;
+ // for each node of the face, compute its signed distance to the plane
+ vector<double> dist( nbNodes + 1);
+ for ( i = 0; i < nbNodes; ++i )
+ {
+ gp_Vec n2p( xyz[i], point );
+ dist[i] = n2p * plnNorm;
+ }
+ dist.back() = dist.front();
+ // find the closest intersection
+ int iClosest = -1;
+ double rClosest, distClosest = 1e100;;
+ gp_Pnt pClosest;
+ for ( i = 0; i < nbNodes; ++i )
+ {
+ double r;
+ if ( fabs( dist[i]) < tol )
+ r = 0.;
+ else if ( fabs( dist[i+1]) < tol )
+ r = 1.;
+ else if ( dist[i] * dist[i+1] < 0 )
+ r = dist[i] / ( dist[i] - dist[i+1] );
+ else
+ continue; // no intersection
+ gp_Pnt pInt = xyz[i] * (1.-r) + xyz[i+1] * r;
+ gp_Vec p2int ( point, pInt);
+ if ( p2int * ray > -tol ) // right half-space
+ {
+ double intDist = p2int.SquareMagnitude();
+ if ( intDist < distClosest )
+ {
+ iClosest = i;
+ rClosest = r;
+ pClosest = pInt;
+ distClosest = intDist;
+ }
+ }
+ }
+ if ( iClosest < 0 )
+ return true; // no intesections - out
+
+ // analyse transition
+ gp_Vec edge( xyz[iClosest], xyz[iClosest+1] );
+ gp_Vec edgeNorm = -( edge ^ faceNorm ); // normal to intersected edge pointing out of face
+ gp_Vec p2int ( point, pClosest );
+ bool out = (edgeNorm * p2int) < -tol;
+ if ( rClosest > 0. && rClosest < 1. ) // not node intersection
+ return out;
+
+ // ray pass through a face node; analyze transition through an adjacent edge
+ gp_Pnt p1 = xyz[ (rClosest == 0.) ? ((iClosest+nbNodes-1) % nbNodes) : (iClosest+1) ];
+ gp_Pnt p2 = xyz[ (rClosest == 0.) ? iClosest : ((iClosest+2) % nbNodes) ];
+ gp_Vec edgeAdjacent( p1, p2 );
+ gp_Vec edgeNorm2 = -( edgeAdjacent ^ faceNorm );
+ bool out2 = (edgeNorm2 * p2int) < -tol;
+
+ bool covexCorner = ( edgeNorm * edgeAdjacent * (rClosest==1. ? 1. : -1.)) < 0;
+ return covexCorner ? (out || out2) : (out && out2);
+ }
+ if ( element->GetType() == SMDSAbs_Edge ) // --------------------------------------------------
+ {
+ // point is out of edge if it is NOT ON any straight part of edge
+ // (we consider quadratic edge as being composed of two straight parts)
+ for ( i = 1; i < nbNodes; ++i )
+ {
+ gp_Vec edge( xyz[i-1], xyz[i]);
+ gp_Vec n1p ( xyz[i-1], point);
+ double dist = ( edge ^ n1p ).Magnitude() / edge.Magnitude();
+ if ( dist > tol )
+ continue;
+ gp_Vec n2p( xyz[i], point );
+ if ( fabs( edge.Magnitude() - n1p.Magnitude() - n2p.Magnitude()) > tol )
+ continue;
+ return false; // point is ON this part
+ }
+ return true;
+ }
+ // Node or 0D element -------------------------------------------------------------------------
+ {
+ gp_Vec n2p ( xyz[0], point );
+ return n2p.Magnitude() <= tol;
+ }
+ return true;
+}
+
+//=======================================================================
+
+namespace
+{
+ // Position of a point relative to a segment
+ // . .
+ // . LEFT .
+ // . .
+ // VERTEX 1 o----ON-----> VERTEX 2
+ // . .
+ // . RIGHT .
+ // . .
+ enum PositionName { POS_LEFT = 1, POS_VERTEX = 2, POS_RIGHT = 4, //POS_ON = 8,
+ POS_ALL = POS_LEFT | POS_RIGHT | POS_VERTEX };
+ struct PointPos
+ {
+ PositionName _name;
+ int _index; // index of vertex or segment
+
+ PointPos( PositionName n, int i=-1 ): _name(n), _index(i) {}
+ bool operator < (const PointPos& other ) const
+ {
+ if ( _name == other._name )
+ return ( _index < 0 || other._index < 0 ) ? false : _index < other._index;
+ return _name < other._name;
+ }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Return of a point relative to a segment
+ * \param point2D - the point to analyze position of
+ * \param xyVec - end points of segments
+ * \param index0 - 0-based index of the first point of segment
+ * \param posToFindOut - flags of positions to detect
+ * \retval PointPos - point position
+ */
+ //================================================================================
+
+ PointPos getPointPosition( const gp_XY& point2D,
+ const gp_XY* segEnds,
+ const int index0 = 0,
+ const int posToFindOut = POS_ALL)
+ {
+ const gp_XY& p1 = segEnds[ index0 ];
+ const gp_XY& p2 = segEnds[ index0+1 ];
+ const gp_XY grad = p2 - p1;
+
+ if ( posToFindOut & POS_VERTEX )
+ {
+ // check if the point2D is at "vertex 1" zone
+ gp_XY pp1[2] = { p1, gp_XY( p1.X() - grad.Y(),
+ p1.Y() + grad.X() ) };
+ if ( getPointPosition( point2D, pp1, 0, POS_LEFT|POS_RIGHT )._name == POS_LEFT )
+ return PointPos( POS_VERTEX, index0 );
+
+ // check if the point2D is at "vertex 2" zone
+ gp_XY pp2[2] = { p2, gp_XY( p2.X() - grad.Y(),
+ p2.Y() + grad.X() ) };
+ if ( getPointPosition( point2D, pp2, 0, POS_LEFT|POS_RIGHT )._name == POS_RIGHT )
+ return PointPos( POS_VERTEX, index0 + 1);
+ }
+ double edgeEquation =
+ ( point2D.X() - p1.X() ) * grad.Y() - ( point2D.Y() - p1.Y() ) * grad.X();
+ return PointPos( edgeEquation < 0 ? POS_LEFT : POS_RIGHT, index0 );
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Return minimal distance from a point to a face
+ *
+ * Currently we ignore non-planarity and 2nd order of face
+ */
+//=======================================================================
+
+double SMESH_MeshEditor::GetDistance( const SMDS_MeshFace* face,
+ const gp_Pnt& point )
+{
+ double badDistance = -1;
+ if ( !face ) return badDistance;
+
+ // coordinates of nodes (medium nodes, if any, ignored)
+ typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
+ vector<gp_XYZ> xyz( TXyzIterator( face->nodesIterator()), TXyzIterator() );
+ xyz.resize( face->NbCornerNodes()+1 );
+
+ // transformation to get xyz[0] lies on the origin, xyz[1] lies on the Z axis,
+ // and xyz[2] lies in the XZ plane. This is to pass to 2D space on XZ plane.
+ gp_Trsf trsf;
+ gp_Vec OZ ( xyz[0], xyz[1] );
+ gp_Vec OX ( xyz[0], xyz[2] );
+ if ( OZ.Magnitude() < std::numeric_limits<double>::min() )
+ {
+ if ( xyz.size() < 4 ) return badDistance;
+ OZ = gp_Vec ( xyz[0], xyz[2] );
+ OX = gp_Vec ( xyz[0], xyz[3] );
+ }
+ gp_Ax3 tgtCS;
+ try {
+ tgtCS = gp_Ax3( xyz[0], OZ, OX );
+ }
+ catch ( Standard_Failure ) {
+ return badDistance;
+ }
+ trsf.SetTransformation( tgtCS );
+
+ // move all the nodes to 2D
+ vector<gp_XY> xy( xyz.size() );
+ for ( size_t i = 0;i < xyz.size()-1; ++i )
+ {
+ gp_XYZ p3d = xyz[i];
+ trsf.Transforms( p3d );
+ xy[i].SetCoord( p3d.X(), p3d.Z() );
+ }
+ xyz.back() = xyz.front();
+ xy.back() = xy.front();
+
+ // // move the point in 2D
+ gp_XYZ tmpPnt = point.XYZ();
+ trsf.Transforms( tmpPnt );
+ gp_XY point2D( tmpPnt.X(), tmpPnt.Z() );
+
+ // loop on segments of the face to analyze point position ralative to the face
+ set< PointPos > pntPosSet;
+ for ( size_t i = 1; i < xy.size(); ++i )
+ {
+ PointPos pos = getPointPosition( point2D, &xy[0], i-1 );
+ pntPosSet.insert( pos );
+ }
+
+ // compute distance
+ PointPos pos = *pntPosSet.begin();
+ // cout << "Face " << face->GetID() << " DIST: ";
+ switch ( pos._name )
+ {
+ case POS_LEFT: {
+ // point is most close to a segment
+ gp_Vec p0p1( point, xyz[ pos._index ] );
+ gp_Vec p1p2( xyz[ pos._index ], xyz[ pos._index+1 ]); // segment vector
+ p1p2.Normalize();
+ double projDist = p0p1 * p1p2; // distance projected to the segment
+ gp_Vec projVec = p1p2 * projDist;
+ gp_Vec distVec = p0p1 - projVec;
+ // cout << distVec.Magnitude() << ", SEG " << face->GetNode(pos._index)->GetID()
+ // << " - " << face->GetNodeWrap(pos._index+1)->GetID() << endl;
+ return distVec.Magnitude();
+ }
+ case POS_RIGHT: {
+ // point is inside the face
+ double distToFacePlane = tmpPnt.Y();
+ // cout << distToFacePlane << ", INSIDE " << endl;
+ return Abs( distToFacePlane );
+ }
+ case POS_VERTEX: {
+ // point is most close to a node
+ gp_Vec distVec( point, xyz[ pos._index ]);
+ // cout << distVec.Magnitude() << " VERTEX " << face->GetNode(pos._index)->GetID() << endl;
+ return distVec.Magnitude();
+ }
+ }
+ return badDistance;
+}
+
+//=======================================================================
+//function : SimplifyFace
+//purpose :
+//=======================================================================
+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;
+ //MESSAGE("node to keep " << nToKeep->GetID());
+ for ( ++nIt; nIt != nodes.end(); nIt++ ) {
+ const SMDS_MeshNode* nToRemove = *nIt;
+ nodeNodeMap.insert( TNodeNodeMap::value_type( nToRemove, nToKeep ));
+ if ( nToRemove != nToKeep ) {
+ //MESSAGE(" node to remove " << nToRemove->GetID());
+ rmNodeIds.push_back( nToRemove->GetID() );
+ AddToSameGroups( nToKeep, nToRemove, aMesh );
+ }
+
+ 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_MeshElement*>::iterator eIt = elems.begin();
+ for ( ; eIt != elems.end(); eIt++ ) {
+ const SMDS_MeshElement* elem = *eIt;
+ //MESSAGE(" ---- inverse elem on node to remove " << elem->GetID());
+ int nbNodes = elem->NbNodes();
+ int aShapeId = FindShape( elem );
+
+ set<const SMDS_MeshNode*> nodeSet;
+ vector< const SMDS_MeshNode*> curNodes( nbNodes ), uniqueNodes( nbNodes );
+ int iUnique = 0, iCur = 0, nbRepl = 0;
+ vector<int> iRepl( nbNodes );
+
+ // 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();
+ //MESSAGE("nbNodes nbUniqueNodes " << nbNodes << " " << nbUniqueNodes);
+ if ( nbNodes != nbUniqueNodes ) { // some nodes stick
+ // Polygons and Polyhedral volumes
+ if (elem->IsPoly()) {
+
+ if (elem->GetType() == SMDSAbs_Face) {
+ // Polygon
+ vector<const SMDS_MeshNode *> face_nodes (nbNodes);
+ int inode = 0;
+ for (; inode < nbNodes; inode++) {
+ face_nodes[inode] = curNodes[inode];
+ }
+
+ vector<const SMDS_MeshNode *> polygons_nodes;
+ vector<int> quantities;
+ int nbNew = SimplifyFace(face_nodes, polygons_nodes, quantities);
+ if (nbNew > 0) {
+ inode = 0;
+ for (int iface = 0; iface < nbNew; iface++) {
+ int nbNodes = quantities[iface];
+ vector<const SMDS_MeshNode *> poly_nodes (nbNodes);
+ for (int ii = 0; ii < nbNodes; ii++, inode++) {
+ poly_nodes[ii] = polygons_nodes[inode];
+ }
+ SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
+ myLastCreatedElems.Append(newElem);
+ if (aShapeId)
+ aMesh->SetMeshElementOnShape(newElem, aShapeId);
+ }
+
+ MESSAGE("ChangeElementNodes MergeNodes Polygon");
+ //aMesh->ChangeElementNodes(elem, &polygons_nodes[inode], quantities[nbNew - 1]);
+ vector<const SMDS_MeshNode *> polynodes(polygons_nodes.begin()+inode,polygons_nodes.end());
+ int quid =0;
+ if (nbNew > 0) quid = nbNew - 1;
+ vector<int> newquant(quantities.begin()+quid, quantities.end());
+ const SMDS_MeshElement* newElem = 0;
+ newElem = aMesh->AddPolyhedralVolume(polynodes, newquant);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ rmElemIds.push_back(elem->GetID());
+ }
+ else {
+ rmElemIds.push_back(elem->GetID());
+ }
+
+ }
+ 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)
+ {
+ MESSAGE("ChangeElementNodes MergeNodes Polyhedron");
+ //aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
+ const SMDS_MeshElement* newElem = 0;
+ 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 elem remains valid after sticking nodes
+ if (elem->IsPoly() && elem->GetType() == SMDSAbs_Volume)
+ {
+ // Change nodes of polyedre
+ 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 (nbFaces);
+
+ for (int iface = 1; iface <= nbFaces; iface++) {
+ int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
+ quantities[iface - 1] = nbFaceNodes;
+
+ for (inode = 1; inode <= nbFaceNodes; inode++) {
+ const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
+
+ TNodeNodeMap::iterator nnIt = nodeNodeMap.find( curNode );
+ if (nnIt != nodeNodeMap.end()) { // curNode sticks
+ curNode = (*nnIt).second;
+ }
+ poly_nodes.push_back(curNode);
+ }
+ }
+ aMesh->ChangePolyhedronNodes( elem, poly_nodes, quantities );
+ }
+ }
+ else // replace non-polyhedron elements
+ {
+ const SMDSAbs_ElementType etyp = elem->GetType();
+ const int elemId = elem->GetID();
+ const bool isPoly = (elem->GetEntityType() == SMDSEntity_Polygon);
+ uniqueNodes.resize(nbUniqueNodes);
+
+ SMESHDS_SubMesh * sm = aShapeId > 0 ? aMesh->MeshElements(aShapeId) : 0;
+
+ aMesh->RemoveFreeElement(elem, sm, /*fromGroups=*/false);
+ SMDS_MeshElement* newElem = this->AddElement(uniqueNodes, etyp, isPoly, elemId);
+ 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 );
+
+}
+
+
+// ========================================================
+// 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(set<const SMDS_MeshElement*> & theElements,
+ TListOfListOfElementsID & theGroupsOfElementsID)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ typedef set<const SMDS_MeshElement*> TElemsSet;
+ typedef map< SortableElement, int > TMapOfNodeSet;
+ typedef list<int> TGroupOfElems;
+
+ TElemsSet elems;
+ if ( theElements.empty() )
+ { // get all elements in the mesh
+ SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
+ while ( eIt->more() )
+ elems.insert( elems.end(), eIt->next());
+ }
+ else
+ elems = theElements;
+
+ vector< TGroupOfElems > arrayOfGroups;
+ TGroupOfElems groupOfElems;
+ TMapOfNodeSet mapOfNodeSet;
+
+ TElemsSet::iterator elemIt = elems.begin();
+ for ( int i = 0, j=0; elemIt != elems.end(); ++elemIt, ++j ) {
+ const SMDS_MeshElement* curElem = *elemIt;
+ SortableElement SE(curElem);
+ int ind = -1;
+ // check uniqueness
+ pair< TMapOfNodeSet::iterator, bool> pp = mapOfNodeSet.insert(make_pair(SE, i));
+ if( !(pp.second) ) {
+ TMapOfNodeSet::iterator& itSE = pp.first;
+ ind = (*itSE).second;
+ arrayOfGroups[ind].push_back(curElem->GetID());
+ }
+ else {
+ groupOfElems.clear();
+ groupOfElems.push_back(curElem->GetID());
+ arrayOfGroups.push_back(groupOfElems);
+ i++;
+ }
+ }
+
+ vector< TGroupOfElems >::iterator groupIt = arrayOfGroups.begin();
+ for ( ; groupIt != arrayOfGroups.end(); ++groupIt ) {
+ groupOfElems = *groupIt;
+ if ( groupOfElems.size() > 1 ) {
+ groupOfElems.sort();
+ theGroupsOfElementsID.push_back(groupOfElems);
+ }
+ }
+}
+
+//=======================================================================
+//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()
+{
+ set<const SMDS_MeshElement*> aMeshElements; /* empty input -
+ to merge equal elements in the whole mesh */
+ TListOfListOfElementsID aGroupsOfElementsID;
+ FindEqualElements(aMeshElements, aGroupsOfElementsID);
+ MergeElements(aGroupsOfElementsID);
+}
+
+//=======================================================================
+//function : FindFaceInSet
+//purpose : Return a face having linked nodes n1 and n2 and which is
+// - not in avoidSet,
+// - in elemSet provided that !elemSet.empty()
+// i1 and i2 optionally returns indices of n1 and n2
+//=======================================================================
+
+const SMDS_MeshElement*
+SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ const TIDSortedElemSet& elemSet,
+ const TIDSortedElemSet& avoidSet,
+ int* n1ind,
+ int* n2ind)
+
+{
+ int i1, i2;
+ const SMDS_MeshElement* face = 0;
+
+ SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face);
+ //MESSAGE("n1->GetInverseElementIterator(SMDSAbs_Face) " << invElemIt);
+ while ( invElemIt->more() && !face ) // loop on inverse faces of n1
+ {
+ //MESSAGE("in while ( invElemIt->more() && !face )");
+ const SMDS_MeshElement* elem = invElemIt->next();
+ if (avoidSet.count( elem ))
+ continue;
+ if ( !elemSet.empty() && !elemSet.count( elem ))
+ continue;
+ // index of n1
+ i1 = elem->GetNodeIndex( n1 );
+ // find a n2 linked to n1
+ int nbN = elem->IsQuadratic() ? elem->NbNodes()/2 : elem->NbNodes();
+ for ( int di = -1; di < 2 && !face; di += 2 )
+ {
+ i2 = (i1+di+nbN) % nbN;
+ if ( elem->GetNode( i2 ) == n2 )
+ face = elem;
+ }
+ if ( !face && elem->IsQuadratic())
+ {
+ // analysis for quadratic elements using all nodes
+ 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();
+ const SMDS_MeshNode* prevN = cast2Node( anIter->next() );
+ for ( i1 = -1, i2 = 0; anIter->more() && !face; i1++, i2++ )
+ {
+ const SMDS_MeshNode* n = cast2Node( anIter->next() );
+ if ( n1 == prevN && n2 == n )
+ {
+ face = elem;
+ }
+ else if ( n2 == prevN && n1 == n )
+ {
+ face = elem; swap( i1, i2 );
+ }
+ prevN = n;
+ }
+ }
+ }
+ if ( n1ind ) *n1ind = i1;
+ if ( n2ind ) *n2ind = i2;
+ return face;
+}
+
+//=======================================================================
+//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_MeshEditor::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 ( ! 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 ;
+
+ // 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;
+
+ 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() && !isFLN ) {
+ const SMDS_MeshNode* n = cast2Node(anIter->next());
+ poly_nodes[iNode++] = n;
+ if (n == nodes[il1]) {
+ isFLN = true;
+ }
+ }
+ // 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 ( anIter->more() ) {
+ poly_nodes[iNode++] = cast2Node(anIter->next());
+ }
+ }
+ else {
+ SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
+ while ( nodeIt->more() && !isFLN ) {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ poly_nodes[iNode++] = n;
+ if (n == nodes[il1]) {
+ isFLN = true;
+ }
+ }
+ // 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;
+ }
+ }
+
+ // edit or replace the face
+ SMESHDS_Mesh *aMesh = GetMeshDS();
+
+ if (theFace->IsPoly()) {
+ aMesh->ChangePolygonNodes(theFace, poly_nodes);
+ }
+ else {
+ int aShapeId = FindShape( theFace );
+
+ SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+
+ aMesh->RemoveElement(theFace);
+ }
+ return;
+ }
+
+ SMESHDS_Mesh *aMesh = GetMeshDS();
+ 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;
+ 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
+ int aShapeId = FindShape( theFace );
+
+ i1 = 0; i2 = 1;
+ for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ ) {
+ SMDS_MeshElement* newElem = 0;
+ if ( iSplit == iBestQuad )
+ newElem = aMesh->AddFace (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ i3 ],
+ nodes[ i4 ]);
+ else
+ newElem = aMesh->AddFace (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ iSplit < iBestQuad ? i4 : i3 ]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ }
+
+ // change nodes of theFace
+ const SMDS_MeshNode* newNodes[ 4 ];
+ newNodes[ 0 ] = linkNodes[ i1 ];
+ newNodes[ 1 ] = linkNodes[ i2 ];
+ newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
+ newNodes[ 3 ] = nodes[ i4 ];
+ //aMesh->ChangeElementNodes( theFace, newNodes, iSplit == iBestQuad ? 4 : 3 );
+ const SMDS_MeshElement* newElem = 0;
+ if (iSplit == iBestQuad)
+ newElem = aMesh->AddFace( newNodes[0], newNodes[1], newNodes[2], newNodes[3] );
+ else
+ newElem = aMesh->AddFace( newNodes[0], newNodes[1], newNodes[2] );
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+} // 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 aShapeId = FindShape( theFace );
+
+ int n1,n2,n3;
+ if(nbFaceNodes==6) { // quadratic triangle
+ SMDS_MeshElement* newElem =
+ aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ if(theFace->IsMediumNode(nodes[il1])) {
+ // create quadrangle
+ newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[5]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ n1 = 1;
+ n2 = 2;
+ n3 = 3;
+ }
+ else {
+ // create quadrangle
+ newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[5]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ n1 = 0;
+ n2 = 1;
+ n3 = 5;
+ }
+ }
+ else { // nbFaceNodes==8 - quadratic quadrangle
+ SMDS_MeshElement* newElem =
+ aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ newElem = aMesh->AddFace(nodes[5],nodes[6],nodes[7]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ newElem = aMesh->AddFace(nodes[5],nodes[7],nodes[3]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ if(theFace->IsMediumNode(nodes[il1])) {
+ // create quadrangle
+ newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[7]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ n1 = 1;
+ n2 = 2;
+ n3 = 3;
+ }
+ else {
+ // create quadrangle
+ newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[7]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ n1 = 0;
+ n2 = 1;
+ n3 = 7;
+ }
+ }
+ // create needed triangles using n1,n2,n3 and inserted nodes
+ int nbn = 2 + aNodesToInsert.size();
+ //const SMDS_MeshNode* aNodes[nbn];
+ 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++) {
+ SMDS_MeshElement* newElem =
+ aMesh->AddFace(aNodes[i-1],aNodes[i],nodes[n3]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ }
+ }
+ // remove old face
+ aMesh->RemoveElement(theFace);
+}
+
+//=======================================================================
+//function : UpdateVolumes
+//purpose :
+//=======================================================================
+void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode* theBetweenNode1,
+ const SMDS_MeshNode* theBetweenNode2,
+ list<const SMDS_MeshNode*>& theNodesToInsert)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ 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 or update the volume
+ SMESHDS_Mesh *aMesh = GetMeshDS();
+
+ if (elem->IsPoly()) {
+ aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
+
+ }
+ else {
+ int aShapeId = FindShape( elem );
+
+ SMDS_MeshElement* newElem =
+ aMesh->AddPolyhedralVolume(poly_nodes, quantities);
+ myLastCreatedElems.Append(newElem);
+ if (aShapeId && newElem)
+ aMesh->SetMeshElementOnShape(newElem, aShapeId);
+
+ 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 submesh to quadratic
+ * \return int - nb of checked elements
+ */
+//=======================================================================
+
+int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm,
+ SMESH_MesherHelper& theHelper,
+ const bool theForce3d)
+{
+ 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 || elem->IsQuadratic() ) continue;
+
+ // get elem data needed to re-create it
+ //
+ const int id = elem->GetID();
+ const int nbNodes = elem->NbNodes();
+ const SMDSAbs_ElementType aType = elem->GetType();
+ const SMDSAbs_EntityType aGeomType = elem->GetEntityType();
+ nodes.assign(elem->begin_nodes(), elem->end_nodes());
+ if ( aGeomType == SMDSEntity_Polyhedra )
+ nbNodeInFaces = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
+ else if ( aGeomType == SMDSEntity_Hexagonal_Prism )
+ volumeToPolyhedron( elem, nodes, nbNodeInFaces );
+
+ // remove a linear element
+ GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
+
+ 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);
+ continue;
+ }
+ 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:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], id, theForce3d);
+ break;
+ case SMDSEntity_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 )
+ theSm->AddElement( NewElem );
+ }
+ return nbElem;
+}
+
+//=======================================================================
+//function : ConvertToQuadratic
+//purpose :
+//=======================================================================
+
+void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
+{
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+
+ SMESH_MesherHelper aHelper(*myMesh);
+ aHelper.SetIsQuadratic( true );
+
+ 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);
+ }
+ }
+ }
+ }
+ int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
+ if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
+ {
+ SMESHDS_SubMesh *smDS = 0;
+ SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator();
+ while(aEdgeItr->more())
+ {
+ const SMDS_MeshEdge* edge = aEdgeItr->next();
+ if(edge && !edge->IsQuadratic())
+ {
+ int id = edge->GetID();
+ //MESSAGE("edge->GetID() " << id);
+ 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());
+ }
+ }
+ SMDS_FaceIteratorPtr aFaceItr = meshDS->facesIterator();
+ while(aFaceItr->more())
+ {
+ const SMDS_MeshFace* face = aFaceItr->next();
+ if(!face || face->IsQuadratic() ) continue;
+
+ const int id = face->GetID();
+ const SMDSAbs_EntityType type = face->GetEntityType();
+ vector<const SMDS_MeshNode *> nodes ( face->begin_nodes(), face->end_nodes());
+
+ meshDS->RemoveFreeElement(face, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshFace * NewFace = 0;
+ switch( type )
+ {
+ case SMDSEntity_Triangle:
+ NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
+ break;
+ case SMDSEntity_Quadrangle:
+ NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ default:
+ NewFace = aHelper.AddPolygonalFace(nodes, id, theForce3d);
+ }
+ ReplaceElemInGroups( face, NewFace, GetMeshDS());
+ }
+ vector<int> nbNodeInFaces;
+ SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator();
+ while(aVolumeItr->more())
+ {
+ const SMDS_MeshVolume* volume = aVolumeItr->next();
+ if(!volume || volume->IsQuadratic() ) continue;
+
+ const int id = volume->GetID();
+ const SMDSAbs_EntityType type = volume->GetEntityType();
+ vector<const SMDS_MeshNode *> nodes (volume->begin_nodes(), volume->end_nodes());
+ if ( type == SMDSEntity_Polyhedra )
+ nbNodeInFaces = static_cast<const SMDS_VtkVolume* >(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:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+ break;
+ case SMDSEntity_Pyramid:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+ nodes[3], nodes[4], id, theForce3d);
+ break;
+ case SMDSEntity_Penta:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+ nodes[3], nodes[4], nodes[5], id, theForce3d);
+ 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);
+ }
+}
+
+//================================================================================
+/*!
+ * \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)
+{
+ 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();
+ if ( e->IsQuadratic() )
+ {
+ quadAdjacentElems[ e->GetType() ].insert( e );
+ continue;
+ }
+ if ( e->GetType() >= elemType )
+ {
+ continue; // same type of more complex linear element
+ }
+
+ if ( !checkedAdjacentElems[ e->GetType() ].insert( e ).second )
+ continue; // e is already checked
+
+ // check nodes
+ bool allIn = true;
+ SMDS_ElemIteratorPtr nodeIt = e->nodesIterator();
+ while ( nodeIt->more() && allIn )
+ allIn = allNodes.count( cast2Node( nodeIt->next() ));
+ if ( allIn )
+ theElements.insert(e );
+ }
+ }
+
+ SMESH_MesherHelper helper(*myMesh);
+ helper.SetIsQuadratic( true );
+
+ // 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 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;
+ if( elem->IsQuadratic() || elem->NbNodes() < 2 || elem->IsPoly() )
+ continue;
+
+ const int id = elem->GetID();
+ const SMDSAbs_ElementType type = elem->GetType();
+ vector<const SMDS_MeshNode *> nodes ( elem->begin_nodes(), elem->end_nodes());
+
+ if ( !smDS || !smDS->Contains( elem ))
+ smDS = meshDS->MeshElements( elem->getshapeId() );
+ meshDS->RemoveFreeElement(elem, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshElement * newElem = 0;
+ switch( nodes.size() )
+ {
+ 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 );
+ }
+
+ if ( !theForce3d && !getenv("NO_FixQuadraticElements"))
+ { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
+ helper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
+ helper.FixQuadraticElements( 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();
+
+ while( theItr->more() )
+ {
+ const SMDS_MeshElement* elem = theItr->next();
+ nbElem++;
+ if( elem && elem->IsQuadratic())
+ {
+ int id = elem->GetID();
+ int nbCornerNodes = elem->NbCornerNodes();
+ SMDSAbs_ElementType aType = elem->GetType();
+
+ vector<const SMDS_MeshNode *> nodes( elem->begin_nodes(), elem->end_nodes() );
+
+ //remove a quadratic element
+ if ( !theSm || !theSm->Contains( elem ))
+ theSm = meshDS->MeshElements( elem->getshapeId() );
+ meshDS->RemoveFreeElement( elem, theSm, /*fromGroups=*/false );
+
+ // remove medium nodes
+ for ( unsigned 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, aType, false, id );
+ 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
+ typedef SMDS_SetIterator<const SMDS_MeshElement*, TIDSortedElemSet::iterator> TSetIterator;
+ SMDS_ElemIteratorPtr elemIt( new TSetIterator( theElements.begin(), theElements.end() ));
+ 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_Algo::GetCommonNodes( e, eComplex ).size();
+ if ( nbCommonNodes == e->NbNodes())
+ {
+ complexFound = true;
+ type = SMDSAbs_NbElementTypes; // to quit from the outer loop
+ break;
+ }
+ }
+ }
+ }
+ if ( !complexFound )
+ moreElemsToConvert.insert( e );
+ }
+ }
+ }
+ elemIt = SMDS_ElemIteratorPtr
+ (new TSetIterator( moreElemsToConvert.begin(), moreElemsToConvert.end() ));
+ 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)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ MESSAGE ("::::SewSideElements()");
+ 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 algoritm 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 );