+ 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;
+}
+
+//=======================================================================
+//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
+ iRepl[ nbRepl++ ] = iCur;
+ }
+ curNodes[ iCur ] = n;
+ bool isUnique = nodeSet.insert( n ).second;
+ if ( isUnique )
+ uniqueNodes[ iUnique++ ] = n;
+ 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;
+ }
+
+ // 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 == 6 ) {
+ //////////////////////// ---> 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;
+ iUnique = 2; // reverse a tetrahedron bottom
+ for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
+ if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
+ nbStick++;
+ else if ( iUnique >= 0 )
+ uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
+ }
+ 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 == 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 )
+ break;
+ } // HEXAHEDRON
+
+ default:
+ isOk = false;
+ } // switch ( nbNodes )
+
+ } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
+
+ if ( isOk ) {
+ 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 {
+ //int elemId = elem->GetID();
+ //MESSAGE("Change regular element or polygon " << elemId);
+ SMDSAbs_ElementType etyp = elem->GetType();
+ uniqueNodes.resize(nbUniqueNodes);
+ SMDS_MeshElement* newElem = 0;
+ if (elem->GetEntityType() == SMDSEntity_Polygon)
+ newElem = this->AddElement(uniqueNodes, etyp, true);
+ else
+ newElem = this->AddElement(uniqueNodes, etyp, false);
+ if (newElem)
+ {
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ }
+ aMesh->RemoveElement(elem);
+ }
+ }
+ else {
+ // Remove invalid regular element or invalid polygon
+ //MESSAGE("Remove invalid " << elem->GetID());
+ 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);
+ }
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Convert elements contained in a submesh to quadratic
+ * \retval 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;
+ SMDS_ElemIteratorPtr ElemItr = theSm->GetElements();
+ while(ElemItr->more())
+ {
+ nbElem++;
+ const SMDS_MeshElement* elem = ElemItr->next();
+ if( !elem || elem->IsQuadratic() ) continue;
+
+ int id = elem->GetID();
+ //MESSAGE("elem " << id);
+ id = 0; // get a free number for new elements
+ int nbNodes = elem->NbNodes();
+ SMDSAbs_ElementType aType = elem->GetType();
+
+ vector<const SMDS_MeshNode *> nodes (elem->begin_nodes(), elem->end_nodes());
+ if ( elem->GetEntityType() == SMDSEntity_Polyhedra )
+ nbNodeInFaces = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
+
+ 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(nbNodes)
+ {
+ case 4:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ case 5:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], id, theForce3d);
+ break;
+ case 6:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], id, theForce3d);
+ break;
+ case 8:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+ break;
+ default:
+ NewElem = theHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
+ }
+ break;
+ }
+ default :
+ continue;
+ }
+ ReplaceElemInGroups( elem, NewElem, GetMeshDS());
+ if( NewElem )
+ theSm->AddElement( NewElem );
+
+ GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
+ }
+// if (!GetMeshDS()->isCompacted())
+// GetMeshDS()->compactMesh();
+ 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;
+
+ int id = face->GetID();
+ int nbNodes = face->NbNodes();
+ vector<const SMDS_MeshNode *> nodes ( face->begin_nodes(), face->end_nodes());
+
+ meshDS->RemoveFreeElement(face, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshFace * NewFace = 0;
+ switch(nbNodes)
+ {
+ case 3:
+ NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
+ break;
+ case 4:
+ 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;
+
+ int id = volume->GetID();
+ int nbNodes = volume->NbNodes();
+ vector<const SMDS_MeshNode *> nodes (volume->begin_nodes(), volume->end_nodes());
+ if ( volume->GetEntityType() == SMDSEntity_Polyhedra )
+ nbNodeInFaces = static_cast<const SMDS_VtkVolume* >(volume)->GetQuantities();
+
+ meshDS->RemoveFreeElement(volume, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshVolume * NewVolume = 0;
+ switch(nbNodes)
+ {
+ case 4:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+ nodes[3], id, theForce3d );
+ break;
+ case 5:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+ nodes[3], nodes[4], id, theForce3d);
+ break;
+ case 6:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+ nodes[3], nodes[4], nodes[5], id, theForce3d);
+ break;
+ case 8:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+ break;
+ default:
+ NewVolume = aHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
+ }
+ ReplaceElemInGroups(volume, NewVolume, meshDS);
+ }
+ }
+
+ if ( !theForce3d && !getenv("NO_FixQuadraticElements"))
+ { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
+ aHelper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
+ aHelper.FixQuadraticElements();
+ }
+ if (!GetMeshDS()->isCompacted())
+ GetMeshDS()->compactMesh();
+}
+
+//=======================================================================
+/*!
+ * \brief Convert quadratic elements to linear ones and remove quadratic nodes
+ * \retval 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();
+ const bool notFromGroups = false;
+
+ while( theItr->more() )
+ {
+ const SMDS_MeshElement* elem = theItr->next();
+ nbElem++;
+ if( elem && elem->IsQuadratic())
+ {
+ int id = elem->GetID();
+ int nbNodes = elem->NbNodes();
+ vector<const SMDS_MeshNode *> nodes, mediumNodes;
+ nodes.reserve( nbNodes );
+ mediumNodes.reserve( nbNodes );
+
+ for(int i = 0; i < nbNodes; i++)
+ {
+ const SMDS_MeshNode* n = elem->GetNode(i);
+
+ if( elem->IsMediumNode( n ) )
+ mediumNodes.push_back( n );
+ else
+ nodes.push_back( n );
+ }
+ if( nodes.empty() ) continue;
+ SMDSAbs_ElementType aType = elem->GetType();
+
+ //remove old quadratic element
+ meshDS->RemoveFreeElement( elem, theSm, notFromGroups );
+
+ SMDS_MeshElement * NewElem = AddElement( nodes, aType, false, id );
+ ReplaceElemInGroups(elem, NewElem, meshDS);
+ if( theSm && NewElem )
+ theSm->AddElement( NewElem );
+
+ // remove medium nodes
+ vector<const SMDS_MeshNode*>::iterator nIt = mediumNodes.begin();
+ for ( ; nIt != mediumNodes.end(); ++nIt ) {
+ const SMDS_MeshNode* n = *nIt;
+ if ( n->NbInverseElements() == 0 ) {
+ if ( n->getshapeId() != theShapeID )
+ meshDS->RemoveFreeNode( n, meshDS->MeshElements
+ ( n->getshapeId() ));
+ else
+ meshDS->RemoveFreeNode( n, theSm );
+ }
+ }
+ }
+ }
+ 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;
+}
+
+//=======================================================================
+//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
+ set<const SMDS_MeshElement*> faceSet1, faceSet2;
+ set<const SMDS_MeshElement*> volSet1, volSet2;
+ set<const SMDS_MeshNode*> nodeSet1, nodeSet2;
+ set<const SMDS_MeshElement*> * faceSetPtr[] = { &faceSet1, &faceSet2 };
+ set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
+ 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 ];
+ set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
+ set<const SMDS_MeshElement*> * volSet = volSetPtr [ iSide ];
+ set<const SMDS_MeshElement*>::iterator vIt;
+ TIDSortedElemSet::iterator eIt;
+ set<const SMDS_MeshNode*>::iterator nIt;
+
+ // check that given nodes belong to given elements
+ const SMDS_MeshNode* n1 = ( iSide == 0 ) ? theFirstNode1 : theFirstNode2;
+ const SMDS_MeshNode* n2 = ( iSide == 0 ) ? theSecondNode1 : theSecondNode2;
+ int firstIndex = -1, secondIndex = -1;
+ for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
+ const SMDS_MeshElement* elem = *eIt;
+ if ( firstIndex < 0 ) firstIndex = elem->GetNodeIndex( n1 );
+ if ( secondIndex < 0 ) secondIndex = elem->GetNodeIndex( n2 );
+ if ( firstIndex > -1 && secondIndex > -1 ) break;
+ }
+ if ( firstIndex < 0 || secondIndex < 0 ) {
+ // we can simply return until temporary faces created
+ return (iSide == 0 ) ? SEW_BAD_SIDE1_NODES : SEW_BAD_SIDE2_NODES;
+ }
+
+ // -----------------------------------------------------------
+ // 1a. Collect nodes of existing faces
+ // and build set of face nodes in order to detect missing
+ // faces corresponding to sides of volumes
+ // -----------------------------------------------------------
+
+ set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
+
+ // loop on the given element of a side
+ for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
+ //const SMDS_MeshElement* elem = *eIt;
+ const SMDS_MeshElement* elem = *eIt;
+ if ( elem->GetType() == SMDSAbs_Face ) {
+ faceSet->insert( elem );
+ set <const SMDS_MeshNode*> faceNodeSet;
+ SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+ while ( nodeIt->more() ) {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ nodeSet->insert( n );
+ faceNodeSet.insert( n );
+ }
+ setOfFaceNodeSet.insert( faceNodeSet );
+ }
+ else if ( elem->GetType() == SMDSAbs_Volume )
+ volSet->insert( elem );
+ }
+ // ------------------------------------------------------------------------------
+ // 1b. Complete set of faces: find missing faces whose nodes are in set of nodes
+ // ------------------------------------------------------------------------------
+
+ for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
+ SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() ) { // loop on faces sharing a node
+ const SMDS_MeshElement* f = fIt->next();
+ if ( faceSet->find( f ) == faceSet->end() ) {
+ // check if all nodes are in nodeSet and
+ // complete setOfFaceNodeSet if they are
+ set <const SMDS_MeshNode*> faceNodeSet;
+ SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
+ bool allInSet = true;
+ while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ if ( nodeSet->find( n ) == nodeSet->end() )
+ allInSet = false;
+ else
+ faceNodeSet.insert( n );
+ }
+ if ( allInSet ) {
+ faceSet->insert( f );
+ setOfFaceNodeSet.insert( faceNodeSet );
+ }
+ }
+ }
+ }
+
+ // -------------------------------------------------------------------------
+ // 1c. Create temporary faces representing sides of volumes if correspondent
+ // face does not exist
+ // -------------------------------------------------------------------------
+
+ if ( !volSet->empty() ) {
+ //int nodeSetSize = nodeSet->size();
+
+ // loop on given volumes
+ for ( vIt = volSet->begin(); vIt != volSet->end(); vIt++ ) {
+ SMDS_VolumeTool vol (*vIt);
+ // loop on volume faces: find free faces
+ // --------------------------------------
+ list<const SMDS_MeshElement* > freeFaceList;
+ for ( iFace = 0; iFace < vol.NbFaces(); iFace++ ) {
+ if ( !vol.IsFreeFace( iFace ))
+ continue;
+ // check if there is already a face with same nodes in a face set
+ const SMDS_MeshElement* aFreeFace = 0;
+ const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iFace );
+ int nbNodes = vol.NbFaceNodes( iFace );
+ set <const SMDS_MeshNode*> faceNodeSet;
+ vol.GetFaceNodes( iFace, faceNodeSet );
+ bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
+ if ( isNewFace ) {
+ // no such a face is given but it still can exist, check it
+ if ( nbNodes == 3 ) {
+ aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2] );
+ }
+ else if ( nbNodes == 4 ) {
+ aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
+ }
+ else {
+ vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
+ aFreeFace = aMesh->FindFace(poly_nodes);
+ }
+ }
+ if ( !aFreeFace ) {
+ // create a temporary face
+ if ( nbNodes == 3 ) {
+ //aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
+ aFreeFace = aMesh->AddFace( fNodes[0],fNodes[1],fNodes[2] );
+ }
+ else if ( nbNodes == 4 ) {
+ //aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
+ aFreeFace = aMesh->AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
+ }
+ else {
+ vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
+ //aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
+ aFreeFace = aMesh->AddPolygonalFace(poly_nodes);
+ }
+ }
+ if ( aFreeFace ) {
+ freeFaceList.push_back( aFreeFace );
+ tempFaceList.push_back( aFreeFace );
+ }
+
+ } // loop on faces of a volume
+
+ // choose one of several free faces
+ // --------------------------------------
+ if ( freeFaceList.size() > 1 ) {
+ // choose a face having max nb of nodes shared by other elems of a side
+ int maxNbNodes = -1/*, nbExcludedFaces = 0*/;
+ list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
+ while ( fIt != freeFaceList.end() ) { // loop on free faces
+ int nbSharedNodes = 0;
+ SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
+ while ( nodeIt->more() ) { // loop on free face nodes
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
+ while ( invElemIt->more() ) {
+ const SMDS_MeshElement* e = invElemIt->next();
+ if ( faceSet->find( e ) != faceSet->end() )
+ nbSharedNodes++;
+ if ( elemSet->find( e ) != elemSet->end() )
+ nbSharedNodes++;
+ }
+ }
+ if ( nbSharedNodes >= maxNbNodes ) {
+ maxNbNodes = nbSharedNodes;
+ fIt++;
+ }
+ else
+ freeFaceList.erase( fIt++ ); // here fIt++ occurs before erase
+ }
+ if ( freeFaceList.size() > 1 )
+ {
+ // could not choose one face, use another way
+ // choose a face most close to the bary center of the opposite side
+ gp_XYZ aBC( 0., 0., 0. );
+ set <const SMDS_MeshNode*> addedNodes;
+ TIDSortedElemSet * elemSet2 = elemSetPtr[ 1 - iSide ];
+ eIt = elemSet2->begin();
+ for ( eIt = elemSet2->begin(); eIt != elemSet2->end(); eIt++ ) {
+ SMDS_ElemIteratorPtr nodeIt = (*eIt)->nodesIterator();
+ while ( nodeIt->more() ) { // loop on free face nodes
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ if ( addedNodes.insert( n ).second )
+ aBC += gp_XYZ( n->X(),n->Y(),n->Z() );
+ }
+ }
+ aBC /= addedNodes.size();
+ double minDist = DBL_MAX;
+ fIt = freeFaceList.begin();
+ while ( fIt != freeFaceList.end() ) { // loop on free faces
+ double dist = 0;
+ SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
+ while ( nodeIt->more() ) { // loop on free face nodes
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ gp_XYZ p( n->X(),n->Y(),n->Z() );
+ dist += ( aBC - p ).SquareModulus();
+ }
+ if ( dist < minDist ) {
+ minDist = dist;
+ freeFaceList.erase( freeFaceList.begin(), fIt++ );
+ }
+ else
+ fIt = freeFaceList.erase( fIt++ );
+ }
+ }
+ } // choose one of several free faces of a volume
+
+ if ( freeFaceList.size() == 1 ) {
+ const SMDS_MeshElement* aFreeFace = freeFaceList.front();
+ faceSet->insert( aFreeFace );
+ // complete a node set with nodes of a found free face
+ // for ( iNode = 0; iNode < ; iNode++ )
+ // nodeSet->insert( fNodes[ iNode ] );
+ }
+
+ } // loop on volumes of a side
+
+ // // complete a set of faces if new nodes in a nodeSet appeared
+ // // ----------------------------------------------------------
+ // if ( nodeSetSize != nodeSet->size() ) {
+ // for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
+ // SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
+ // while ( fIt->more() ) { // loop on faces sharing a node
+ // const SMDS_MeshElement* f = fIt->next();
+ // if ( faceSet->find( f ) == faceSet->end() ) {
+ // // check if all nodes are in nodeSet and
+ // // complete setOfFaceNodeSet if they are
+ // set <const SMDS_MeshNode*> faceNodeSet;
+ // SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
+ // bool allInSet = true;
+ // while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
+ // const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ // if ( nodeSet->find( n ) == nodeSet->end() )
+ // allInSet = false;
+ // else
+ // faceNodeSet.insert( n );
+ // }
+ // if ( allInSet ) {
+ // faceSet->insert( f );
+ // setOfFaceNodeSet.insert( faceNodeSet );
+ // }
+ // }
+ // }
+ // }
+ // }
+ } // Create temporary faces, if there are volumes given
+ } // loop on sides
+
+ if ( faceSet1.size() != faceSet2.size() ) {
+ // delete temporary faces: they are in reverseElements of actual nodes
+// SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
+// while ( tmpFaceIt->more() )
+// aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
+// list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
+// for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
+// aMesh->RemoveElement(*tmpFaceIt);
+ MESSAGE("Diff nb of faces");
+ return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ }
+
+ // ============================================================
+ // 2. Find nodes to merge:
+ // bind a node to remove to a node to put instead
+ // ============================================================
+
+ TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
+ if ( theFirstNode1 != theFirstNode2 )
+ nReplaceMap.insert( TNodeNodeMap::value_type( theFirstNode1, theFirstNode2 ));
+ if ( theSecondNode1 != theSecondNode2 )
+ nReplaceMap.insert( TNodeNodeMap::value_type( theSecondNode1, theSecondNode2 ));
+
+ LinkID_Gen aLinkID_Gen( GetMeshDS() );
+ set< long > linkIdSet; // links to process
+ linkIdSet.insert( aLinkID_Gen.GetLinkID( theFirstNode1, theSecondNode1 ));
+
+ typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > NLink;
+ list< NLink > linkList[2];
+ linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
+ linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
+ // loop on links in linkList; find faces by links and append links
+ // of the found faces to linkList
+ list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
+ for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
+ NLink link[] = { *linkIt[0], *linkIt[1] };
+ long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
+ if ( linkIdSet.find( linkID ) == linkIdSet.end() )
+ continue;
+
+ // by links, find faces in the face sets,
+ // and find indices of link nodes in the found faces;
+ // in a face set, there is only one or no face sharing a link
+ // ---------------------------------------------------------------
+
+ const SMDS_MeshElement* face[] = { 0, 0 };
+ //const SMDS_MeshNode* faceNodes[ 2 ][ 5 ];
+ vector<const SMDS_MeshNode*> fnodes1(9);
+ vector<const SMDS_MeshNode*> fnodes2(9);
+ //const SMDS_MeshNode* notLinkNodes[ 2 ][ 2 ] = {{ 0, 0 },{ 0, 0 }} ;
+ vector<const SMDS_MeshNode*> notLinkNodes1(6);
+ vector<const SMDS_MeshNode*> notLinkNodes2(6);
+ int iLinkNode[2][2];
+ for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
+ const SMDS_MeshNode* n1 = link[iSide].first;
+ const SMDS_MeshNode* n2 = link[iSide].second;
+ set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
+ set< const SMDS_MeshElement* > fMap;
+ for ( int i = 0; i < 2; i++ ) { // loop on 2 nodes of a link
+ const SMDS_MeshNode* n = i ? n1 : n2; // a node of a link
+ SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() ) { // loop on faces sharing a node
+ const SMDS_MeshElement* f = fIt->next();
+ if (faceSet->find( f ) != faceSet->end() && // f is in face set
+ ! fMap.insert( f ).second ) // f encounters twice
+ {
+ if ( face[ iSide ] ) {
+ MESSAGE( "2 faces per link " );
+ aResult = iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES;
+ break;
+ }
+ face[ iSide ] = f;
+ faceSet->erase( f );
+ // get face nodes and find ones of a link
+ iNode = 0;
+ int nbl = -1;
+ if(f->IsPoly()) {
+ if(iSide==0) {
+ fnodes1.resize(f->NbNodes()+1);
+ notLinkNodes1.resize(f->NbNodes()-2);
+ }
+ else {
+ fnodes2.resize(f->NbNodes()+1);
+ notLinkNodes2.resize(f->NbNodes()-2);
+ }
+ }
+ if(!f->IsQuadratic()) {
+ SMDS_ElemIteratorPtr nIt = f->nodesIterator();
+ while ( nIt->more() ) {
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nIt->next() );
+ if ( n == n1 ) {
+ iLinkNode[ iSide ][ 0 ] = iNode;
+ }
+ else if ( n == n2 ) {
+ iLinkNode[ iSide ][ 1 ] = iNode;
+ }
+ //else if ( notLinkNodes[ iSide ][ 0 ] )
+ // notLinkNodes[ iSide ][ 1 ] = n;
+ //else
+ // notLinkNodes[ iSide ][ 0 ] = n;
+ else {
+ nbl++;
+ if(iSide==0)
+ notLinkNodes1[nbl] = n;
+ //notLinkNodes1.push_back(n);
+ else
+ notLinkNodes2[nbl] = n;
+ //notLinkNodes2.push_back(n);
+ }
+ //faceNodes[ iSide ][ iNode++ ] = n;
+ if(iSide==0) {
+ fnodes1[iNode++] = n;
+ }
+ else {
+ fnodes2[iNode++] = n;
+ }
+ }
+ }
+ else { // f->IsQuadratic()
+ const SMDS_VtkFace* F =
+ dynamic_cast<const SMDS_VtkFace*>(f);
+ if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
+ // use special nodes iterator
+ SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
+ while ( anIter->more() ) {
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( anIter->next() );
+ if ( n == n1 ) {
+ iLinkNode[ iSide ][ 0 ] = iNode;
+ }
+ else if ( n == n2 ) {
+ iLinkNode[ iSide ][ 1 ] = iNode;
+ }
+ else {
+ nbl++;
+ if(iSide==0) {
+ notLinkNodes1[nbl] = n;
+ }
+ else {
+ notLinkNodes2[nbl] = n;
+ }
+ }
+ if(iSide==0) {
+ fnodes1[iNode++] = n;
+ }
+ else {
+ fnodes2[iNode++] = n;
+ }
+ }
+ }
+ //faceNodes[ iSide ][ iNode ] = faceNodes[ iSide ][ 0 ];
+ if(iSide==0) {
+ fnodes1[iNode] = fnodes1[0];
+ }
+ else {
+ fnodes2[iNode] = fnodes1[0];
+ }
+ }