+ // loop on nodes on seam
+ list< const SMDS_MeshNode* >::iterator noSeIt = seamNodes.begin();
+ for ( ; noSeIt != seamNodes.end(); ++noSeIt ) {
+ const SMDS_MeshNode* nSeam = *noSeIt;
+ map< const SMDS_MeshNode*, gp_XY* >::iterator n_uv = uvMap.find( nSeam );
+ if ( n_uv == uvMap.end() )
+ continue;
+ // set the first UV
+ n_uv->second->SetCoord( iPar, uv1.Coord( iPar ));
+ // set the second UV
+ listUV.push_back( *n_uv->second );
+ listUV.back().SetCoord( iPar, uv2.Coord( iPar ));
+ if ( uvMap2.empty() )
+ uvMap2 = uvMap; // copy the uvMap contents
+ uvMap2[ nSeam ] = &listUV.back();
+
+ // collect movable nodes linked to ones on seam in nodesNearSeam
+ SMDS_ElemIteratorPtr eIt = nSeam->GetInverseElementIterator(SMDSAbs_Face);
+ while ( eIt->more() ) {
+ const SMDS_MeshElement* e = eIt->next();
+ int nbUseMap1 = 0, nbUseMap2 = 0;
+ SMDS_ElemIteratorPtr nIt = e->nodesIterator();
+ int nn = 0, nbn = e->NbNodes();
+ if(e->IsQuadratic()) nbn = nbn/2;
+ while ( nn++ < nbn )
+ {
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nIt->next() );
+ if (n == nSeam ||
+ setMovableNodes.find( n ) == setMovableNodes.end() )
+ continue;
+ // add only nodes being closer to uv2 than to uv1
+ // gp_Pnt pMid (0.5 * ( n->X() + nSeam->X() ),
+ // 0.5 * ( n->Y() + nSeam->Y() ),
+ // 0.5 * ( n->Z() + nSeam->Z() ));
+ // gp_XY uv;
+ // getClosestUV( projector, pMid, uv );
+ double x = uvMap[ n ]->Coord( iPar );
+ if ( Abs( uv1.Coord( iPar ) - x ) >
+ Abs( uv2.Coord( iPar ) - x )) {
+ nodesNearSeam.insert( n );
+ nbUseMap2++;
+ }
+ else
+ nbUseMap1++;
+ }
+ // for centroidalSmooth all element nodes must
+ // be on one side of a seam
+ if ( theSmoothMethod == CENTROIDAL && nbUseMap1 && nbUseMap2 ) {
+ SMDS_ElemIteratorPtr nIt = e->nodesIterator();
+ nn = 0;
+ while ( nn++ < nbn ) {
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nIt->next() );
+ setMovableNodes.erase( n );
+ }
+ }
+ }
+ } // loop on nodes on seam
+ } // loop on edge of a face
+ } // if ( !face.IsNull() )
+
+ if ( setMovableNodes.empty() ) {
+ MESSAGE( "Face id : " << *fId << " - NO SMOOTHING: no nodes to move!!!");
+ continue; // goto next face
+ }
+
+ // -------------
+ // SMOOTHING //
+ // -------------
+
+ int it = -1;
+ double maxRatio = -1., maxDisplacement = -1.;
+ set<const SMDS_MeshNode*>::iterator nodeToMove;
+ for ( it = 0; it < theNbIterations; it++ ) {
+ maxDisplacement = 0.;
+ nodeToMove = setMovableNodes.begin();
+ for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) {
+ const SMDS_MeshNode* node = (*nodeToMove);
+ gp_XYZ aPrevPos ( node->X(), node->Y(), node->Z() );
+
+ // smooth
+ bool map2 = ( nodesNearSeam.find( node ) != nodesNearSeam.end() );
+ if ( theSmoothMethod == LAPLACIAN )
+ laplacianSmooth( node, surface, map2 ? uvMap2 : uvMap );
+ else
+ centroidalSmooth( node, surface, map2 ? uvMap2 : uvMap );
+
+ // node displacement
+ gp_XYZ aNewPos ( node->X(), node->Y(), node->Z() );
+ Standard_Real aDispl = (aPrevPos - aNewPos).SquareModulus();
+ if ( aDispl > maxDisplacement )
+ maxDisplacement = aDispl;
+ }
+ // no node movement => exit
+ //if ( maxDisplacement < 1.e-16 ) {
+ if ( maxDisplacement < disttol ) {
+ MESSAGE("-- no node movement --");
+ break;
+ }
+
+ // check elements quality
+ maxRatio = 0;
+ list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
+ for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
+ const SMDS_MeshElement* elem = (*elemIt);
+ if ( !elem || elem->GetType() != SMDSAbs_Face )
+ continue;
+ SMESH::Controls::TSequenceOfXYZ aPoints;
+ if ( aQualityFunc.GetPoints( elem, aPoints )) {
+ double aValue = aQualityFunc.GetValue( aPoints );
+ if ( aValue > maxRatio )
+ maxRatio = aValue;
+ }
+ }
+ if ( maxRatio <= theTgtAspectRatio ) {
+ MESSAGE("-- quality achived --");
+ break;
+ }
+ if (it+1 == theNbIterations) {
+ MESSAGE("-- Iteration limit exceeded --");
+ }
+ } // smoothing iterations
+
+ MESSAGE(" Face id: " << *fId <<
+ " Nb iterstions: " << it <<
+ " Displacement: " << maxDisplacement <<
+ " Aspect Ratio " << maxRatio);
+
+ // ---------------------------------------
+ // new nodes positions are computed,
+ // record movement in DS and set new UV
+ // ---------------------------------------
+ nodeToMove = setMovableNodes.begin();
+ for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) {
+ SMDS_MeshNode* node = const_cast< SMDS_MeshNode* > (*nodeToMove);
+ aMesh->MoveNode( node, node->X(), node->Y(), node->Z() );
+ map< const SMDS_MeshNode*, gp_XY* >::iterator node_uv = uvMap.find( node );
+ if ( node_uv != uvMap.end() ) {
+ gp_XY* uv = node_uv->second;
+ node->SetPosition
+ ( SMDS_PositionPtr( new SMDS_FacePosition( uv->X(), uv->Y() )));
+ }
+ }
+
+ // move medium nodes of quadratic elements
+ if ( isQuadratic )
+ {
+ vector<const SMDS_MeshNode*> nodes;
+ bool checkUV;
+ list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
+ for ( ; elemIt != elemsOnFace.end(); ++elemIt )
+ {
+ const SMDS_MeshElement* QF = *elemIt;
+ if ( QF->IsQuadratic() )
+ {
+ nodes.assign( SMDS_MeshElement::iterator( QF->interlacedNodesElemIterator() ),
+ SMDS_MeshElement::iterator() );
+ nodes.push_back( nodes[0] );
+ gp_Pnt xyz;
+ for (size_t i = 1; i < nodes.size(); i += 2 ) // i points to a medium node
+ {
+ if ( !surface.IsNull() )
+ {
+ gp_XY uv1 = helper.GetNodeUV( face, nodes[i-1], nodes[i+1], &checkUV );
+ gp_XY uv2 = helper.GetNodeUV( face, nodes[i+1], nodes[i-1], &checkUV );
+ gp_XY uv = helper.GetMiddleUV( surface, uv1, uv2 );
+ xyz = surface->Value( uv.X(), uv.Y() );
+ }
+ else {
+ xyz = 0.5 * ( SMESH_TNodeXYZ( nodes[i-1] ) + SMESH_TNodeXYZ( nodes[i+1] ));
+ }
+ if (( SMESH_TNodeXYZ( nodes[i] ) - xyz.XYZ() ).Modulus() > disttol )
+ // we have to move a medium node
+ aMesh->MoveNode( nodes[i], xyz.X(), xyz.Y(), xyz.Z() );
+ }
+ }
+ }
+ }
+
+ } // loop on face ids
+
+}
+
+namespace
+{
+ //=======================================================================
+ //function : isReverse
+ //purpose : Return true if normal of prevNodes is not co-directied with
+ // gp_Vec(prevNodes[iNotSame],nextNodes[iNotSame]).
+ // iNotSame is where prevNodes and nextNodes are different.
+ // If result is true then future volume orientation is OK
+ //=======================================================================
+
+ bool isReverse(const SMDS_MeshElement* face,
+ const vector<const SMDS_MeshNode*>& prevNodes,
+ const vector<const SMDS_MeshNode*>& nextNodes,
+ const int iNotSame)
+ {
+
+ SMESH_TNodeXYZ pP = prevNodes[ iNotSame ];
+ SMESH_TNodeXYZ pN = nextNodes[ iNotSame ];
+ gp_XYZ extrDir( pN - pP ), faceNorm;
+ SMESH_MeshAlgos::FaceNormal( face, faceNorm, /*normalized=*/false );
+
+ return faceNorm * extrDir < 0.0;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Assure that theElemSets[0] holds elements, not nodes
+ */
+ //================================================================================
+
+ void setElemsFirst( TIDSortedElemSet theElemSets[2] )
+ {
+ if ( !theElemSets[0].empty() &&
+ (*theElemSets[0].begin())->GetType() == SMDSAbs_Node )
+ {
+ std::swap( theElemSets[0], theElemSets[1] );
+ }
+ else if ( !theElemSets[1].empty() &&
+ (*theElemSets[1].begin())->GetType() != SMDSAbs_Node )
+ {
+ std::swap( theElemSets[0], theElemSets[1] );
+ }
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Create elements by sweeping an element
+ * \param elem - element to sweep
+ * \param newNodesItVec - nodes generated from each node of the element
+ * \param newElems - generated elements
+ * \param nbSteps - number of sweeping steps
+ * \param srcElements - to append elem for each generated element
+ */
+//=======================================================================
+
+void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement* elem,
+ const vector<TNodeOfNodeListMapItr> & newNodesItVec,
+ list<const SMDS_MeshElement*>& newElems,
+ const size_t nbSteps,
+ SMESH_SequenceOfElemPtr& srcElements)
+{
+ //MESSAGE("sweepElement " << nbSteps);
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ const int nbNodes = elem->NbNodes();
+ const int nbCorners = elem->NbCornerNodes();
+ SMDSAbs_EntityType baseType = elem->GetEntityType(); /* it can change in case of
+ polyhedron creation !!! */
+ // Loop on elem nodes:
+ // find new nodes and detect same nodes indices
+ vector < list< const SMDS_MeshNode* >::const_iterator > itNN( nbNodes );
+ vector<const SMDS_MeshNode*> prevNod( nbNodes );
+ vector<const SMDS_MeshNode*> nextNod( nbNodes );
+ vector<const SMDS_MeshNode*> midlNod( nbNodes );
+
+ int iNode, nbSame = 0, nbDouble = 0, iNotSameNode = 0;
+ vector<int> sames(nbNodes);
+ vector<bool> isSingleNode(nbNodes);
+
+ for ( iNode = 0; iNode < nbNodes; iNode++ ) {
+ TNodeOfNodeListMapItr nnIt = newNodesItVec[ iNode ];
+ const SMDS_MeshNode* node = nnIt->first;
+ const list< const SMDS_MeshNode* > & listNewNodes = nnIt->second;
+ if ( listNewNodes.empty() )
+ return;
+
+ itNN [ iNode ] = listNewNodes.begin();
+ prevNod[ iNode ] = node;
+ nextNod[ iNode ] = listNewNodes.front();
+
+ isSingleNode[iNode] = (listNewNodes.size()==nbSteps); /* medium node of quadratic or
+ corner node of linear */
+ if ( prevNod[ iNode ] != nextNod [ iNode ])
+ nbDouble += !isSingleNode[iNode];
+
+ if( iNode < nbCorners ) { // check corners only
+ if ( prevNod[ iNode ] == nextNod [ iNode ])
+ sames[nbSame++] = iNode;
+ else
+ iNotSameNode = iNode;
+ }
+ }
+
+ if ( nbSame == nbNodes || nbSame > 2) {
+ MESSAGE( " Too many same nodes of element " << elem->GetID() );
+ return;
+ }
+
+ if ( elem->GetType() == SMDSAbs_Face && !isReverse( elem, prevNod, nextNod, iNotSameNode ))
+ {
+ // fix nodes order to have bottom normal external
+ if ( baseType == SMDSEntity_Polygon )
+ {
+ std::reverse( itNN.begin(), itNN.end() );
+ std::reverse( prevNod.begin(), prevNod.end() );
+ std::reverse( midlNod.begin(), midlNod.end() );
+ std::reverse( nextNod.begin(), nextNod.end() );
+ std::reverse( isSingleNode.begin(), isSingleNode.end() );
+ }
+ else
+ {
+ const vector<int>& ind = SMDS_MeshCell::reverseSmdsOrder( baseType, nbNodes );
+ SMDS_MeshCell::applyInterlace( ind, itNN );
+ SMDS_MeshCell::applyInterlace( ind, prevNod );
+ SMDS_MeshCell::applyInterlace( ind, nextNod );
+ SMDS_MeshCell::applyInterlace( ind, midlNod );
+ SMDS_MeshCell::applyInterlace( ind, isSingleNode );
+ if ( nbSame > 0 )
+ {
+ sames[nbSame] = iNotSameNode;
+ for ( int j = 0; j <= nbSame; ++j )
+ for ( size_t i = 0; i < ind.size(); ++i )
+ if ( ind[i] == sames[j] )
+ {
+ sames[j] = i;
+ break;
+ }
+ iNotSameNode = sames[nbSame];
+ }
+ }
+ }
+ else if ( elem->GetType() == SMDSAbs_Edge )
+ {
+ // orient a new face same as adjacent one
+ int i1, i2;
+ const SMDS_MeshElement* e;
+ TIDSortedElemSet dummy;
+ if (( e = SMESH_MeshAlgos::FindFaceInSet( nextNod[0], prevNod[0], dummy,dummy, &i1, &i2 )) ||
+ ( e = SMESH_MeshAlgos::FindFaceInSet( prevNod[1], nextNod[1], dummy,dummy, &i1, &i2 )) ||
+ ( e = SMESH_MeshAlgos::FindFaceInSet( prevNod[0], prevNod[1], dummy,dummy, &i1, &i2 )))
+ {
+ // there is an adjacent face, check order of nodes in it
+ bool sameOrder = ( Abs( i2 - i1 ) == 1 ) ? ( i2 > i1 ) : ( i2 < i1 );
+ if ( sameOrder )
+ {
+ std::swap( itNN[0], itNN[1] );
+ std::swap( prevNod[0], prevNod[1] );
+ std::swap( nextNod[0], nextNod[1] );
+ isSingleNode.swap( isSingleNode[0], isSingleNode[1] );
+ if ( nbSame > 0 )
+ sames[0] = 1 - sames[0];
+ iNotSameNode = 1 - iNotSameNode;
+ }
+ }
+ }
+
+ int iSameNode = 0, iBeforeSame = 0, iAfterSame = 0, iOpposSame = 0;
+ if ( nbSame > 0 ) {
+ iSameNode = sames[ nbSame-1 ];
+ iBeforeSame = ( iSameNode + nbCorners - 1 ) % nbCorners;
+ iAfterSame = ( iSameNode + 1 ) % nbCorners;
+ iOpposSame = ( iSameNode - 2 < 0 ? iSameNode + 2 : iSameNode - 2 );
+ }
+
+ if ( baseType == SMDSEntity_Polygon )
+ {
+ if ( nbNodes == 3 ) baseType = SMDSEntity_Triangle;
+ else if ( nbNodes == 4 ) baseType = SMDSEntity_Quadrangle;
+ }
+ else if ( baseType == SMDSEntity_Quad_Polygon )
+ {
+ if ( nbNodes == 6 ) baseType = SMDSEntity_Quad_Triangle;
+ else if ( nbNodes == 8 ) baseType = SMDSEntity_Quad_Quadrangle;
+ }
+
+ // make new elements
+ for ( size_t iStep = 0; iStep < nbSteps; iStep++ )
+ {
+ // get next nodes
+ for ( iNode = 0; iNode < nbNodes; iNode++ )
+ {
+ midlNod[ iNode ] = isSingleNode[iNode] ? 0 : *itNN[ iNode ]++;
+ nextNod[ iNode ] = *itNN[ iNode ]++;
+ }
+
+ SMDS_MeshElement* aNewElem = 0;
+ /*if(!elem->IsPoly())*/ {
+ switch ( baseType ) {
+ case SMDSEntity_0D:
+ case SMDSEntity_Node: { // sweep NODE
+ if ( nbSame == 0 ) {
+ if ( isSingleNode[0] )
+ aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ] );
+ else
+ aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ], midlNod[ 0 ] );
+ }
+ else
+ return;
+ break;
+ }
+ case SMDSEntity_Edge: { // sweep EDGE
+ if ( nbDouble == 0 )
+ {
+ if ( nbSame == 0 ) // ---> quadrangle
+ aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
+ nextNod[ 1 ], nextNod[ 0 ] );
+ else // ---> triangle
+ aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
+ nextNod[ iNotSameNode ] );
+ }
+ else // ---> polygon
+ {
+ vector<const SMDS_MeshNode*> poly_nodes;
+ poly_nodes.push_back( prevNod[0] );
+ poly_nodes.push_back( prevNod[1] );
+ if ( prevNod[1] != nextNod[1] )
+ {
+ if ( midlNod[1]) poly_nodes.push_back( midlNod[1]);
+ poly_nodes.push_back( nextNod[1] );
+ }
+ if ( prevNod[0] != nextNod[0] )
+ {
+ poly_nodes.push_back( nextNod[0] );
+ if ( midlNod[0]) poly_nodes.push_back( midlNod[0]);
+ }
+ switch ( poly_nodes.size() ) {
+ case 3:
+ aNewElem = aMesh->AddFace( poly_nodes[ 0 ], poly_nodes[ 1 ], poly_nodes[ 2 ]);
+ break;
+ case 4:
+ aNewElem = aMesh->AddFace( poly_nodes[ 0 ], poly_nodes[ 1 ],
+ poly_nodes[ 2 ], poly_nodes[ 3 ]);
+ break;
+ default:
+ aNewElem = aMesh->AddPolygonalFace (poly_nodes);
+ }
+ }
+ break;
+ }
+ case SMDSEntity_Triangle: // TRIANGLE --->
+ {
+ if ( nbDouble > 0 ) break;
+ if ( nbSame == 0 ) // ---> pentahedron
+ aNewElem = aMesh->AddVolume (prevNod[ 0 ], prevNod[ 1 ], prevNod[ 2 ],
+ nextNod[ 0 ], nextNod[ 1 ], nextNod[ 2 ] );
+
+ else if ( nbSame == 1 ) // ---> pyramid
+ aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
+ nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
+ nextNod[ iSameNode ]);
+
+ else // 2 same nodes: ---> tetrahedron
+ aNewElem = aMesh->AddVolume (prevNod[ 0 ], prevNod[ 1 ], prevNod[ 2 ],
+ nextNod[ iNotSameNode ]);
+ break;
+ }
+ case SMDSEntity_Quad_Edge: // sweep quadratic EDGE --->
+ {
+ if ( nbSame == 2 )
+ return;
+ if ( nbDouble+nbSame == 2 )
+ {
+ if(nbSame==0) { // ---> quadratic quadrangle
+ aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[1], nextNod[0],
+ prevNod[2], midlNod[1], nextNod[2], midlNod[0]);
+ }
+ else { //(nbSame==1) // ---> quadratic triangle
+ if(sames[0]==2) {
+ return; // medium node on axis
+ }
+ else if(sames[0]==0)
+ aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[1],
+ prevNod[2], midlNod[1], nextNod[2] );
+ else // sames[0]==1
+ aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[0],
+ prevNod[2], nextNod[2], midlNod[0]);
+ }
+ }
+ else if ( nbDouble == 3 )
+ {
+ if ( nbSame == 0 ) { // ---> bi-quadratic quadrangle
+ aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[1], nextNod[0],
+ prevNod[2], midlNod[1], nextNod[2], midlNod[0], midlNod[2]);
+ }
+ }
+ else
+ return;
+ break;
+ }
+ case SMDSEntity_Quadrangle: { // sweep QUADRANGLE --->
+ if ( nbDouble > 0 ) break;
+
+ if ( nbSame == 0 ) // ---> hexahedron
+ aNewElem = aMesh->AddVolume (prevNod[ 0 ], prevNod[ 1 ], prevNod[ 2 ], prevNod[ 3 ],
+ nextNod[ 0 ], nextNod[ 1 ], nextNod[ 2 ], nextNod[ 3 ]);
+
+ else if ( nbSame == 1 ) { // ---> pyramid + pentahedron
+ aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
+ nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
+ nextNod[ iSameNode ]);
+ newElems.push_back( aNewElem );
+ aNewElem = aMesh->AddVolume (prevNod[ iAfterSame ], prevNod[ iOpposSame ],
+ prevNod[ iBeforeSame ], nextNod[ iAfterSame ],
+ nextNod[ iOpposSame ], nextNod[ iBeforeSame ] );
+ }
+ else if ( nbSame == 2 ) { // ---> pentahedron
+ if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] )
+ // iBeforeSame is same too
+ aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iOpposSame ],
+ nextNod[ iOpposSame ], prevNod[ iSameNode ],
+ prevNod[ iAfterSame ], nextNod[ iAfterSame ]);
+ else
+ // iAfterSame is same too
+ aNewElem = aMesh->AddVolume (prevNod[ iSameNode ], prevNod[ iBeforeSame ],
+ nextNod[ iBeforeSame ], prevNod[ iAfterSame ],
+ prevNod[ iOpposSame ], nextNod[ iOpposSame ]);
+ }
+ break;
+ }
+ case SMDSEntity_Quad_Triangle: // sweep (Bi)Quadratic TRIANGLE --->
+ case SMDSEntity_BiQuad_Triangle: /* ??? */ {
+ if ( nbDouble+nbSame != 3 ) break;
+ if(nbSame==0) {
+ // ---> pentahedron with 15 nodes
+ aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2],
+ nextNod[0], nextNod[1], nextNod[2],
+ prevNod[3], prevNod[4], prevNod[5],
+ nextNod[3], nextNod[4], nextNod[5],
+ midlNod[0], midlNod[1], midlNod[2]);
+ }
+ else if(nbSame==1) {
+ // ---> 2d order pyramid of 13 nodes
+ int apex = iSameNode;
+ int i0 = ( apex + 1 ) % nbCorners;
+ int i1 = ( apex - 1 + nbCorners ) % nbCorners;
+ int i0a = apex + 3;
+ int i1a = i1 + 3;
+ int i01 = i0 + 3;
+ aNewElem = aMesh->AddVolume(prevNod[i1], prevNod[i0],
+ nextNod[i0], nextNod[i1], prevNod[apex],
+ prevNod[i01], midlNod[i0],
+ nextNod[i01], midlNod[i1],
+ prevNod[i1a], prevNod[i0a],
+ nextNod[i0a], nextNod[i1a]);
+ }
+ else if(nbSame==2) {
+ // ---> 2d order tetrahedron of 10 nodes
+ int n1 = iNotSameNode;
+ int n2 = ( n1 + 1 ) % nbCorners;
+ int n3 = ( n1 + nbCorners - 1 ) % nbCorners;
+ int n12 = n1 + 3;
+ int n23 = n2 + 3;
+ int n31 = n3 + 3;
+ aNewElem = aMesh->AddVolume (prevNod[n1], prevNod[n2], prevNod[n3], nextNod[n1],
+ prevNod[n12], prevNod[n23], prevNod[n31],
+ midlNod[n1], nextNod[n12], nextNod[n31]);
+ }
+ break;
+ }
+ case SMDSEntity_Quad_Quadrangle: { // sweep Quadratic QUADRANGLE --->
+ if( nbSame == 0 ) {
+ if ( nbDouble != 4 ) break;
+ // ---> hexahedron with 20 nodes
+ aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2], prevNod[3],
+ nextNod[0], nextNod[1], nextNod[2], nextNod[3],
+ prevNod[4], prevNod[5], prevNod[6], prevNod[7],
+ nextNod[4], nextNod[5], nextNod[6], nextNod[7],
+ midlNod[0], midlNod[1], midlNod[2], midlNod[3]);
+ }
+ else if(nbSame==1) {
+ // ---> pyramid + pentahedron - can not be created since it is needed
+ // additional middle node at the center of face
+ //INFOS( " Sweep for face " << elem->GetID() << " can not be created" );
+ return;
+ }
+ else if( nbSame == 2 ) {
+ if ( nbDouble != 2 ) break;
+ // ---> 2d order Pentahedron with 15 nodes
+ int n1,n2,n4,n5;
+ if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] ) {
+ // iBeforeSame is same too
+ n1 = iBeforeSame;
+ n2 = iOpposSame;
+ n4 = iSameNode;
+ n5 = iAfterSame;
+ }
+ else {
+ // iAfterSame is same too
+ n1 = iSameNode;
+ n2 = iBeforeSame;
+ n4 = iAfterSame;
+ n5 = iOpposSame;
+ }
+ int n12 = n2 + 4;
+ int n45 = n4 + 4;
+ int n14 = n1 + 4;
+ int n25 = n5 + 4;
+ aNewElem = aMesh->AddVolume (prevNod[n1], prevNod[n2], nextNod[n2],
+ prevNod[n4], prevNod[n5], nextNod[n5],
+ prevNod[n12], midlNod[n2], nextNod[n12],
+ prevNod[n45], midlNod[n5], nextNod[n45],
+ prevNod[n14], prevNod[n25], nextNod[n25]);
+ }
+ break;
+ }
+ case SMDSEntity_BiQuad_Quadrangle: { // sweep BiQuadratic QUADRANGLE --->
+
+ if( nbSame == 0 && nbDouble == 9 ) {
+ // ---> tri-quadratic hexahedron with 27 nodes
+ aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2], prevNod[3],
+ nextNod[0], nextNod[1], nextNod[2], nextNod[3],
+ prevNod[4], prevNod[5], prevNod[6], prevNod[7],
+ nextNod[4], nextNod[5], nextNod[6], nextNod[7],
+ midlNod[0], midlNod[1], midlNod[2], midlNod[3],
+ prevNod[8], // bottom center
+ midlNod[4], midlNod[5], midlNod[6], midlNod[7],
+ nextNod[8], // top center
+ midlNod[8]);// elem center
+ }
+ else
+ {
+ return;
+ }
+ break;
+ }
+ case SMDSEntity_Polygon: { // sweep POLYGON
+
+ if ( nbNodes == 6 && nbSame == 0 && nbDouble == 0 ) {
+ // ---> hexagonal prism
+ aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2],
+ prevNod[3], prevNod[4], prevNod[5],
+ nextNod[0], nextNod[1], nextNod[2],
+ nextNod[3], nextNod[4], nextNod[5]);
+ }
+ break;
+ }
+ case SMDSEntity_Ball:
+ return;
+
+ default:
+ break;
+ } // switch ( baseType )
+ } // scope
+
+ if ( !aNewElem && elem->GetType() == SMDSAbs_Face ) // try to create a polyherdal prism
+ {
+ if ( baseType != SMDSEntity_Polygon )
+ {
+ const std::vector<int>& ind = SMDS_MeshCell::interlacedSmdsOrder(baseType,nbNodes);
+ SMDS_MeshCell::applyInterlace( ind, prevNod );
+ SMDS_MeshCell::applyInterlace( ind, nextNod );
+ SMDS_MeshCell::applyInterlace( ind, midlNod );
+ SMDS_MeshCell::applyInterlace( ind, itNN );
+ SMDS_MeshCell::applyInterlace( ind, isSingleNode );
+ baseType = SMDSEntity_Polygon; // WARNING: change baseType !!!!
+ }
+ vector<const SMDS_MeshNode*> polyedre_nodes (nbNodes*2 + 4*nbNodes);
+ vector<int> quantities (nbNodes + 2);
+ polyedre_nodes.clear();
+ quantities.clear();
+
+ // bottom of prism
+ for (int inode = 0; inode < nbNodes; inode++)
+ polyedre_nodes.push_back( prevNod[inode] );
+ quantities.push_back( nbNodes );
+
+ // top of prism
+ polyedre_nodes.push_back( nextNod[0] );
+ for (int inode = nbNodes; inode-1; --inode )
+ polyedre_nodes.push_back( nextNod[inode-1] );
+ quantities.push_back( nbNodes );
+
+ // side faces
+ // 3--6--2
+ // | |
+ // 7 5
+ // | |
+ // 0--4--1
+ const int iQuad = elem->IsQuadratic();
+ for (int iface = 0; iface < nbNodes; iface += 1+iQuad )
+ {
+ const int prevNbNodes = polyedre_nodes.size(); // to detect degenerated face
+ int inextface = (iface+1+iQuad) % nbNodes;
+ int imid = (iface+1) % nbNodes;
+ polyedre_nodes.push_back( prevNod[inextface] ); // 0
+ if ( iQuad ) polyedre_nodes.push_back( prevNod[imid] ); // 4
+ polyedre_nodes.push_back( prevNod[iface] ); // 1
+ if ( prevNod[iface] != nextNod[iface] ) // 1 != 2
+ {
+ if ( midlNod[ iface ]) polyedre_nodes.push_back( midlNod[ iface ]); // 5
+ polyedre_nodes.push_back( nextNod[iface] ); // 2
+ }
+ if ( iQuad ) polyedre_nodes.push_back( nextNod[imid] ); // 6
+ if ( prevNod[inextface] != nextNod[inextface] ) // 0 != 3
+ {
+ polyedre_nodes.push_back( nextNod[inextface] ); // 3
+ if ( midlNod[ inextface ]) polyedre_nodes.push_back( midlNod[ inextface ]);// 7
+ }
+ const int nbFaceNodes = polyedre_nodes.size() - prevNbNodes;
+ if ( nbFaceNodes > 2 )
+ quantities.push_back( nbFaceNodes );
+ else // degenerated face
+ polyedre_nodes.resize( prevNbNodes );
+ }
+ aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
+
+ } // try to create a polyherdal prism
+
+ if ( aNewElem ) {
+ newElems.push_back( aNewElem );
+ myLastCreatedElems.Append(aNewElem);
+ srcElements.Append( elem );
+ }
+
+ // set new prev nodes
+ for ( iNode = 0; iNode < nbNodes; iNode++ )
+ prevNod[ iNode ] = nextNod[ iNode ];
+
+ } // loop on steps
+}
+
+//=======================================================================
+/*!
+ * \brief Create 1D and 2D elements around swept elements
+ * \param mapNewNodes - source nodes and ones generated from them
+ * \param newElemsMap - source elements and ones generated from them
+ * \param elemNewNodesMap - nodes generated from each node of each element
+ * \param elemSet - all swept elements
+ * \param nbSteps - number of sweeping steps
+ * \param srcElements - to append elem for each generated element
+ */
+//=======================================================================
+
+void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes,
+ TTElemOfElemListMap & newElemsMap,
+ TElemOfVecOfNnlmiMap & elemNewNodesMap,
+ TIDSortedElemSet& elemSet,
+ const int nbSteps,
+ SMESH_SequenceOfElemPtr& srcElements)
+{
+ ASSERT( newElemsMap.size() == elemNewNodesMap.size() );
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ // Find nodes belonging to only one initial element - sweep them into edges.
+
+ TNodeOfNodeListMapItr nList = mapNewNodes.begin();
+ for ( ; nList != mapNewNodes.end(); nList++ )
+ {
+ const SMDS_MeshNode* node =
+ static_cast<const SMDS_MeshNode*>( nList->first );
+ if ( newElemsMap.count( node ))
+ continue; // node was extruded into edge
+ SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
+ int nbInitElems = 0;
+ const SMDS_MeshElement* el = 0;
+ SMDSAbs_ElementType highType = SMDSAbs_Edge; // count most complex elements only
+ while ( eIt->more() && nbInitElems < 2 ) {
+ const SMDS_MeshElement* e = eIt->next();
+ SMDSAbs_ElementType type = e->GetType();
+ if ( type == SMDSAbs_Volume ||
+ type < highType ||
+ !elemSet.count(e))
+ continue;
+ if ( type > highType ) {
+ nbInitElems = 0;
+ highType = type;
+ }
+ el = e;
+ ++nbInitElems;
+ }
+ if ( nbInitElems == 1 ) {
+ bool NotCreateEdge = el && el->IsMediumNode(node);
+ if(!NotCreateEdge) {
+ vector<TNodeOfNodeListMapItr> newNodesItVec( 1, nList );
+ list<const SMDS_MeshElement*> newEdges;
+ sweepElement( node, newNodesItVec, newEdges, nbSteps, srcElements );
+ }
+ }
+ }
+
+ // Make a ceiling for each element ie an equal element of last new nodes.
+ // Find free links of faces - make edges and sweep them into faces.
+
+ ElemFeatures polyFace( SMDSAbs_Face, /*isPoly=*/true ), anyFace;
+
+ TTElemOfElemListMap::iterator itElem = newElemsMap.begin();
+ TElemOfVecOfNnlmiMap::iterator itElemNodes = elemNewNodesMap.begin();
+ for ( ; itElem != newElemsMap.end(); itElem++, itElemNodes++ )
+ {
+ const SMDS_MeshElement* elem = itElem->first;
+ vector<TNodeOfNodeListMapItr>& vecNewNodes = itElemNodes->second;
+
+ if(itElem->second.size()==0) continue;
+
+ const bool isQuadratic = elem->IsQuadratic();
+
+ if ( elem->GetType() == SMDSAbs_Edge ) {
+ // create a ceiling edge
+ if ( !isQuadratic ) {
+ if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
+ vecNewNodes[ 1 ]->second.back())) {
+ myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
+ vecNewNodes[ 1 ]->second.back()));
+ srcElements.Append( elem );
+ }
+ }
+ else {
+ if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
+ vecNewNodes[ 1 ]->second.back(),
+ vecNewNodes[ 2 ]->second.back())) {
+ myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
+ vecNewNodes[ 1 ]->second.back(),
+ vecNewNodes[ 2 ]->second.back()));
+ srcElements.Append( elem );
+ }
+ }
+ }
+ if ( elem->GetType() != SMDSAbs_Face )
+ continue;
+
+ bool hasFreeLinks = false;
+
+ TIDSortedElemSet avoidSet;
+ avoidSet.insert( elem );
+
+ set<const SMDS_MeshNode*> aFaceLastNodes;
+ int iNode, nbNodes = vecNewNodes.size();
+ if ( !isQuadratic ) {
+ // loop on the face nodes
+ for ( iNode = 0; iNode < nbNodes; iNode++ ) {
+ aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
+ // look for free links of the face
+ int iNext = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1;
+ const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
+ const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
+ // check if a link n1-n2 is free
+ if ( ! SMESH_MeshAlgos::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
+ hasFreeLinks = true;
+ // make a new edge and a ceiling for a new edge
+ const SMDS_MeshElement* edge;
+ if ( ! ( edge = aMesh->FindEdge( n1, n2 ))) {
+ myLastCreatedElems.Append( edge = aMesh->AddEdge( n1, n2 )); // free link edge
+ srcElements.Append( myLastCreatedElems.Last() );
+ }
+ n1 = vecNewNodes[ iNode ]->second.back();
+ n2 = vecNewNodes[ iNext ]->second.back();
+ if ( !aMesh->FindEdge( n1, n2 )) {
+ myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // new edge ceiling
+ srcElements.Append( edge );
+ }
+ }
+ }
+ }
+ else { // elem is quadratic face
+ int nbn = nbNodes/2;
+ for ( iNode = 0; iNode < nbn; iNode++ ) {
+ aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
+ int iNext = ( iNode + 1 == nbn ) ? 0 : iNode + 1;
+ const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
+ const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
+ const SMDS_MeshNode* n3 = vecNewNodes[ iNode+nbn ]->first;
+ // check if a link is free
+ if ( ! SMESH_MeshAlgos::FindFaceInSet ( n1, n2, elemSet, avoidSet ) &&
+ ! SMESH_MeshAlgos::FindFaceInSet ( n1, n3, elemSet, avoidSet ) &&
+ ! SMESH_MeshAlgos::FindFaceInSet ( n3, n2, elemSet, avoidSet ) ) {
+ hasFreeLinks = true;
+ // make an edge and a ceiling for a new edge
+ // find medium node
+ if ( !aMesh->FindEdge( n1, n2, n3 )) {
+ myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // free link edge
+ srcElements.Append( elem );
+ }
+ n1 = vecNewNodes[ iNode ]->second.back();
+ n2 = vecNewNodes[ iNext ]->second.back();
+ n3 = vecNewNodes[ iNode+nbn ]->second.back();
+ if ( !aMesh->FindEdge( n1, n2, n3 )) {
+ myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // ceiling edge
+ srcElements.Append( elem );
+ }
+ }
+ }
+ for ( iNode = nbn; iNode < nbNodes; iNode++ ) {
+ aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
+ }
+ }
+
+ // sweep free links into faces
+
+ if ( hasFreeLinks ) {
+ list<const SMDS_MeshElement*> & newVolumes = itElem->second;
+ int iVol, volNb, nbVolumesByStep = newVolumes.size() / nbSteps;
+
+ set<const SMDS_MeshNode*> initNodeSet, topNodeSet, faceNodeSet;
+ set<const SMDS_MeshNode*> initNodeSetNoCenter/*, topNodeSetNoCenter*/;
+ for ( iNode = 0; iNode < nbNodes; iNode++ ) {
+ initNodeSet.insert( vecNewNodes[ iNode ]->first );
+ topNodeSet .insert( vecNewNodes[ iNode ]->second.back() );
+ }
+ if ( isQuadratic && nbNodes % 2 ) { // node set for the case of a biquadratic
+ initNodeSetNoCenter = initNodeSet; // swept face and a not biquadratic volume
+ initNodeSetNoCenter.erase( vecNewNodes.back()->first );
+ }
+ for ( volNb = 0; volNb < nbVolumesByStep; volNb++ ) {
+ list<const SMDS_MeshElement*>::iterator v = newVolumes.begin();
+ std::advance( v, volNb );
+ // find indices of free faces of a volume and their source edges
+ list< int > freeInd;
+ list< const SMDS_MeshElement* > srcEdges; // source edges of free faces
+ SMDS_VolumeTool vTool( *v, /*ignoreCentralNodes=*/false );
+ int iF, nbF = vTool.NbFaces();
+ for ( iF = 0; iF < nbF; iF ++ ) {
+ if (vTool.IsFreeFace( iF ) &&
+ vTool.GetFaceNodes( iF, faceNodeSet ) &&
+ initNodeSet != faceNodeSet) // except an initial face
+ {
+ if ( nbSteps == 1 && faceNodeSet == topNodeSet )
+ continue;
+ if ( faceNodeSet == initNodeSetNoCenter )
+ continue;
+ freeInd.push_back( iF );
+ // find source edge of a free face iF
+ vector<const SMDS_MeshNode*> commonNodes; // shared by the initial and free faces
+ vector<const SMDS_MeshNode*>::iterator lastCommom;
+ commonNodes.resize( nbNodes, 0 );
+ lastCommom = std::set_intersection( faceNodeSet.begin(), faceNodeSet.end(),
+ initNodeSet.begin(), initNodeSet.end(),
+ commonNodes.begin());
+ if ( std::distance( commonNodes.begin(), lastCommom ) == 3 )
+ srcEdges.push_back(aMesh->FindEdge (commonNodes[0],commonNodes[1],commonNodes[2]));
+ else
+ srcEdges.push_back(aMesh->FindEdge (commonNodes[0],commonNodes[1]));
+#ifdef _DEBUG_
+ if ( !srcEdges.back() )
+ {
+ cout << "SMESH_MeshEditor::makeWalls(), no source edge found for a free face #"
+ << iF << " of volume #" << vTool.ID() << endl;
+ }
+#endif
+ }
+ }
+ if ( freeInd.empty() )
+ continue;
+
+ // create wall faces for all steps;
+ // if such a face has been already created by sweep of edge,
+ // assure that its orientation is OK
+ for ( int iStep = 0; iStep < nbSteps; iStep++ )
+ {
+ vTool.Set( *v, /*ignoreCentralNodes=*/false );
+ vTool.SetExternalNormal();
+ const int nextShift = vTool.IsForward() ? +1 : -1;
+ list< int >::iterator ind = freeInd.begin();
+ list< const SMDS_MeshElement* >::iterator srcEdge = srcEdges.begin();
+ for ( ; ind != freeInd.end(); ++ind, ++srcEdge ) // loop on free faces
+ {
+ const SMDS_MeshNode** nodes = vTool.GetFaceNodes( *ind );
+ int nbn = vTool.NbFaceNodes( *ind );
+ const SMDS_MeshElement * f = 0;
+ if ( nbn == 3 ) ///// triangle
+ {
+ f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]);
+ if ( !f ||
+ nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ]) + nextShift ))
+ {
+ const SMDS_MeshNode* newOrder[3] = { nodes[ 1 - nextShift ],
+ nodes[ 1 ],
+ nodes[ 1 + nextShift ] };
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace( newOrder[ 0 ], newOrder[ 1 ],
+ newOrder[ 2 ] ));
+ }
+ }
+ else if ( nbn == 4 ) ///// quadrangle
+ {
+ f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]);
+ if ( !f ||
+ nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ]) + nextShift ))
+ {
+ const SMDS_MeshNode* newOrder[4] = { nodes[ 0 ], nodes[ 2-nextShift ],
+ nodes[ 2 ], nodes[ 2+nextShift ] };
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace( newOrder[ 0 ], newOrder[ 1 ],
+ newOrder[ 2 ], newOrder[ 3 ]));
+ }
+ }
+ else if ( nbn == 6 && isQuadratic ) /////// quadratic triangle
+ {
+ f = aMesh->FindFace( nodes[0], nodes[2], nodes[4], nodes[1], nodes[3], nodes[5] );
+ if ( !f ||
+ nodes[2] != f->GetNodeWrap( f->GetNodeIndex( nodes[0] ) + 2*nextShift ))
+ {
+ const SMDS_MeshNode* newOrder[6] = { nodes[2 - 2*nextShift],
+ nodes[2],
+ nodes[2 + 2*nextShift],
+ nodes[3 - 2*nextShift],
+ nodes[3],
+ nodes[3 + 2*nextShift]};
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace( newOrder[ 0 ],
+ newOrder[ 1 ],
+ newOrder[ 2 ],
+ newOrder[ 3 ],
+ newOrder[ 4 ],
+ newOrder[ 5 ] ));
+ }
+ }
+ else if ( nbn == 8 && isQuadratic ) /////// quadratic quadrangle
+ {
+ f = aMesh->FindFace( nodes[0], nodes[2], nodes[4], nodes[6],
+ nodes[1], nodes[3], nodes[5], nodes[7] );
+ if ( !f ||
+ nodes[ 2 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 2*nextShift ))
+ {
+ const SMDS_MeshNode* newOrder[8] = { nodes[0],
+ nodes[4 - 2*nextShift],
+ nodes[4],
+ nodes[4 + 2*nextShift],
+ nodes[1],
+ nodes[5 - 2*nextShift],
+ nodes[5],
+ nodes[5 + 2*nextShift] };
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace(newOrder[ 0 ], newOrder[ 1 ],
+ newOrder[ 2 ], newOrder[ 3 ],
+ newOrder[ 4 ], newOrder[ 5 ],
+ newOrder[ 6 ], newOrder[ 7 ]));
+ }
+ }
+ else if ( nbn == 9 && isQuadratic ) /////// bi-quadratic quadrangle
+ {
+ f = aMesh->FindElement( vector<const SMDS_MeshNode*>( nodes, nodes+nbn ),
+ SMDSAbs_Face, /*noMedium=*/false);
+ if ( !f ||
+ nodes[ 2 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 2*nextShift ))
+ {
+ const SMDS_MeshNode* newOrder[9] = { nodes[0],
+ nodes[4 - 2*nextShift],
+ nodes[4],
+ nodes[4 + 2*nextShift],
+ nodes[1],
+ nodes[5 - 2*nextShift],
+ nodes[5],
+ nodes[5 + 2*nextShift],
+ nodes[8] };
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace(newOrder[ 0 ], newOrder[ 1 ],
+ newOrder[ 2 ], newOrder[ 3 ],
+ newOrder[ 4 ], newOrder[ 5 ],
+ newOrder[ 6 ], newOrder[ 7 ],
+ newOrder[ 8 ]));
+ }
+ }
+ else //////// polygon
+ {
+ vector<const SMDS_MeshNode*> polygon_nodes ( nodes, nodes+nbn );
+ const SMDS_MeshFace * f = aMesh->FindFace( polygon_nodes );
+ if ( !f ||
+ nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + nextShift ))
+ {
+ if ( !vTool.IsForward() )
+ std::reverse( polygon_nodes.begin(), polygon_nodes.end());
+ if ( f )
+ aMesh->ChangeElementNodes( f, &polygon_nodes[0], nbn );
+ else
+ AddElement( polygon_nodes, polyFace.SetQuad( (*v)->IsQuadratic() ));
+ }
+ }
+
+ while ( srcElements.Length() < myLastCreatedElems.Length() )
+ srcElements.Append( *srcEdge );
+
+ } // loop on free faces
+
+ // go to the next volume
+ iVol = 0;
+ while ( iVol++ < nbVolumesByStep ) v++;
+
+ } // loop on steps
+ } // loop on volumes of one step
+ } // sweep free links into faces
+
+ // Make a ceiling face with a normal external to a volume
+
+ // use SMDS_VolumeTool to get a correctly ordered nodes of a ceiling face
+ SMDS_VolumeTool lastVol( itElem->second.back(), /*ignoreCentralNodes=*/false );
+ int iF = lastVol.GetFaceIndex( aFaceLastNodes );
+
+ if ( iF < 0 && isQuadratic && nbNodes % 2 ) { // remove a central node of biquadratic
+ aFaceLastNodes.erase( vecNewNodes.back()->second.back() );
+ iF = lastVol.GetFaceIndex( aFaceLastNodes );
+ }
+ if ( iF >= 0 )
+ {
+ lastVol.SetExternalNormal();
+ const SMDS_MeshNode** nodes = lastVol.GetFaceNodes( iF );
+ const int nbn = lastVol.NbFaceNodes( iF );
+ vector<const SMDS_MeshNode*> nodeVec( nodes, nodes+nbn );
+ if ( !hasFreeLinks ||
+ !aMesh->FindElement( nodeVec, SMDSAbs_Face, /*noMedium=*/false) )
+ {
+ const vector<int>& interlace =
+ SMDS_MeshCell::interlacedSmdsOrder( elem->GetEntityType(), nbn );
+ SMDS_MeshCell::applyInterlaceRev( interlace, nodeVec );
+
+ AddElement( nodeVec, anyFace.Init( elem ));
+
+ while ( srcElements.Length() < myLastCreatedElems.Length() )
+ srcElements.Append( elem );
+ }
+ }
+ } // loop on swept elements
+}
+
+//=======================================================================
+//function : RotationSweep
+//purpose :
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::RotationSweep(TIDSortedElemSet theElemSets[2],
+ const gp_Ax1& theAxis,
+ const double theAngle,
+ const int theNbSteps,
+ const double theTol,
+ const bool theMakeGroups,
+ const bool theMakeWalls)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+
+ MESSAGE( "RotationSweep()");
+ gp_Trsf aTrsf;
+ aTrsf.SetRotation( theAxis, theAngle );
+ gp_Trsf aTrsf2;
+ aTrsf2.SetRotation( theAxis, theAngle/2. );
+
+ gp_Lin aLine( theAxis );
+ double aSqTol = theTol * theTol;
+
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ TNodeOfNodeListMap mapNewNodes;
+ TElemOfVecOfNnlmiMap mapElemNewNodes;
+ TTElemOfElemListMap newElemsMap;
+
+ const bool isQuadraticMesh = bool( myMesh->NbEdges(ORDER_QUADRATIC) +
+ myMesh->NbFaces(ORDER_QUADRATIC) +
+ myMesh->NbVolumes(ORDER_QUADRATIC) );
+ // loop on theElemSets
+ setElemsFirst( theElemSets );
+ TIDSortedElemSet::iterator itElem;
+ for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+ {
+ TIDSortedElemSet& theElems = theElemSets[ is2ndSet ];
+ for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
+ const SMDS_MeshElement* elem = *itElem;
+ if ( !elem || elem->GetType() == SMDSAbs_Volume )
+ continue;
+ vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
+ newNodesItVec.reserve( elem->NbNodes() );
+
+ // loop on elem nodes
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() )
+ {
+ const SMDS_MeshNode* node = cast2Node( itN->next() );
+
+ gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
+ double coord[3];
+ aXYZ.Coord( coord[0], coord[1], coord[2] );
+ bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
+
+ // check if a node has been already sweeped
+ TNodeOfNodeListMapItr nIt =
+ mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
+ list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
+ if ( listNewNodes.empty() )
+ {
+ // check if we are to create medium nodes between corner ones
+ bool needMediumNodes = false;
+ if ( isQuadraticMesh )
+ {
+ SMDS_ElemIteratorPtr it = node->GetInverseElementIterator();
+ while (it->more() && !needMediumNodes )
+ {
+ const SMDS_MeshElement* invElem = it->next();
+ if ( invElem != elem && !theElems.count( invElem )) continue;
+ needMediumNodes = ( invElem->IsQuadratic() && !invElem->IsMediumNode(node) );
+ if ( !needMediumNodes && invElem->GetEntityType() == SMDSEntity_BiQuad_Quadrangle )
+ needMediumNodes = true;
+ }
+ }
+
+ // make new nodes
+ const SMDS_MeshNode * newNode = node;
+ for ( int i = 0; i < theNbSteps; i++ ) {
+ if ( !isOnAxis ) {
+ if ( needMediumNodes ) // create a medium node
+ {
+ aTrsf2.Transforms( coord[0], coord[1], coord[2] );
+ newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+ myLastCreatedNodes.Append(newNode);
+ srcNodes.Append( node );
+ listNewNodes.push_back( newNode );
+ aTrsf2.Transforms( coord[0], coord[1], coord[2] );
+ }
+ else {
+ aTrsf.Transforms( coord[0], coord[1], coord[2] );
+ }
+ // create a corner node
+ newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+ myLastCreatedNodes.Append(newNode);
+ srcNodes.Append( node );
+ listNewNodes.push_back( newNode );
+ }
+ else {
+ listNewNodes.push_back( newNode );
+ // if ( needMediumNodes )
+ // listNewNodes.push_back( newNode );
+ }
+ }
+ }
+ newNodesItVec.push_back( nIt );
+ }
+ // make new elements
+ sweepElement( elem, newNodesItVec, newElemsMap[elem], theNbSteps, srcElems );
+ }
+ }
+
+ if ( theMakeWalls )
+ makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElemSets[0], theNbSteps, srcElems );
+
+ PGroupIDs newGroupIDs;
+ if ( theMakeGroups )
+ newGroupIDs = generateGroups( srcNodes, srcElems, "rotated");
+
+ return newGroupIDs;
+}
+
+//=======================================================================
+//function : ExtrusParam
+//purpose : standard construction
+//=======================================================================
+
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Vec& theStep,
+ const int theNbSteps,
+ const int theFlags,
+ const double theTolerance):
+ myDir( theStep ),
+ myFlags( theFlags ),
+ myTolerance( theTolerance ),
+ myElemsToUse( NULL )
+{
+ mySteps = new TColStd_HSequenceOfReal;
+ const double stepSize = theStep.Magnitude();
+ for (int i=1; i<=theNbSteps; i++ )
+ mySteps->Append( stepSize );
+
+ if (( theFlags & EXTRUSION_FLAG_SEW ) &&
+ ( theTolerance > 0 ))
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDirAndSew;
+ }
+ else
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDir;
+ }
+}
+
+//=======================================================================
+//function : ExtrusParam
+//purpose : steps are given explicitly
+//=======================================================================
+
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Dir& theDir,
+ Handle(TColStd_HSequenceOfReal) theSteps,
+ const int theFlags,
+ const double theTolerance):
+ myDir( theDir ),
+ mySteps( theSteps ),
+ myFlags( theFlags ),
+ myTolerance( theTolerance ),
+ myElemsToUse( NULL )
+{
+ if (( theFlags & EXTRUSION_FLAG_SEW ) &&
+ ( theTolerance > 0 ))
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDirAndSew;
+ }
+ else
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDir;
+ }
+}
+
+//=======================================================================
+//function : ExtrusParam
+//purpose : for extrusion by normal
+//=======================================================================
+
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const double theStepSize,
+ const int theNbSteps,
+ const int theFlags,
+ const int theDim ):
+ myDir( 1,0,0 ),
+ mySteps( new TColStd_HSequenceOfReal ),
+ myFlags( theFlags ),
+ myTolerance( 0 ),
+ myElemsToUse( NULL )
+{
+ for (int i = 0; i < theNbSteps; i++ )
+ mySteps->Append( theStepSize );
+
+ if ( theDim == 1 )
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByNormal1D;
+ }
+ else
+ {
+ myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByNormal2D;
+ }
+}
+
+//=======================================================================
+//function : ExtrusParam::SetElementsToUse
+//purpose : stores elements to use for extrusion by normal, depending on
+// state of EXTRUSION_FLAG_USE_INPUT_ELEMS_ONLY flag
+//=======================================================================
+
+void SMESH_MeshEditor::ExtrusParam::SetElementsToUse( const TIDSortedElemSet& elems )
+{
+ myElemsToUse = ToUseInpElemsOnly() ? & elems : 0;
+}
+
+//=======================================================================
+//function : ExtrusParam::beginStepIter
+//purpose : prepare iteration on steps
+//=======================================================================
+
+void SMESH_MeshEditor::ExtrusParam::beginStepIter( bool withMediumNodes )
+{
+ myWithMediumNodes = withMediumNodes;
+ myNextStep = 1;
+ myCurSteps.clear();
+}
+//=======================================================================
+//function : ExtrusParam::moreSteps
+//purpose : are there more steps?
+//=======================================================================
+
+bool SMESH_MeshEditor::ExtrusParam::moreSteps()
+{
+ return myNextStep <= mySteps->Length() || !myCurSteps.empty();
+}
+//=======================================================================
+//function : ExtrusParam::nextStep
+//purpose : returns the next step
+//=======================================================================
+
+double SMESH_MeshEditor::ExtrusParam::nextStep()
+{
+ double res = 0;
+ if ( !myCurSteps.empty() )
+ {
+ res = myCurSteps.back();
+ myCurSteps.pop_back();
+ }
+ else if ( myNextStep <= mySteps->Length() )
+ {
+ myCurSteps.push_back( mySteps->Value( myNextStep ));
+ ++myNextStep;
+ if ( myWithMediumNodes )
+ {
+ myCurSteps.back() /= 2.;
+ myCurSteps.push_back( myCurSteps.back() );
+ }
+ res = nextStep();
+ }
+ return res;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesByDir
+//purpose : create nodes for standard extrusion
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByDir( SMESHDS_Mesh* mesh,
+ const SMDS_MeshNode* srcNode,
+ std::list<const SMDS_MeshNode*> & newNodes,
+ const bool makeMediumNodes)
+{
+ gp_XYZ p = SMESH_TNodeXYZ( srcNode );
+
+ int nbNodes = 0;
+ for ( beginStepIter( makeMediumNodes ); moreSteps(); ++nbNodes ) // loop on steps
+ {
+ p += myDir.XYZ() * nextStep();
+ const SMDS_MeshNode * newNode = mesh->AddNode( p.X(), p.Y(), p.Z() );
+ newNodes.push_back( newNode );
+ }
+ return nbNodes;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesByDirAndSew
+//purpose : create nodes for standard extrusion with sewing
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByDirAndSew( SMESHDS_Mesh* mesh,
+ const SMDS_MeshNode* srcNode,
+ std::list<const SMDS_MeshNode*> & newNodes,
+ const bool makeMediumNodes)
+{
+ gp_XYZ P1 = SMESH_TNodeXYZ( srcNode );
+
+ int nbNodes = 0;
+ for ( beginStepIter( makeMediumNodes ); moreSteps(); ++nbNodes ) // loop on steps
+ {
+ P1 += myDir.XYZ() * nextStep();
+
+ // try to search in sequence of existing nodes
+ // if myNodes.Length()>0 we 'nave to use given sequence
+ // else - use all nodes of mesh
+ const SMDS_MeshNode * node = 0;
+ if ( myNodes.Length() > 0 ) {
+ int i;
+ for(i=1; i<=myNodes.Length(); i++) {
+ gp_XYZ P2 = SMESH_TNodeXYZ( myNodes.Value(i) );
+ if (( P1 - P2 ).SquareModulus() < myTolerance * myTolerance )
+ {
+ node = myNodes.Value(i);
+ break;
+ }
+ }
+ }
+ else {
+ SMDS_NodeIteratorPtr itn = mesh->nodesIterator();
+ while(itn->more()) {
+ SMESH_TNodeXYZ P2( itn->next() );
+ if (( P1 - P2 ).SquareModulus() < myTolerance * myTolerance )
+ {
+ node = P2._node;
+ break;
+ }
+ }
+ }
+
+ if ( !node )
+ node = mesh->AddNode( P1.X(), P1.Y(), P1.Z() );
+
+ newNodes.push_back( node );
+
+ } // loop on steps
+
+ return nbNodes;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesByNormal2D
+//purpose : create nodes for extrusion using normals of faces
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByNormal2D( SMESHDS_Mesh* mesh,
+ const SMDS_MeshNode* srcNode,
+ std::list<const SMDS_MeshNode*> & newNodes,
+ const bool makeMediumNodes)
+{
+ const bool alongAvgNorm = ( myFlags & EXTRUSION_FLAG_BY_AVG_NORMAL );
+
+ gp_XYZ p = SMESH_TNodeXYZ( srcNode );
+
+ // get normals to faces sharing srcNode
+ vector< gp_XYZ > norms, baryCenters;
+ gp_XYZ norm, avgNorm( 0,0,0 );
+ SMDS_ElemIteratorPtr faceIt = srcNode->GetInverseElementIterator( SMDSAbs_Face );
+ while ( faceIt->more() )
+ {
+ const SMDS_MeshElement* face = faceIt->next();
+ if ( myElemsToUse && !myElemsToUse->count( face ))
+ continue;
+ if ( SMESH_MeshAlgos::FaceNormal( face, norm, /*normalized=*/true ))
+ {
+ norms.push_back( norm );
+ avgNorm += norm;
+ if ( !alongAvgNorm )
+ {
+ gp_XYZ bc(0,0,0);
+ int nbN = 0;
+ for ( SMDS_ElemIteratorPtr nIt = face->nodesIterator(); nIt->more(); ++nbN )
+ bc += SMESH_TNodeXYZ( nIt->next() );
+ baryCenters.push_back( bc / nbN );
+ }
+ }
+ }
+
+ if ( norms.empty() ) return 0;
+
+ double normSize = avgNorm.Modulus();
+ if ( normSize < std::numeric_limits<double>::min() )
+ return 0;
+
+ if ( myFlags & EXTRUSION_FLAG_BY_AVG_NORMAL ) // extrude along avgNorm
+ {
+ myDir = avgNorm;
+ return makeNodesByDir( mesh, srcNode, newNodes, makeMediumNodes );
+ }
+
+ avgNorm /= normSize;
+
+ int nbNodes = 0;
+ for ( beginStepIter( makeMediumNodes ); moreSteps(); ++nbNodes ) // loop on steps
+ {
+ gp_XYZ pNew = p;
+ double stepSize = nextStep();
+
+ if ( norms.size() > 1 )
+ {
+ for ( size_t iF = 0; iF < norms.size(); ++iF ) // loop on faces
+ {
+ // translate plane of a face
+ baryCenters[ iF ] += norms[ iF ] * stepSize;
+
+ // find point of intersection of the face plane located at baryCenters[ iF ]
+ // and avgNorm located at pNew
+ double d = -( norms[ iF ] * baryCenters[ iF ]); // d of plane equation ax+by+cz+d=0
+ double dot = ( norms[ iF ] * avgNorm );
+ if ( dot < std::numeric_limits<double>::min() )
+ dot = stepSize * 1e-3;
+ double step = -( norms[ iF ] * pNew + d ) / dot;
+ pNew += step * avgNorm;
+ }
+ }
+ else
+ {
+ pNew += stepSize * avgNorm;
+ }
+ p = pNew;
+
+ const SMDS_MeshNode * newNode = mesh->AddNode( p.X(), p.Y(), p.Z() );
+ newNodes.push_back( newNode );
+ }
+ return nbNodes;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesByNormal1D
+//purpose : create nodes for extrusion using normals of edges
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByNormal1D( SMESHDS_Mesh* mesh,
+ const SMDS_MeshNode* srcNode,
+ std::list<const SMDS_MeshNode*> & newNodes,
+ const bool makeMediumNodes)
+{
+ throw SALOME_Exception("Extrusion 1D by Normal not implemented");
+ return 0;
+}
+
+//=======================================================================
+//function : ExtrusionSweep
+//purpose :
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet theElems[2],
+ const gp_Vec& theStep,
+ const int theNbSteps,
+ TTElemOfElemListMap& newElemsMap,
+ const int theFlags,
+ const double theTolerance)
+{
+ ExtrusParam aParams( theStep, theNbSteps, theFlags, theTolerance );
+ return ExtrusionSweep( theElems, aParams, newElemsMap );
+}
+
+
+//=======================================================================
+//function : ExtrusionSweep
+//purpose :
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet theElemSets[2],
+ ExtrusParam& theParams,
+ TTElemOfElemListMap& newElemsMap)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+
+ //SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ setElemsFirst( theElemSets );
+ const int nbSteps = theParams.NbSteps();
+ theParams.SetElementsToUse( theElemSets[0] );
+
+ TNodeOfNodeListMap mapNewNodes;
+ //TNodeOfNodeVecMap mapNewNodes;
+ TElemOfVecOfNnlmiMap mapElemNewNodes;
+ //TElemOfVecOfMapNodesMap mapElemNewNodes;
+
+ const bool isQuadraticMesh = bool( myMesh->NbEdges(ORDER_QUADRATIC) +
+ myMesh->NbFaces(ORDER_QUADRATIC) +
+ myMesh->NbVolumes(ORDER_QUADRATIC) );
+ // loop on theElems
+ TIDSortedElemSet::iterator itElem;
+ for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+ {
+ TIDSortedElemSet& theElems = theElemSets[ is2ndSet ];
+ for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
+ {
+ // check element type
+ const SMDS_MeshElement* elem = *itElem;
+ if ( !elem || elem->GetType() == SMDSAbs_Volume )
+ continue;
+
+ const size_t nbNodes = elem->NbNodes();
+ vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
+ newNodesItVec.reserve( nbNodes );
+
+ // loop on elem nodes
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() )
+ {
+ // check if a node has been already sweeped
+ const SMDS_MeshNode* node = cast2Node( itN->next() );
+ TNodeOfNodeListMap::iterator nIt =
+ mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
+ list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
+ if ( listNewNodes.empty() )
+ {
+ // make new nodes
+
+ // check if we are to create medium nodes between corner ones
+ bool needMediumNodes = false;
+ if ( isQuadraticMesh )
+ {
+ SMDS_ElemIteratorPtr it = node->GetInverseElementIterator();
+ while (it->more() && !needMediumNodes )
+ {
+ const SMDS_MeshElement* invElem = it->next();
+ if ( invElem != elem && !theElems.count( invElem )) continue;
+ needMediumNodes = ( invElem->IsQuadratic() && !invElem->IsMediumNode(node) );
+ if ( !needMediumNodes && invElem->GetEntityType() == SMDSEntity_BiQuad_Quadrangle )
+ needMediumNodes = true;
+ }
+ }
+ // create nodes for all steps
+ if ( theParams.MakeNodes( GetMeshDS(), node, listNewNodes, needMediumNodes ))
+ {
+ list<const SMDS_MeshNode*>::iterator newNodesIt = listNewNodes.begin();
+ for ( ; newNodesIt != listNewNodes.end(); ++newNodesIt )
+ {
+ myLastCreatedNodes.Append( *newNodesIt );
+ srcNodes.Append( node );
+ }
+ }
+ else
+ {
+ break; // newNodesItVec will be shorter than nbNodes
+ }
+ }
+ newNodesItVec.push_back( nIt );
+ }
+ // make new elements
+ if ( newNodesItVec.size() == nbNodes )
+ sweepElement( elem, newNodesItVec, newElemsMap[elem], nbSteps, srcElems );
+ }
+ }
+
+ if ( theParams.ToMakeBoundary() ) {
+ makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElemSets[0], nbSteps, srcElems );
+ }
+ PGroupIDs newGroupIDs;
+ if ( theParams.ToMakeGroups() )
+ newGroupIDs = generateGroups( srcNodes, srcElems, "extruded");
+
+ return newGroupIDs;
+}
+
+//=======================================================================
+//function : ExtrusionAlongTrack
+//purpose :
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet theElements[2],
+ SMESH_subMesh* theTrack,
+ const SMDS_MeshNode* theN1,
+ const bool theHasAngles,
+ list<double>& theAngles,
+ const bool theLinearVariation,
+ const bool theHasRefPoint,
+ const gp_Pnt& theRefPoint,
+ const bool theMakeGroups)
+{
+ MESSAGE("ExtrusionAlongTrack");
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ int aNbE;
+ std::list<double> aPrms;
+ TIDSortedElemSet::iterator itElem;
+
+ gp_XYZ aGC;
+ TopoDS_Edge aTrackEdge;
+ TopoDS_Vertex aV1, aV2;
+
+ SMDS_ElemIteratorPtr aItE;
+ SMDS_NodeIteratorPtr aItN;
+ SMDSAbs_ElementType aTypeE;
+
+ TNodeOfNodeListMap mapNewNodes;
+
+ // 1. Check data
+ aNbE = theElements[0].size() + theElements[1].size();
+ // nothing to do
+ if ( !aNbE )
+ return EXTR_NO_ELEMENTS;
+
+ // 1.1 Track Pattern
+ ASSERT( theTrack );
+
+ SMESHDS_SubMesh* pSubMeshDS = theTrack->GetSubMeshDS();
+ if ( !pSubMeshDS )
+ return ExtrusionAlongTrack( theElements, theTrack->GetFather(), theN1,
+ theHasAngles, theAngles, theLinearVariation,
+ theHasRefPoint, theRefPoint, theMakeGroups );
+
+ aItE = pSubMeshDS->GetElements();
+ while ( aItE->more() ) {
+ const SMDS_MeshElement* pE = aItE->next();
+ aTypeE = pE->GetType();
+ // Pattern must contain links only
+ if ( aTypeE != SMDSAbs_Edge )
+ return EXTR_PATH_NOT_EDGE;
+ }
+
+ list<SMESH_MeshEditor_PathPoint> fullList;
+
+ const TopoDS_Shape& aS = theTrack->GetSubShape();
+ // Sub-shape for the Pattern must be an Edge or Wire
+ if( aS.ShapeType() == TopAbs_EDGE ) {
+ aTrackEdge = TopoDS::Edge( aS );
+ // the Edge must not be degenerated
+ if ( SMESH_Algo::isDegenerated( aTrackEdge ) )
+ return EXTR_BAD_PATH_SHAPE;
+ TopExp::Vertices( aTrackEdge, aV1, aV2 );
+ aItN = theTrack->GetFather()->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN1 = aItN->next();
+ aItN = theTrack->GetFather()->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN2 = aItN->next();
+ // starting node must be aN1 or aN2
+ if ( !( aN1 == theN1 || aN2 == theN1 ) )
+ return EXTR_BAD_STARTING_NODE;
+ aItN = pSubMeshDS->GetNodes();
+ while ( aItN->more() ) {
+ const SMDS_MeshNode* pNode = aItN->next();
+ const SMDS_EdgePosition* pEPos =
+ static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
+ double aT = pEPos->GetUParameter();
+ aPrms.push_back( aT );
+ }
+ //Extrusion_Error err =
+ MakeEdgePathPoints(aPrms, aTrackEdge, (aN1==theN1), fullList);
+ } else if( aS.ShapeType() == TopAbs_WIRE ) {
+ list< SMESH_subMesh* > LSM;
+ TopTools_SequenceOfShape Edges;
+ SMESH_subMeshIteratorPtr itSM = theTrack->getDependsOnIterator(false,true);
+ while(itSM->more()) {
+ SMESH_subMesh* SM = itSM->next();
+ LSM.push_back(SM);
+ const TopoDS_Shape& aS = SM->GetSubShape();
+ Edges.Append(aS);
+ }
+ list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
+ int startNid = theN1->GetID();
+ TColStd_MapOfInteger UsedNums;
+
+ int NbEdges = Edges.Length();
+ int i = 1;
+ for(; i<=NbEdges; i++) {
+ int k = 0;
+ list< SMESH_subMesh* >::iterator itLSM = LSM.begin();
+ for(; itLSM!=LSM.end(); itLSM++) {
+ k++;
+ if(UsedNums.Contains(k)) continue;
+ aTrackEdge = TopoDS::Edge( Edges.Value(k) );
+ SMESH_subMesh* locTrack = *itLSM;
+ SMESHDS_SubMesh* locMeshDS = locTrack->GetSubMeshDS();
+ TopExp::Vertices( aTrackEdge, aV1, aV2 );
+ aItN = locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN1 = aItN->next();
+ aItN = locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS()->GetNodes();
+ const SMDS_MeshNode* aN2 = aItN->next();
+ // starting node must be aN1 or aN2
+ if ( !( aN1->GetID() == startNid || aN2->GetID() == startNid ) ) continue;
+ // 2. Collect parameters on the track edge
+ aPrms.clear();
+ aItN = locMeshDS->GetNodes();
+ while ( aItN->more() ) {
+ const SMDS_MeshNode* pNode = aItN->next();
+ const SMDS_EdgePosition* pEPos =
+ static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
+ double aT = pEPos->GetUParameter();
+ aPrms.push_back( aT );
+ }
+ list<SMESH_MeshEditor_PathPoint> LPP;
+ //Extrusion_Error err =
+ MakeEdgePathPoints(aPrms, aTrackEdge,(aN1->GetID()==startNid), LPP);
+ LLPPs.push_back(LPP);
+ UsedNums.Add(k);
+ // update startN for search following egde
+ if( aN1->GetID() == startNid ) startNid = aN2->GetID();
+ else startNid = aN1->GetID();
+ break;
+ }
+ }
+ list< list<SMESH_MeshEditor_PathPoint> >::iterator itLLPP = LLPPs.begin();
+ list<SMESH_MeshEditor_PathPoint> firstList = *itLLPP;
+ list<SMESH_MeshEditor_PathPoint>::iterator itPP = firstList.begin();
+ for(; itPP!=firstList.end(); itPP++) {
+ fullList.push_back( *itPP );
+ }
+ SMESH_MeshEditor_PathPoint PP1 = fullList.back();
+ fullList.pop_back();
+ itLLPP++;
+ for(; itLLPP!=LLPPs.end(); itLLPP++) {
+ list<SMESH_MeshEditor_PathPoint> currList = *itLLPP;
+ itPP = currList.begin();
+ SMESH_MeshEditor_PathPoint PP2 = currList.front();
+ gp_Dir D1 = PP1.Tangent();
+ gp_Dir D2 = PP2.Tangent();
+ gp_Dir Dnew( gp_Vec( (D1.X()+D2.X())/2, (D1.Y()+D2.Y())/2,
+ (D1.Z()+D2.Z())/2 ) );
+ PP1.SetTangent(Dnew);
+ fullList.push_back(PP1);
+ itPP++;
+ for(; itPP!=firstList.end(); itPP++) {
+ fullList.push_back( *itPP );
+ }
+ PP1 = fullList.back();
+ fullList.pop_back();
+ }
+ // if wire not closed
+ fullList.push_back(PP1);
+ // else ???
+ }
+ else {
+ return EXTR_BAD_PATH_SHAPE;
+ }
+
+ return MakeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
+ theHasRefPoint, theRefPoint, theMakeGroups);
+}
+
+
+//=======================================================================
+//function : ExtrusionAlongTrack
+//purpose :
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet theElements[2],
+ SMESH_Mesh* theTrack,
+ const SMDS_MeshNode* theN1,
+ const bool theHasAngles,
+ list<double>& theAngles,
+ const bool theLinearVariation,
+ const bool theHasRefPoint,
+ const gp_Pnt& theRefPoint,
+ const bool theMakeGroups)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ int aNbE;
+ std::list<double> aPrms;
+ TIDSortedElemSet::iterator itElem;
+
+ gp_XYZ aGC;
+ TopoDS_Edge aTrackEdge;
+ TopoDS_Vertex aV1, aV2;
+
+ SMDS_ElemIteratorPtr aItE;
+ SMDS_NodeIteratorPtr aItN;
+ SMDSAbs_ElementType aTypeE;
+
+ TNodeOfNodeListMap mapNewNodes;
+
+ // 1. Check data
+ aNbE = theElements[0].size() + theElements[1].size();
+ // nothing to do
+ if ( !aNbE )
+ return EXTR_NO_ELEMENTS;
+
+ // 1.1 Track Pattern
+ ASSERT( theTrack );
+
+ SMESHDS_Mesh* pMeshDS = theTrack->GetMeshDS();
+
+ aItE = pMeshDS->elementsIterator();
+ while ( aItE->more() ) {
+ const SMDS_MeshElement* pE = aItE->next();
+ aTypeE = pE->GetType();
+ // Pattern must contain links only
+ if ( aTypeE != SMDSAbs_Edge )
+ return EXTR_PATH_NOT_EDGE;
+ }
+
+ list<SMESH_MeshEditor_PathPoint> fullList;
+
+ const TopoDS_Shape& aS = theTrack->GetShapeToMesh();
+
+ if ( !theTrack->HasShapeToMesh() ) {
+ //Mesh without shape
+ const SMDS_MeshNode* currentNode = NULL;
+ const SMDS_MeshNode* prevNode = theN1;
+ std::vector<const SMDS_MeshNode*> aNodesList;
+ aNodesList.push_back(theN1);
+ int nbEdges = 0, conn=0;
+ const SMDS_MeshElement* prevElem = NULL;
+ const SMDS_MeshElement* currentElem = NULL;
+ int totalNbEdges = theTrack->NbEdges();
+ SMDS_ElemIteratorPtr nIt;
+
+ //check start node
+ if( !theTrack->GetMeshDS()->Contains(theN1) ) {
+ return EXTR_BAD_STARTING_NODE;
+ }
+
+ conn = nbEdgeConnectivity(theN1);
+ if( conn != 1 )
+ return EXTR_PATH_NOT_EDGE;
+
+ aItE = theN1->GetInverseElementIterator();
+ prevElem = aItE->next();
+ currentElem = prevElem;
+ //Get all nodes
+ if(totalNbEdges == 1 ) {
+ nIt = currentElem->nodesIterator();
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ if(currentNode == prevNode)
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ aNodesList.push_back(currentNode);
+ } else {
+ nIt = currentElem->nodesIterator();
+ while( nIt->more() ) {
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ if(currentNode == prevNode)
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ aNodesList.push_back(currentNode);
+
+ //case of the closed mesh
+ if(currentNode == theN1) {
+ nbEdges++;
+ break;
+ }
+
+ conn = nbEdgeConnectivity(currentNode);
+ if(conn > 2) {
+ return EXTR_PATH_NOT_EDGE;
+ }else if( conn == 1 && nbEdges > 0 ) {
+ //End of the path
+ nbEdges++;
+ break;
+ }else {
+ prevNode = currentNode;
+ aItE = currentNode->GetInverseElementIterator();
+ currentElem = aItE->next();
+ if( currentElem == prevElem)
+ currentElem = aItE->next();
+ nIt = currentElem->nodesIterator();
+ prevElem = currentElem;
+ nbEdges++;
+ }
+ }
+ }
+
+ if(nbEdges != totalNbEdges)
+ return EXTR_PATH_NOT_EDGE;
+
+ TopTools_SequenceOfShape Edges;
+ list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
+ int startNid = theN1->GetID();
+ for ( size_t i = 1; i < aNodesList.size(); i++ )
+ {
+ gp_Pnt p1 = SMESH_TNodeXYZ( aNodesList[i-1] );
+ gp_Pnt p2 = SMESH_TNodeXYZ( aNodesList[i] );
+ TopoDS_Edge e = BRepBuilderAPI_MakeEdge( p1, p2 );
+ list<SMESH_MeshEditor_PathPoint> LPP;
+ aPrms.clear();
+ MakeEdgePathPoints(aPrms, e, (aNodesList[i-1]->GetID()==startNid), LPP);
+ LLPPs.push_back(LPP);
+ if ( aNodesList[i-1]->GetID() == startNid ) startNid = aNodesList[i ]->GetID();
+ else startNid = aNodesList[i-1]->GetID();
+ }
+
+ list< list<SMESH_MeshEditor_PathPoint> >::iterator itLLPP = LLPPs.begin();
+ list<SMESH_MeshEditor_PathPoint> firstList = *itLLPP;
+ list<SMESH_MeshEditor_PathPoint>::iterator itPP = firstList.begin();
+ for(; itPP!=firstList.end(); itPP++) {
+ fullList.push_back( *itPP );
+ }
+
+ SMESH_MeshEditor_PathPoint PP1 = fullList.back();
+ SMESH_MeshEditor_PathPoint PP2;
+ fullList.pop_back();
+ itLLPP++;
+ for(; itLLPP!=LLPPs.end(); itLLPP++) {
+ list<SMESH_MeshEditor_PathPoint> currList = *itLLPP;
+ itPP = currList.begin();
+ PP2 = currList.front();
+ gp_Dir D1 = PP1.Tangent();
+ gp_Dir D2 = PP2.Tangent();
+ gp_Dir Dnew( 0.5 * ( D1.XYZ() + D2.XYZ() ));
+ PP1.SetTangent(Dnew);
+ fullList.push_back(PP1);
+ itPP++;
+ for(; itPP!=currList.end(); itPP++) {
+ fullList.push_back( *itPP );
+ }
+ PP1 = fullList.back();
+ fullList.pop_back();
+ }
+ fullList.push_back(PP1);
+
+ } // Sub-shape for the Pattern must be an Edge or Wire
+ else if ( aS.ShapeType() == TopAbs_EDGE )
+ {
+ aTrackEdge = TopoDS::Edge( aS );
+ // the Edge must not be degenerated
+ if ( SMESH_Algo::isDegenerated( aTrackEdge ) )
+ return EXTR_BAD_PATH_SHAPE;
+ TopExp::Vertices( aTrackEdge, aV1, aV2 );
+ const SMDS_MeshNode* aN1 = SMESH_Algo::VertexNode( aV1, pMeshDS );
+ const SMDS_MeshNode* aN2 = SMESH_Algo::VertexNode( aV2, pMeshDS );
+ // starting node must be aN1 or aN2
+ if ( !( aN1 == theN1 || aN2 == theN1 ) )
+ return EXTR_BAD_STARTING_NODE;
+ aItN = pMeshDS->nodesIterator();
+ while ( aItN->more() ) {
+ const SMDS_MeshNode* pNode = aItN->next();
+ if( pNode==aN1 || pNode==aN2 ) continue;
+ const SMDS_EdgePosition* pEPos =
+ static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
+ double aT = pEPos->GetUParameter();
+ aPrms.push_back( aT );
+ }
+ //Extrusion_Error err =
+ MakeEdgePathPoints(aPrms, aTrackEdge, (aN1==theN1), fullList);
+ }
+ else if( aS.ShapeType() == TopAbs_WIRE ) {
+ list< SMESH_subMesh* > LSM;
+ TopTools_SequenceOfShape Edges;
+ TopExp_Explorer eExp(aS, TopAbs_EDGE);
+ for(; eExp.More(); eExp.Next()) {
+ TopoDS_Edge E = TopoDS::Edge( eExp.Current() );
+ if( SMESH_Algo::isDegenerated(E) ) continue;
+ SMESH_subMesh* SM = theTrack->GetSubMesh(E);
+ if(SM) {
+ LSM.push_back(SM);
+ Edges.Append(E);
+ }
+ }
+ list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
+ TopoDS_Vertex aVprev;
+ TColStd_MapOfInteger UsedNums;
+ int NbEdges = Edges.Length();
+ int i = 1;
+ for(; i<=NbEdges; i++) {
+ int k = 0;
+ list< SMESH_subMesh* >::iterator itLSM = LSM.begin();
+ for(; itLSM!=LSM.end(); itLSM++) {
+ k++;
+ if(UsedNums.Contains(k)) continue;
+ aTrackEdge = TopoDS::Edge( Edges.Value(k) );
+ SMESH_subMesh* locTrack = *itLSM;
+ SMESHDS_SubMesh* locMeshDS = locTrack->GetSubMeshDS();
+ TopExp::Vertices( aTrackEdge, aV1, aV2 );
+ bool aN1isOK = false, aN2isOK = false;
+ if ( aVprev.IsNull() ) {
+ // if previous vertex is not yet defined, it means that we in the beginning of wire
+ // and we have to find initial vertex corresponding to starting node theN1
+ const SMDS_MeshNode* aN1 = SMESH_Algo::VertexNode( aV1, pMeshDS );
+ const SMDS_MeshNode* aN2 = SMESH_Algo::VertexNode( aV2, pMeshDS );
+ // starting node must be aN1 or aN2
+ aN1isOK = ( aN1 && aN1 == theN1 );
+ aN2isOK = ( aN2 && aN2 == theN1 );
+ }
+ else {
+ // we have specified ending vertex of the previous edge on the previous iteration
+ // and we have just to check that it corresponds to any vertex in current segment
+ aN1isOK = aVprev.IsSame( aV1 );
+ aN2isOK = aVprev.IsSame( aV2 );
+ }
+ if ( !aN1isOK && !aN2isOK ) continue;
+ // 2. Collect parameters on the track edge
+ aPrms.clear();
+ aItN = locMeshDS->GetNodes();
+ while ( aItN->more() ) {
+ const SMDS_MeshNode* pNode = aItN->next();
+ const SMDS_EdgePosition* pEPos =
+ static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
+ double aT = pEPos->GetUParameter();
+ aPrms.push_back( aT );
+ }
+ list<SMESH_MeshEditor_PathPoint> LPP;
+ //Extrusion_Error err =
+ MakeEdgePathPoints(aPrms, aTrackEdge, aN1isOK, LPP);
+ LLPPs.push_back(LPP);
+ UsedNums.Add(k);
+ // update startN for search following egde
+ if ( aN1isOK ) aVprev = aV2;
+ else aVprev = aV1;
+ break;
+ }
+ }
+ list< list<SMESH_MeshEditor_PathPoint> >::iterator itLLPP = LLPPs.begin();
+ list<SMESH_MeshEditor_PathPoint>& firstList = *itLLPP;
+ fullList.splice( fullList.end(), firstList );
+
+ SMESH_MeshEditor_PathPoint PP1 = fullList.back();
+ fullList.pop_back();
+ itLLPP++;
+ for(; itLLPP!=LLPPs.end(); itLLPP++) {
+ list<SMESH_MeshEditor_PathPoint>& currList = *itLLPP;
+ SMESH_MeshEditor_PathPoint PP2 = currList.front();
+ gp_Dir D1 = PP1.Tangent();
+ gp_Dir D2 = PP2.Tangent();
+ gp_Dir Dnew( D1.XYZ() + D2.XYZ() );
+ PP1.SetTangent(Dnew);
+ fullList.push_back(PP1);
+ fullList.splice( fullList.end(), currList, ++currList.begin(), currList.end() );
+ PP1 = fullList.back();
+ fullList.pop_back();
+ }
+ // if wire not closed
+ fullList.push_back(PP1);
+ // else ???
+ }
+ else {
+ return EXTR_BAD_PATH_SHAPE;
+ }
+
+ return MakeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
+ theHasRefPoint, theRefPoint, theMakeGroups);
+}
+
+
+//=======================================================================
+//function : MakeEdgePathPoints
+//purpose : auxilary for ExtrusionAlongTrack
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::MakeEdgePathPoints(std::list<double>& aPrms,
+ const TopoDS_Edge& aTrackEdge,
+ bool FirstIsStart,
+ list<SMESH_MeshEditor_PathPoint>& LPP)
+{
+ Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2;
+ aTolVec=1.e-7;
+ aTolVec2=aTolVec*aTolVec;
+ double aT1, aT2;
+ TopoDS_Vertex aV1, aV2;
+ TopExp::Vertices( aTrackEdge, aV1, aV2 );
+ aT1=BRep_Tool::Parameter( aV1, aTrackEdge );
+ aT2=BRep_Tool::Parameter( aV2, aTrackEdge );
+ // 2. Collect parameters on the track edge
+ aPrms.push_front( aT1 );
+ aPrms.push_back( aT2 );
+ // sort parameters
+ aPrms.sort();
+ if( FirstIsStart ) {
+ if ( aT1 > aT2 ) {
+ aPrms.reverse();
+ }
+ }
+ else {
+ if ( aT2 > aT1 ) {
+ aPrms.reverse();
+ }
+ }
+ // 3. Path Points
+ SMESH_MeshEditor_PathPoint aPP;
+ Handle(Geom_Curve) aC3D = BRep_Tool::Curve( aTrackEdge, aTx1, aTx2 );
+ std::list<double>::iterator aItD = aPrms.begin();
+ for(; aItD != aPrms.end(); ++aItD) {
+ double aT = *aItD;
+ gp_Pnt aP3D;
+ gp_Vec aVec;
+ aC3D->D1( aT, aP3D, aVec );
+ aL2 = aVec.SquareMagnitude();
+ if ( aL2 < aTolVec2 )
+ return EXTR_CANT_GET_TANGENT;
+ gp_Dir aTgt( FirstIsStart ? aVec : -aVec );
+ aPP.SetPnt( aP3D );
+ aPP.SetTangent( aTgt );
+ aPP.SetParameter( aT );
+ LPP.push_back(aPP);
+ }
+ return EXTR_OK;
+}
+
+
+//=======================================================================
+//function : MakeExtrElements
+//purpose : auxilary for ExtrusionAlongTrack
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet theElemSets[2],
+ list<SMESH_MeshEditor_PathPoint>& fullList,
+ const bool theHasAngles,
+ list<double>& theAngles,
+ const bool theLinearVariation,
+ const bool theHasRefPoint,
+ const gp_Pnt& theRefPoint,
+ const bool theMakeGroups)
+{
+ const int aNbTP = fullList.size();
+
+ // Angles
+ if( theHasAngles && !theAngles.empty() && theLinearVariation )
+ LinearAngleVariation(aNbTP-1, theAngles);
+
+ // fill vector of path points with angles
+ vector<SMESH_MeshEditor_PathPoint> aPPs;
+ list<SMESH_MeshEditor_PathPoint>::iterator itPP = fullList.begin();
+ list<double>::iterator itAngles = theAngles.begin();
+ aPPs.push_back( *itPP++ );
+ for( ; itPP != fullList.end(); itPP++) {
+ aPPs.push_back( *itPP );
+ if ( theHasAngles && itAngles != theAngles.end() )
+ aPPs.back().SetAngle( *itAngles++ );
+ }
+
+ TNodeOfNodeListMap mapNewNodes;
+ TElemOfVecOfNnlmiMap mapElemNewNodes;
+ TTElemOfElemListMap newElemsMap;
+ TIDSortedElemSet::iterator itElem;
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+
+ // 3. Center of rotation aV0
+ gp_Pnt aV0 = theRefPoint;
+ if ( !theHasRefPoint )
+ {
+ gp_XYZ aGC( 0.,0.,0. );
+ TIDSortedElemSet newNodes;
+
+ for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+ {
+ TIDSortedElemSet& theElements = theElemSets[ is2ndSet ];
+ itElem = theElements.begin();
+ for ( ; itElem != theElements.end(); itElem++ )
+ {
+ const SMDS_MeshElement* elem = *itElem;
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() ) {
+ const SMDS_MeshElement* node = itN->next();
+ if ( newNodes.insert( node ).second )
+ aGC += SMESH_TNodeXYZ( node );
+ }
+ }
+ }
+ aGC /= newNodes.size();
+ aV0.SetXYZ( aGC );
+ } // if (!theHasRefPoint) {
+
+ // 4. Processing the elements
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+ list<const SMDS_MeshNode*> emptyList;
+
+ setElemsFirst( theElemSets );
+ for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+ {
+ TIDSortedElemSet& theElements = theElemSets[ is2ndSet ];
+ for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ )
+ {
+ const SMDS_MeshElement* elem = *itElem;
+
+ vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
+ newNodesItVec.reserve( elem->NbNodes() );
+
+ // loop on elem nodes
+ int nodeIndex = -1;
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() )
+ {
+ ++nodeIndex;
+ // check if a node has been already processed
+ const SMDS_MeshNode* node = cast2Node( itN->next() );
+ TNodeOfNodeListMap::iterator nIt = mapNewNodes.insert( make_pair( node, emptyList )).first;
+ list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
+ if ( listNewNodes.empty() )
+ {
+ // make new nodes
+ Standard_Real aAngle1x, aAngleT1T0, aTolAng;
+ gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x;
+ gp_Ax1 anAx1, anAxT1T0;
+ gp_Dir aDT1x, aDT0x, aDT1T0;
+
+ aTolAng=1.e-4;
+
+ aV0x = aV0;
+ aPN0 = SMESH_TNodeXYZ( node );
+
+ const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
+ aP0x = aPP0.Pnt();
+ aDT0x= aPP0.Tangent();
+
+ for ( int j = 1; j < aNbTP; ++j ) {
+ const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j];
+ aP1x = aPP1.Pnt();
+ aDT1x = aPP1.Tangent();
+ aAngle1x = aPP1.Angle();
+
+ gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
+ // Translation
+ gp_Vec aV01x( aP0x, aP1x );
+ aTrsf.SetTranslation( aV01x );
+
+ // traslated point
+ aV1x = aV0x.Transformed( aTrsf );
+ aPN1 = aPN0.Transformed( aTrsf );
+
+ // rotation 1 [ T1,T0 ]
+ aAngleT1T0=-aDT1x.Angle( aDT0x );
+ if (fabs(aAngleT1T0) > aTolAng)
+ {
+ aDT1T0=aDT1x^aDT0x;
+ anAxT1T0.SetLocation( aV1x );
+ anAxT1T0.SetDirection( aDT1T0 );
+ aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 );
+
+ aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
+ }
+
+ // rotation 2
+ if ( theHasAngles ) {
+ anAx1.SetLocation( aV1x );
+ anAx1.SetDirection( aDT1x );
+ aTrsfRot.SetRotation( anAx1, aAngle1x );
+
+ aPN1 = aPN1.Transformed( aTrsfRot );
+ }
+
+ // make new node
+ if ( elem->IsQuadratic() && !elem->IsMediumNode(node) )
+ {
+ // create additional node
+ gp_XYZ midP = 0.5 * ( aPN1.XYZ() + aPN0.XYZ() );
+ const SMDS_MeshNode* newNode = aMesh->AddNode( midP.X(), midP.Y(), midP.Z() );
+ myLastCreatedNodes.Append(newNode);
+ srcNodes.Append( node );
+ listNewNodes.push_back( newNode );
+ }
+ const SMDS_MeshNode* newNode = aMesh->AddNode( aPN1.X(), aPN1.Y(), aPN1.Z() );
+ myLastCreatedNodes.Append(newNode);
+ srcNodes.Append( node );
+ listNewNodes.push_back( newNode );
+
+ aPN0 = aPN1;
+ aP0x = aP1x;
+ aV0x = aV1x;
+ aDT0x = aDT1x;
+ }
+ }
+ else if( elem->IsQuadratic() && !elem->IsMediumNode(node) )
+ {
+ // if current elem is quadratic and current node is not medium
+ // we have to check - may be it is needed to insert additional nodes
+ list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
+ if ((int) listNewNodes.size() == aNbTP-1 )
+ {
+ vector<const SMDS_MeshNode*> aNodes(2*(aNbTP-1));
+ gp_XYZ P(node->X(), node->Y(), node->Z());
+ list< const SMDS_MeshNode* >::iterator it = listNewNodes.begin();
+ int i;
+ for(i=0; i<aNbTP-1; i++) {
+ const SMDS_MeshNode* N = *it;
+ double x = ( N->X() + P.X() )/2.;
+ double y = ( N->Y() + P.Y() )/2.;
+ double z = ( N->Z() + P.Z() )/2.;
+ const SMDS_MeshNode* newN = aMesh->AddNode(x,y,z);
+ srcNodes.Append( node );
+ myLastCreatedNodes.Append(newN);
+ aNodes[2*i] = newN;
+ aNodes[2*i+1] = N;
+ P = gp_XYZ(N->X(),N->Y(),N->Z());
+ }
+ listNewNodes.clear();
+ for(i=0; i<2*(aNbTP-1); i++) {
+ listNewNodes.push_back(aNodes[i]);
+ }
+ }
+ }
+
+ newNodesItVec.push_back( nIt );
+ }
+
+ // make new elements
+ sweepElement( elem, newNodesItVec, newElemsMap[elem], aNbTP-1, srcElems );
+ }
+ }
+
+ makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElemSets[0], aNbTP-1, srcElems );
+
+ if ( theMakeGroups )
+ generateGroups( srcNodes, srcElems, "extruded");
+
+ return EXTR_OK;
+}
+
+
+//=======================================================================
+//function : LinearAngleVariation
+//purpose : auxilary for ExtrusionAlongTrack
+//=======================================================================
+void SMESH_MeshEditor::LinearAngleVariation(const int nbSteps,
+ list<double>& Angles)
+{
+ int nbAngles = Angles.size();
+ if( nbSteps > nbAngles ) {
+ vector<double> theAngles(nbAngles);
+ list<double>::iterator it = Angles.begin();
+ int i = -1;
+ for(; it!=Angles.end(); it++) {
+ i++;
+ theAngles[i] = (*it);
+ }
+ list<double> res;
+ double rAn2St = double( nbAngles ) / double( nbSteps );
+ double angPrev = 0, angle;
+ for ( int iSt = 0; iSt < nbSteps; ++iSt ) {
+ double angCur = rAn2St * ( iSt+1 );
+ double angCurFloor = floor( angCur );
+ double angPrevFloor = floor( angPrev );
+ if ( angPrevFloor == angCurFloor )
+ angle = rAn2St * theAngles[ int( angCurFloor ) ];
+ else {
+ int iP = int( angPrevFloor );
+ double angPrevCeil = ceil(angPrev);
+ angle = ( angPrevCeil - angPrev ) * theAngles[ iP ];
+
+ int iC = int( angCurFloor );
+ if ( iC < nbAngles )
+ angle += ( angCur - angCurFloor ) * theAngles[ iC ];
+
+ iP = int( angPrevCeil );
+ while ( iC-- > iP )
+ angle += theAngles[ iC ];
+ }
+ res.push_back(angle);
+ angPrev = angCur;
+ }
+ Angles.clear();
+ it = res.begin();
+ for(; it!=res.end(); it++)
+ Angles.push_back( *it );
+ }
+}
+
+
+//================================================================================
+/*!
+ * \brief Move or copy theElements applying theTrsf to their nodes
+ * \param theElems - elements to transform, if theElems is empty then apply to all mesh nodes
+ * \param theTrsf - transformation to apply
+ * \param theCopy - if true, create translated copies of theElems
+ * \param theMakeGroups - if true and theCopy, create translated groups
+ * \param theTargetMesh - mesh to copy translated elements into
+ * \return SMESH_MeshEditor::PGroupIDs - list of ids of created groups
+ */
+//================================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
+ const gp_Trsf& theTrsf,
+ const bool theCopy,
+ const bool theMakeGroups,
+ SMESH_Mesh* theTargetMesh)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ bool needReverse = false;
+ string groupPostfix;
+ switch ( theTrsf.Form() ) {
+ case gp_PntMirror:
+ MESSAGE("gp_PntMirror");
+ needReverse = true;
+ groupPostfix = "mirrored";
+ break;
+ case gp_Ax1Mirror:
+ MESSAGE("gp_Ax1Mirror");
+ groupPostfix = "mirrored";
+ break;
+ case gp_Ax2Mirror:
+ MESSAGE("gp_Ax2Mirror");
+ needReverse = true;
+ groupPostfix = "mirrored";
+ break;
+ case gp_Rotation:
+ MESSAGE("gp_Rotation");
+ groupPostfix = "rotated";
+ break;
+ case gp_Translation:
+ MESSAGE("gp_Translation");
+ groupPostfix = "translated";
+ break;
+ case gp_Scale:
+ MESSAGE("gp_Scale");
+ groupPostfix = "scaled";
+ break;
+ case gp_CompoundTrsf: // different scale by axis
+ MESSAGE("gp_CompoundTrsf");
+ groupPostfix = "scaled";
+ break;
+ default:
+ MESSAGE("default");
+ needReverse = false;
+ groupPostfix = "transformed";
+ }
+
+ SMESHDS_Mesh* aTgtMesh = theTargetMesh ? theTargetMesh->GetMeshDS() : 0;
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ SMESH_MeshEditor targetMeshEditor( theTargetMesh );
+ SMESH_MeshEditor* editor = theTargetMesh ? & targetMeshEditor : theCopy ? this : 0;
+ SMESH_MeshEditor::ElemFeatures elemType;
+
+ // map old node to new one
+ TNodeNodeMap nodeMap;
+
+ // elements sharing moved nodes; those of them which have all
+ // nodes mirrored but are not in theElems are to be reversed
+ TIDSortedElemSet inverseElemSet;
+
+ // source elements for each generated one
+ SMESH_SequenceOfElemPtr srcElems, srcNodes;
+
+ // issue 021015: EDF 1578 SMESH: Free nodes are removed when translating a mesh
+ TIDSortedElemSet orphanNode;
+
+ if ( theElems.empty() ) // transform the whole mesh
+ {
+ // add all elements
+ SMDS_ElemIteratorPtr eIt = aMesh->elementsIterator();
+ while ( eIt->more() ) theElems.insert( eIt->next() );
+ // add orphan nodes
+ SMDS_NodeIteratorPtr nIt = aMesh->nodesIterator();
+ while ( nIt->more() )
+ {
+ const SMDS_MeshNode* node = nIt->next();
+ if ( node->NbInverseElements() == 0)
+ orphanNode.insert( node );
+ }
+ }
+
+ // loop on elements to transform nodes : first orphan nodes then elems
+ TIDSortedElemSet::iterator itElem;
+ TIDSortedElemSet *elements[] = { &orphanNode, &theElems };
+ for (int i=0; i<2; i++)
+ for ( itElem = elements[i]->begin(); itElem != elements[i]->end(); itElem++ )
+ {
+ const SMDS_MeshElement* elem = *itElem;
+ if ( !elem )
+ continue;
+
+ // loop on elem nodes
+ double coord[3];
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() )
+ {
+ const SMDS_MeshNode* node = cast2Node( itN->next() );
+ // check if a node has been already transformed
+ pair<TNodeNodeMap::iterator,bool> n2n_isnew =
+ nodeMap.insert( make_pair ( node, node ));
+ if ( !n2n_isnew.second )
+ continue;
+
+ node->GetXYZ( coord );
+ theTrsf.Transforms( coord[0], coord[1], coord[2] );
+ if ( theTargetMesh ) {
+ const SMDS_MeshNode * newNode = aTgtMesh->AddNode( coord[0], coord[1], coord[2] );
+ n2n_isnew.first->second = newNode;
+ myLastCreatedNodes.Append(newNode);
+ srcNodes.Append( node );
+ }
+ else if ( theCopy ) {
+ const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+ n2n_isnew.first->second = newNode;
+ myLastCreatedNodes.Append(newNode);
+ srcNodes.Append( node );
+ }
+ else {
+ aMesh->MoveNode( node, coord[0], coord[1], coord[2] );
+ // node position on shape becomes invalid
+ const_cast< SMDS_MeshNode* > ( node )->SetPosition
+ ( SMDS_SpacePosition::originSpacePosition() );
+ }
+
+ // keep inverse elements
+ if ( !theCopy && !theTargetMesh && needReverse ) {
+ SMDS_ElemIteratorPtr invElemIt = node->GetInverseElementIterator();
+ while ( invElemIt->more() ) {
+ const SMDS_MeshElement* iel = invElemIt->next();
+ inverseElemSet.insert( iel );
+ }
+ }
+ }
+ } // loop on elems in { &orphanNode, &theElems };
+
+ // either create new elements or reverse mirrored ones
+ if ( !theCopy && !needReverse && !theTargetMesh )
+ return PGroupIDs();
+
+ theElems.insert( inverseElemSet.begin(),inverseElemSet.end() );
+
+ // Replicate or reverse elements
+
+ std::vector<int> iForw;
+ vector<const SMDS_MeshNode*> nodes;
+ for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
+ {
+ const SMDS_MeshElement* elem = *itElem;
+ if ( !elem ) continue;
+
+ SMDSAbs_GeometryType geomType = elem->GetGeomType();
+ size_t nbNodes = elem->NbNodes();
+ if ( geomType == SMDSGeom_NONE ) continue; // node
+
+ nodes.resize( nbNodes );
+
+ if ( geomType == SMDSGeom_POLYHEDRA ) // ------------------ polyhedral volume
+ {
+ const SMDS_VtkVolume* aPolyedre = dynamic_cast<const SMDS_VtkVolume*>( elem );
+ if (!aPolyedre)
+ continue;
+ nodes.clear();
+ bool allTransformed = true;
+ int nbFaces = aPolyedre->NbFaces();
+ for (int iface = 1; iface <= nbFaces && allTransformed; iface++)
+ {
+ int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
+ for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++)
+ {
+ const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode);
+ TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
+ if ( nodeMapIt == nodeMap.end() )
+ allTransformed = false; // not all nodes transformed
+ else
+ nodes.push_back((*nodeMapIt).second);
+ }
+ if ( needReverse && allTransformed )
+ std::reverse( nodes.end() - nbFaceNodes, nodes.end() );
+ }
+ if ( !allTransformed )
+ continue; // not all nodes transformed
+ }
+ else // ----------------------- the rest element types
+ {
+ while ( iForw.size() < nbNodes ) iForw.push_back( iForw.size() );
+ const vector<int>& iRev = SMDS_MeshCell::reverseSmdsOrder( elem->GetEntityType(), nbNodes );
+ const vector<int>& i = needReverse ? iRev : iForw;
+
+ // find transformed nodes
+ size_t iNode = 0;
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() ) {
+ const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
+ TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
+ if ( nodeMapIt == nodeMap.end() )
+ break; // not all nodes transformed
+ nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
+ }
+ if ( iNode != nbNodes )
+ continue; // not all nodes transformed
+ }
+
+ if ( editor ) {
+ // copy in this or a new mesh
+ if ( editor->AddElement( nodes, elemType.Init( elem, /*basicOnly=*/false )))
+ srcElems.Append( elem );
+ }
+ else {
+ // reverse element as it was reversed by transformation
+ if ( nbNodes > 2 )
+ aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes );
+ }
+
+ } // loop on elements
+
+ if ( editor && editor != this )
+ myLastCreatedElems = editor->myLastCreatedElems;
+
+ PGroupIDs newGroupIDs;
+
+ if ( ( theMakeGroups && theCopy ) ||
+ ( theMakeGroups && theTargetMesh ) )
+ newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh, false );
+
+ return newGroupIDs;
+}
+
+//=======================================================================
+/*!
+ * \brief Create groups of elements made during transformation
+ * \param nodeGens - nodes making corresponding myLastCreatedNodes
+ * \param elemGens - elements making corresponding myLastCreatedElems
+ * \param postfix - to append to names of new groups
+ * \param targetMesh - mesh to create groups in
+ * \param topPresent - is there "top" elements that are created by sweeping
+ */
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
+ const SMESH_SequenceOfElemPtr& elemGens,
+ const std::string& postfix,
+ SMESH_Mesh* targetMesh,
+ const bool topPresent)
+{
+ PGroupIDs newGroupIDs( new list<int> );
+ SMESH_Mesh* mesh = targetMesh ? targetMesh : GetMesh();
+
+ // Sort existing groups by types and collect their names
+
+ // containers to store an old group and generated new ones;
+ // 1st new group is for result elems of different type than a source one;
+ // 2nd new group is for same type result elems ("top" group at extrusion)
+ using boost::tuple;
+ using boost::make_tuple;
+ typedef tuple< SMESHDS_GroupBase*, SMESHDS_Group*, SMESHDS_Group* > TOldNewGroup;
+ vector< list< TOldNewGroup > > groupsByType( SMDSAbs_NbElementTypes );
+ vector< TOldNewGroup* > orderedOldNewGroups; // in order of old groups
+ // group names
+ set< string > groupNames;
+
+ SMESH_Mesh::GroupIteratorPtr groupIt = GetMesh()->GetGroups();
+ if ( !groupIt->more() ) return newGroupIDs;
+
+ int newGroupID = mesh->GetGroupIds().back()+1;
+ while ( groupIt->more() )
+ {
+ SMESH_Group * group = groupIt->next();
+ if ( !group ) continue;
+ SMESHDS_GroupBase* groupDS = group->GetGroupDS();
+ if ( !groupDS || groupDS->IsEmpty() ) continue;
+ groupNames.insert ( group->GetName() );
+ groupDS->SetStoreName( group->GetName() );
+ const SMDSAbs_ElementType type = groupDS->GetType();
+ SMESHDS_Group* newGroup = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(), type );
+ SMESHDS_Group* newTopGroup = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(), type );
+ groupsByType[ type ].push_back( make_tuple( groupDS, newGroup, newTopGroup ));
+ orderedOldNewGroups.push_back( & groupsByType[ type ].back() );
+ }
+
+ // Loop on nodes and elements to add them in new groups
+
+ vector< const SMDS_MeshElement* > resultElems;
+ for ( int isNodes = 0; isNodes < 2; ++isNodes )
+ {
+ const SMESH_SequenceOfElemPtr& gens = isNodes ? nodeGens : elemGens;
+ const SMESH_SequenceOfElemPtr& elems = isNodes ? myLastCreatedNodes : myLastCreatedElems;
+ if ( gens.Length() != elems.Length() )
+ throw SALOME_Exception("SMESH_MeshEditor::generateGroups(): invalid args");
+
+ // loop on created elements
+ for (int iElem = 1; iElem <= elems.Length(); ++iElem )
+ {
+ const SMDS_MeshElement* sourceElem = gens( iElem );
+ if ( !sourceElem ) {
+ MESSAGE("generateGroups(): NULL source element");
+ continue;
+ }
+ list< TOldNewGroup > & groupsOldNew = groupsByType[ sourceElem->GetType() ];
+ if ( groupsOldNew.empty() ) { // no groups of this type at all
+ while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
+ ++iElem; // skip all elements made by sourceElem
+ continue;
+ }
+ // collect all elements made by the iElem-th sourceElem
+ resultElems.clear();
+ if ( const SMDS_MeshElement* resElem = elems( iElem ))
+ if ( resElem != sourceElem )
+ resultElems.push_back( resElem );
+ while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
+ if ( const SMDS_MeshElement* resElem = elems( ++iElem ))
+ if ( resElem != sourceElem )
+ resultElems.push_back( resElem );
+
+ const SMDS_MeshElement* topElem = 0;
+ if ( isNodes ) // there must be a top element
+ {
+ topElem = resultElems.back();
+ resultElems.pop_back();
+ }
+ else
+ {
+ vector< const SMDS_MeshElement* >::reverse_iterator resElemIt = resultElems.rbegin();
+ for ( ; resElemIt != resultElems.rend() ; ++resElemIt )
+ if ( (*resElemIt)->GetType() == sourceElem->GetType() )
+ {
+ topElem = *resElemIt;
+ *resElemIt = 0; // erase *resElemIt
+ break;
+ }
+ }
+ // add resultElems to groups originted from ones the sourceElem belongs to
+ list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end();
+ for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew )
+ {
+ SMESHDS_GroupBase* oldGroup = gOldNew->get<0>();
+ if ( oldGroup->Contains( sourceElem )) // sourceElem is in oldGroup
+ {
+ // fill in a new group
+ SMDS_MeshGroup & newGroup = gOldNew->get<1>()->SMDSGroup();
+ vector< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
+ for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt )
+ if ( *resElemIt )
+ newGroup.Add( *resElemIt );
+
+ // fill a "top" group
+ if ( topElem )
+ {
+ SMDS_MeshGroup & newTopGroup = gOldNew->get<2>()->SMDSGroup();
+ newTopGroup.Add( topElem );
+ }
+ }
+ }
+ } // loop on created elements
+ }// loop on nodes and elements
+
+ // Create new SMESH_Groups from SMESHDS_Groups and remove empty SMESHDS_Groups
+
+ list<int> topGrouIds;
+ for ( size_t i = 0; i < orderedOldNewGroups.size(); ++i )
+ {
+ SMESHDS_GroupBase* oldGroupDS = orderedOldNewGroups[i]->get<0>();
+ SMESHDS_Group* newGroups[2] = { orderedOldNewGroups[i]->get<1>(),
+ orderedOldNewGroups[i]->get<2>() };
+ for ( int is2nd = 0; is2nd < 2; ++is2nd )
+ {
+ SMESHDS_Group* newGroupDS = newGroups[ is2nd ];
+ if ( newGroupDS->IsEmpty() )
+ {
+ mesh->GetMeshDS()->RemoveGroup( newGroupDS );
+ }
+ else
+ {
+ // set group type
+ newGroupDS->SetType( newGroupDS->GetElements()->next()->GetType() );
+
+ // make a name
+ const bool isTop = ( topPresent &&
+ newGroupDS->GetType() == oldGroupDS->GetType() &&
+ is2nd );
+
+ string name = oldGroupDS->GetStoreName();
+ { // remove trailing whitespaces (issue 22599)
+ size_t size = name.size();
+ while ( size > 1 && isspace( name[ size-1 ]))
+ --size;
+ if ( size != name.size() )
+ {
+ name.resize( size );
+ oldGroupDS->SetStoreName( name.c_str() );
+ }
+ }
+ if ( !targetMesh ) {
+ string suffix = ( isTop ? "top": postfix.c_str() );
+ name += "_";
+ name += suffix;
+ int nb = 1;
+ while ( !groupNames.insert( name ).second ) // name exists
+ name = SMESH_Comment( oldGroupDS->GetStoreName() ) << "_" << suffix << "_" << nb++;
+ }
+ else if ( isTop ) {
+ name += "_top";
+ }
+ newGroupDS->SetStoreName( name.c_str() );
+
+ // make a SMESH_Groups
+ mesh->AddGroup( newGroupDS );
+ if ( isTop )
+ topGrouIds.push_back( newGroupDS->GetID() );
+ else
+ newGroupIDs->push_back( newGroupDS->GetID() );
+ }
+ }
+ }
+ newGroupIDs->splice( newGroupIDs->end(), topGrouIds );
+
+ return newGroupIDs;
+}
+
+//================================================================================
+/*!
+ * * \brief Return list of group of nodes close to each other within theTolerance
+ * * Search among theNodes or in the whole mesh if theNodes is empty using
+ * * an Octree algorithm
+ * \param [in,out] theNodes - the nodes to treat
+ * \param [in] theTolerance - the tolerance
+ * \param [out] theGroupsOfNodes - the result groups of coincident nodes
+ * \param [in] theSeparateCornersAndMedium - if \c true, in quadratic mesh puts
+ * corner and medium nodes in separate groups
+ */
+//================================================================================
+
+void SMESH_MeshEditor::FindCoincidentNodes (TIDSortedNodeSet & theNodes,
+ const double theTolerance,
+ TListOfListOfNodes & theGroupsOfNodes,
+ bool theSeparateCornersAndMedium)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ if ( myMesh->NbEdges ( ORDER_QUADRATIC ) +
+ myMesh->NbFaces ( ORDER_QUADRATIC ) +
+ myMesh->NbVolumes( ORDER_QUADRATIC ) == 0 )
+ theSeparateCornersAndMedium = false;
+
+ TIDSortedNodeSet& corners = theNodes;
+ TIDSortedNodeSet medium;
+
+ if ( theNodes.empty() ) // get all nodes in the mesh
+ {
+ TIDSortedNodeSet* nodes[2] = { &corners, &medium };
+ SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator(/*idInceasingOrder=*/true);
+ if ( theSeparateCornersAndMedium )
+ while ( nIt->more() )
+ {
+ const SMDS_MeshNode* n = nIt->next();
+ TIDSortedNodeSet* & nodeSet = nodes[ SMESH_MesherHelper::IsMedium( n )];
+ nodeSet->insert( nodeSet->end(), n );
+ }
+ else
+ while ( nIt->more() )
+ theNodes.insert( theNodes.end(),nIt->next() );
+ }
+ else if ( theSeparateCornersAndMedium ) // separate corners from medium nodes
+ {
+ TIDSortedNodeSet::iterator nIt = corners.begin();
+ while ( nIt != corners.end() )
+ if ( SMESH_MesherHelper::IsMedium( *nIt ))
+ {
+ medium.insert( medium.end(), *nIt );
+ corners.erase( nIt++ );
+ }
+ else
+ {
+ ++nIt;
+ }
+ }
+
+ if ( !corners.empty() )
+ SMESH_OctreeNode::FindCoincidentNodes ( corners, &theGroupsOfNodes, theTolerance );
+ if ( !medium.empty() )
+ SMESH_OctreeNode::FindCoincidentNodes ( medium, &theGroupsOfNodes, theTolerance );
+}
+
+//=======================================================================
+//function : SimplifyFace
+//purpose : split a chain of nodes into several closed chains
+//=======================================================================
+
+int SMESH_MeshEditor::SimplifyFace (const vector<const SMDS_MeshNode *>& faceNodes,
+ vector<const SMDS_MeshNode *>& poly_nodes,
+ vector<int>& quantities) const
+{
+ int nbNodes = faceNodes.size();
+
+ if (nbNodes < 3)
+ return 0;
+
+ set<const SMDS_MeshNode*> nodeSet;
+
+ // get simple seq of nodes
+ vector<const SMDS_MeshNode*> simpleNodes( nbNodes );
+ int iSimple = 0;
+
+ simpleNodes[iSimple++] = faceNodes[0];
+ for (int iCur = 1; iCur < nbNodes; iCur++) {
+ if (faceNodes[iCur] != simpleNodes[iSimple - 1]) {
+ simpleNodes[iSimple++] = faceNodes[iCur];
+ nodeSet.insert( faceNodes[iCur] );
+ }
+ }
+ int nbUnique = nodeSet.size();
+ int nbSimple = iSimple;
+ if (simpleNodes[nbSimple - 1] == simpleNodes[0]) {
+ nbSimple--;
+ iSimple--;
+ }
+
+ if (nbUnique < 3)
+ return 0;
+
+ // separate loops
+ int nbNew = 0;
+ bool foundLoop = (nbSimple > nbUnique);
+ while (foundLoop) {
+ foundLoop = false;
+ set<const SMDS_MeshNode*> loopSet;
+ for (iSimple = 0; iSimple < nbSimple && !foundLoop; iSimple++) {
+ const SMDS_MeshNode* n = simpleNodes[iSimple];
+ if (!loopSet.insert( n ).second) {
+ foundLoop = true;
+
+ // separate loop
+ int iC = 0, curLast = iSimple;
+ for (; iC < curLast; iC++) {
+ if (simpleNodes[iC] == n) break;
+ }
+ int loopLen = curLast - iC;
+ if (loopLen > 2) {
+ // create sub-element
+ nbNew++;
+ quantities.push_back(loopLen);
+ for (; iC < curLast; iC++) {
+ poly_nodes.push_back(simpleNodes[iC]);
+ }
+ }
+ // shift the rest nodes (place from the first loop position)
+ for (iC = curLast + 1; iC < nbSimple; iC++) {
+ simpleNodes[iC - loopLen] = simpleNodes[iC];
+ }
+ nbSimple -= loopLen;
+ iSimple -= loopLen;
+ }
+ } // for (iSimple = 0; iSimple < nbSimple; iSimple++)
+ } // while (foundLoop)
+
+ if (iSimple > 2) {
+ nbNew++;
+ quantities.push_back(iSimple);
+ for (int i = 0; i < iSimple; i++)
+ poly_nodes.push_back(simpleNodes[i]);
+ }
+
+ return nbNew;
+}
+
+//=======================================================================
+//function : MergeNodes
+//purpose : In each group, the cdr of nodes are substituted by the first one
+// in all elements.
+//=======================================================================
+
+void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
+{
+ MESSAGE("MergeNodes");
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ TNodeNodeMap nodeNodeMap; // node to replace - new node
+ set<const SMDS_MeshElement*> elems; // all elements with changed nodes
+ list< int > rmElemIds, rmNodeIds;
+
+ // Fill nodeNodeMap and elems
+
+ TListOfListOfNodes::iterator grIt = theGroupsOfNodes.begin();
+ for ( ; grIt != theGroupsOfNodes.end(); grIt++ )
+ {
+ list<const SMDS_MeshNode*>& nodes = *grIt;
+ list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
+ const SMDS_MeshNode* nToKeep = *nIt;
+ for ( ++nIt; nIt != nodes.end(); nIt++ )
+ {
+ const SMDS_MeshNode* nToRemove = *nIt;
+ nodeNodeMap.insert( make_pair( nToRemove, nToKeep ));
+ if ( nToRemove != nToKeep )
+ {
+ rmNodeIds.push_back( nToRemove->GetID() );
+ AddToSameGroups( nToKeep, nToRemove, aMesh );
+ // set _alwaysComputed to a sub-mesh of VERTEX to enable mesh computing
+ // after MergeNodes() w/o creating node in place of merged ones.
+ const SMDS_PositionPtr& pos = nToRemove->GetPosition();
+ if ( pos && pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
+ if ( SMESH_subMesh* sm = myMesh->GetSubMeshContaining( nToRemove->getshapeId() ))
+ sm->SetIsAlwaysComputed( true );
+ }
+ SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
+ while ( invElemIt->more() ) {
+ const SMDS_MeshElement* elem = invElemIt->next();
+ elems.insert(elem);
+ }
+ }
+ }
+ // Change element nodes or remove an element
+
+ set<const SMDS_MeshNode*> nodeSet;
+ vector< const SMDS_MeshNode*> curNodes, uniqueNodes;
+ vector<int> iRepl;
+ ElemFeatures elemType;
+
+ set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
+ for ( ; eIt != elems.end(); eIt++ )
+ {
+ const SMDS_MeshElement* elem = *eIt;
+ const int nbNodes = elem->NbNodes();
+ const int aShapeId = FindShape( elem );
+
+ nodeSet.clear();
+ curNodes.resize( nbNodes );
+ uniqueNodes.resize( nbNodes );
+ iRepl.resize( nbNodes );
+ int iUnique = 0, iCur = 0, nbRepl = 0;
+
+ // get new seq of nodes
+ SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+ while ( itN->more() )
+ {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( itN->next() );
+
+ TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
+ if ( nnIt != nodeNodeMap.end() ) { // n sticks
+ n = (*nnIt).second;
+ { ////////// BUG 0020185: begin
+ bool stopRecur = false;
+ set<const SMDS_MeshNode*> nodesRecur;
+ nodesRecur.insert(n);
+ while (!stopRecur) {
+ TNodeNodeMap::iterator nnIt_i = nodeNodeMap.find( n );
+ if ( nnIt_i != nodeNodeMap.end() ) { // n sticks
+ n = (*nnIt_i).second;
+ if (!nodesRecur.insert(n).second) {
+ // error: recursive dependancy
+ stopRecur = true;
+ }
+ }
+ else
+ stopRecur = true;
+ }
+ } ////////// BUG 0020185: end
+ }
+ curNodes[ iCur ] = n;
+ bool isUnique = nodeSet.insert( n ).second;
+ if ( isUnique )
+ uniqueNodes[ iUnique++ ] = n;
+ else
+ iRepl[ nbRepl++ ] = iCur;
+ iCur++;
+ }
+
+ // Analyse element topology after replacement
+
+ bool isOk = true;
+ int nbUniqueNodes = nodeSet.size();
+ if ( nbNodes != nbUniqueNodes ) // some nodes stick
+ {
+ if (elem->IsPoly()) // Polygons and Polyhedral volumes
+ {
+ if (elem->GetType() == SMDSAbs_Face) // Polygon
+ {
+ elemType.Init( elem );
+ const bool isQuad = elemType.myIsQuad;
+ if ( isQuad )
+ SMDS_MeshCell::applyInterlace // interlace medium and corner nodes
+ ( SMDS_MeshCell::interlacedSmdsOrder( SMDSEntity_Quad_Polygon, nbNodes ), curNodes );
+
+ // a polygon can divide into several elements
+ vector<const SMDS_MeshNode *> polygons_nodes;
+ vector<int> quantities;
+ int nbNew = SimplifyFace( curNodes, polygons_nodes, quantities );
+ if (nbNew > 0)
+ {
+ vector<const SMDS_MeshNode *> face_nodes;
+ int inode = 0;
+ for (int iface = 0; iface < nbNew; iface++)
+ {
+ int nbNewNodes = quantities[iface];
+ face_nodes.assign( polygons_nodes.begin() + inode,
+ polygons_nodes.begin() + inode + nbNewNodes );
+ inode += nbNewNodes;
+ if ( isQuad ) // check if a result elem is a valid quadratic polygon
+ {
+ bool isValid = ( nbNewNodes % 2 == 0 );
+ for ( int i = 0; i < nbNewNodes && isValid; ++i )
+ isValid = ( elem->IsMediumNode( face_nodes[i]) == bool( i % 2 ));
+ elemType.SetQuad( isValid );
+ if ( isValid ) // put medium nodes after corners
+ SMDS_MeshCell::applyInterlaceRev
+ ( SMDS_MeshCell::interlacedSmdsOrder( SMDSEntity_Quad_Polygon,
+ nbNewNodes ), face_nodes );
+ }
+ elemType.SetPoly(( nbNewNodes / ( elemType.myIsQuad + 1 ) > 4 ));
+
+ SMDS_MeshElement* newElem = AddElement( face_nodes, elemType );
+ if ( aShapeId )
+ aMesh->SetMeshElementOnShape(newElem, aShapeId);
+ }
+ }
+ rmElemIds.push_back(elem->GetID());
+
+ } // Polygon
+
+ else if (elem->GetType() == SMDSAbs_Volume) // Polyhedral volume
+ {
+ if (nbUniqueNodes < 4) {
+ rmElemIds.push_back(elem->GetID());
+ }
+ else {
+ // each face has to be analyzed in order to check volume validity
+ const SMDS_VtkVolume* aPolyedre = dynamic_cast<const SMDS_VtkVolume*>( elem );
+ if (aPolyedre)
+ {
+ int nbFaces = aPolyedre->NbFaces();
+
+ vector<const SMDS_MeshNode *> poly_nodes;
+ vector<int> quantities;
+
+ for (int iface = 1; iface <= nbFaces; iface++) {
+ int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
+ vector<const SMDS_MeshNode *> faceNodes (nbFaceNodes);
+
+ for (int inode = 1; inode <= nbFaceNodes; inode++) {
+ const SMDS_MeshNode * faceNode = aPolyedre->GetFaceNode(iface, inode);
+ TNodeNodeMap::iterator nnIt = nodeNodeMap.find(faceNode);
+ if (nnIt != nodeNodeMap.end()) { // faceNode sticks
+ faceNode = (*nnIt).second;
+ }
+ faceNodes[inode - 1] = faceNode;
+ }
+
+ SimplifyFace(faceNodes, poly_nodes, quantities);
+ }
+
+ if (quantities.size() > 3) {
+ // to be done: remove coincident faces
+ }
+
+ if (quantities.size() > 3)
+ {
+ const SMDS_MeshElement* newElem =
+ aMesh->AddPolyhedralVolume(poly_nodes, quantities);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId && newElem )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ rmElemIds.push_back(elem->GetID());
+ }
+ }
+ else {
+ rmElemIds.push_back(elem->GetID());
+ }
+ }
+ }
+ else {
+ }
+
+ continue;
+ } // poly element
+
+ // Regular elements
+ // TODO not all the possible cases are solved. Find something more generic?
+ switch ( nbNodes ) {
+ case 2: ///////////////////////////////////// EDGE
+ isOk = false; break;
+ case 3: ///////////////////////////////////// TRIANGLE
+ isOk = false; break;
+ case 4:
+ if ( elem->GetType() == SMDSAbs_Volume ) // TETRAHEDRON
+ isOk = false;
+ else { //////////////////////////////////// QUADRANGLE
+ if ( nbUniqueNodes < 3 )
+ isOk = false;
+ else if ( nbRepl == 2 && iRepl[ 1 ] - iRepl[ 0 ] == 2 )
+ isOk = false; // opposite nodes stick
+ //MESSAGE("isOk " << isOk);
+ }
+ break;
+ case 6: ///////////////////////////////////// PENTAHEDRON
+ if ( nbUniqueNodes == 4 ) {
+ // ---------------------------------> tetrahedron
+ if (nbRepl == 3 &&
+ iRepl[ 0 ] > 2 && iRepl[ 1 ] > 2 && iRepl[ 2 ] > 2 ) {
+ // all top nodes stick: reverse a bottom
+ uniqueNodes[ 0 ] = curNodes [ 1 ];
+ uniqueNodes[ 1 ] = curNodes [ 0 ];
+ }
+ else if (nbRepl == 3 &&
+ iRepl[ 0 ] < 3 && iRepl[ 1 ] < 3 && iRepl[ 2 ] < 3 ) {
+ // all bottom nodes stick: set a top before
+ uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
+ uniqueNodes[ 0 ] = curNodes [ 3 ];
+ uniqueNodes[ 1 ] = curNodes [ 4 ];
+ uniqueNodes[ 2 ] = curNodes [ 5 ];
+ }
+ else if (nbRepl == 4 &&
+ iRepl[ 2 ] - iRepl [ 0 ] == 3 && iRepl[ 3 ] - iRepl [ 1 ] == 3 ) {
+ // a lateral face turns into a line: reverse a bottom
+ uniqueNodes[ 0 ] = curNodes [ 1 ];
+ uniqueNodes[ 1 ] = curNodes [ 0 ];
+ }
+ else
+ isOk = false;
+ }
+ else if ( nbUniqueNodes == 5 ) {
+ // PENTAHEDRON --------------------> 2 tetrahedrons
+ if ( nbRepl == 2 && iRepl[ 1 ] - iRepl [ 0 ] == 3 ) {
+ // a bottom node sticks with a linked top one
+ // 1.
+ SMDS_MeshElement* newElem =
+ aMesh->AddVolume(curNodes[ 3 ],
+ curNodes[ 4 ],
+ curNodes[ 5 ],
+ curNodes[ iRepl[ 0 ] == 2 ? 1 : 2 ]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ // 2. : reverse a bottom
+ uniqueNodes[ 0 ] = curNodes [ 1 ];
+ uniqueNodes[ 1 ] = curNodes [ 0 ];
+ nbUniqueNodes = 4;
+ }
+ else
+ isOk = false;
+ }
+ else
+ isOk = false;
+ break;
+ case 8: {
+ if(elem->IsQuadratic()) { // Quadratic quadrangle
+ // 1 5 2
+ // +---+---+
+ // | |
+ // | |
+ // 4+ +6
+ // | |
+ // | |
+ // +---+---+
+ // 0 7 3
+ isOk = false;
+ if(nbRepl==2) {
+ MESSAGE("nbRepl=2: " << iRepl[0] << " " << iRepl[1]);
+ }
+ if(nbRepl==3) {
+ MESSAGE("nbRepl=3: " << iRepl[0] << " " << iRepl[1] << " " << iRepl[2]);
+ nbUniqueNodes = 6;
+ if( iRepl[0]==0 && iRepl[1]==1 && iRepl[2]==4 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[2];
+ uniqueNodes[2] = curNodes[3];
+ uniqueNodes[3] = curNodes[5];
+ uniqueNodes[4] = curNodes[6];
+ uniqueNodes[5] = curNodes[7];
+ isOk = true;
+ }
+ if( iRepl[0]==0 && iRepl[1]==3 && iRepl[2]==7 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[1];
+ uniqueNodes[2] = curNodes[2];
+ uniqueNodes[3] = curNodes[4];
+ uniqueNodes[4] = curNodes[5];
+ uniqueNodes[5] = curNodes[6];
+ isOk = true;
+ }
+ if( iRepl[0]==0 && iRepl[1]==4 && iRepl[2]==7 ) {
+ uniqueNodes[0] = curNodes[1];
+ uniqueNodes[1] = curNodes[2];
+ uniqueNodes[2] = curNodes[3];
+ uniqueNodes[3] = curNodes[5];
+ uniqueNodes[4] = curNodes[6];
+ uniqueNodes[5] = curNodes[0];
+ isOk = true;
+ }
+ if( iRepl[0]==1 && iRepl[1]==2 && iRepl[2]==5 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[1];
+ uniqueNodes[2] = curNodes[3];
+ uniqueNodes[3] = curNodes[4];
+ uniqueNodes[4] = curNodes[6];
+ uniqueNodes[5] = curNodes[7];
+ isOk = true;
+ }
+ if( iRepl[0]==1 && iRepl[1]==4 && iRepl[2]==5 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[2];
+ uniqueNodes[2] = curNodes[3];
+ uniqueNodes[3] = curNodes[1];
+ uniqueNodes[4] = curNodes[6];
+ uniqueNodes[5] = curNodes[7];
+ isOk = true;
+ }
+ if( iRepl[0]==2 && iRepl[1]==3 && iRepl[2]==6 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[1];
+ uniqueNodes[2] = curNodes[2];
+ uniqueNodes[3] = curNodes[4];
+ uniqueNodes[4] = curNodes[5];
+ uniqueNodes[5] = curNodes[7];
+ isOk = true;
+ }
+ if( iRepl[0]==2 && iRepl[1]==5 && iRepl[2]==6 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[1];
+ uniqueNodes[2] = curNodes[3];
+ uniqueNodes[3] = curNodes[4];
+ uniqueNodes[4] = curNodes[2];
+ uniqueNodes[5] = curNodes[7];
+ isOk = true;
+ }
+ if( iRepl[0]==3 && iRepl[1]==6 && iRepl[2]==7 ) {
+ uniqueNodes[0] = curNodes[0];
+ uniqueNodes[1] = curNodes[1];
+ uniqueNodes[2] = curNodes[2];
+ uniqueNodes[3] = curNodes[4];
+ uniqueNodes[4] = curNodes[5];
+ uniqueNodes[5] = curNodes[3];
+ isOk = true;
+ }
+ }
+ if(nbRepl==4) {
+ MESSAGE("nbRepl=4: " << iRepl[0] << " " << iRepl[1] << " " << iRepl[2] << " " << iRepl[3]);
+ }
+ if(nbRepl==5) {
+ MESSAGE("nbRepl=5: " << iRepl[0] << " " << iRepl[1] << " " << iRepl[2] << " " << iRepl[3] << " " << iRepl[4]);
+ }
+ break;
+ }
+ //////////////////////////////////// HEXAHEDRON
+ isOk = false;
+ SMDS_VolumeTool hexa (elem);
+ hexa.SetExternalNormal();
+ if ( nbUniqueNodes == 4 && nbRepl == 4 ) {
+ //////////////////////// HEX ---> 1 tetrahedron
+ for ( int iFace = 0; iFace < 6; iFace++ ) {
+ const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
+ if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
+ curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
+ curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
+ // one face turns into a point ...
+ int iOppFace = hexa.GetOppFaceIndex( iFace );
+ ind = hexa.GetFaceNodesIndices( iOppFace );
+ int nbStick = 0;
+ for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
+ if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
+ nbStick++;
+ }
+ if ( nbStick == 1 ) {
+ // ... and the opposite one - into a triangle.
+ // set a top node
+ ind = hexa.GetFaceNodesIndices( iFace );
+ uniqueNodes[ 3 ] = curNodes[ind[ 0 ]];
+ isOk = true;
+ }
+ break;
+ }
+ }
+ }
+ else if ( nbUniqueNodes == 6 && nbRepl == 2 ) {
+ //////////////////////// HEX ---> 1 prism
+ int nbTria = 0, iTria[3];
+ const int *ind; // indices of face nodes
+ // look for triangular faces
+ for ( int iFace = 0; iFace < 6 && nbTria < 3; iFace++ ) {
+ ind = hexa.GetFaceNodesIndices( iFace );
+ TIDSortedNodeSet faceNodes;
+ for ( iCur = 0; iCur < 4; iCur++ )
+ faceNodes.insert( curNodes[ind[iCur]] );
+ if ( faceNodes.size() == 3 )
+ iTria[ nbTria++ ] = iFace;
+ }
+ // check if triangles are opposite
+ if ( nbTria == 2 && iTria[0] == hexa.GetOppFaceIndex( iTria[1] ))
+ {
+ isOk = true;
+ // set nodes of the bottom triangle
+ ind = hexa.GetFaceNodesIndices( iTria[ 0 ]);
+ vector<int> indB;
+ for ( iCur = 0; iCur < 4; iCur++ )
+ if ( ind[iCur] != iRepl[0] && ind[iCur] != iRepl[1])
+ indB.push_back( ind[iCur] );
+ if ( !hexa.IsForward() )
+ std::swap( indB[0], indB[2] );
+ for ( iCur = 0; iCur < 3; iCur++ )
+ uniqueNodes[ iCur ] = curNodes[indB[iCur]];
+ // set nodes of the top triangle
+ const int *indT = hexa.GetFaceNodesIndices( iTria[ 1 ]);
+ for ( iCur = 0; iCur < 3; ++iCur )
+ for ( int j = 0; j < 4; ++j )
+ if ( hexa.IsLinked( indB[ iCur ], indT[ j ] ))
+ {
+ uniqueNodes[ iCur + 3 ] = curNodes[ indT[ j ]];
+ break;
+ }
+ }
+ break;
+ }
+ else if (nbUniqueNodes == 5 && nbRepl == 4 ) {
+ //////////////////// HEXAHEDRON ---> 2 tetrahedrons
+ for ( int iFace = 0; iFace < 6; iFace++ ) {
+ const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
+ if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
+ curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
+ curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
+ // one face turns into a point ...
+ int iOppFace = hexa.GetOppFaceIndex( iFace );
+ ind = hexa.GetFaceNodesIndices( iOppFace );
+ int nbStick = 0;
+ iUnique = 2; // reverse a tetrahedron 1 bottom
+ for ( iCur = 0; iCur < 4 && nbStick == 0; iCur++ ) {
+ if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
+ nbStick++;
+ else if ( iUnique >= 0 )
+ uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
+ }
+ if ( nbStick == 0 ) {
+ // ... and the opposite one is a quadrangle
+ // set a top node
+ const int* indTop = hexa.GetFaceNodesIndices( iFace );
+ uniqueNodes[ 3 ] = curNodes[indTop[ 0 ]];
+ nbUniqueNodes = 4;
+ // tetrahedron 2
+ SMDS_MeshElement* newElem =
+ aMesh->AddVolume(curNodes[ind[ 0 ]],
+ curNodes[ind[ 3 ]],
+ curNodes[ind[ 2 ]],
+ curNodes[indTop[ 0 ]]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ isOk = true;
+ }
+ break;
+ }
+ }
+ }
+ else if ( nbUniqueNodes == 6 && nbRepl == 4 ) {
+ ////////////////// HEXAHEDRON ---> 2 tetrahedrons or 1 prism
+ // find indices of quad and tri faces
+ int iQuadFace[ 6 ], iTriFace[ 6 ], nbQuad = 0, nbTri = 0, iFace;
+ for ( iFace = 0; iFace < 6; iFace++ ) {
+ const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
+ nodeSet.clear();
+ for ( iCur = 0; iCur < 4; iCur++ )
+ nodeSet.insert( curNodes[ind[ iCur ]] );
+ nbUniqueNodes = nodeSet.size();
+ if ( nbUniqueNodes == 3 )
+ iTriFace[ nbTri++ ] = iFace;
+ else if ( nbUniqueNodes == 4 )
+ iQuadFace[ nbQuad++ ] = iFace;
+ }
+ if (nbQuad == 2 && nbTri == 4 &&
+ hexa.GetOppFaceIndex( iQuadFace[ 0 ] ) == iQuadFace[ 1 ]) {
+ // 2 opposite quadrangles stuck with a diagonal;
+ // sample groups of merged indices: (0-4)(2-6)
+ // --------------------------------------------> 2 tetrahedrons
+ const int *ind1 = hexa.GetFaceNodesIndices( iQuadFace[ 0 ]); // indices of quad1 nodes
+ const int *ind2 = hexa.GetFaceNodesIndices( iQuadFace[ 1 ]);
+ int i0, i1d, i2, i3d, i0t, i2t; // d-daigonal, t-top
+ if (curNodes[ind1[ 0 ]] == curNodes[ind2[ 0 ]] &&
+ curNodes[ind1[ 2 ]] == curNodes[ind2[ 2 ]]) {
+ // stuck with 0-2 diagonal
+ i0 = ind1[ 3 ];
+ i1d = ind1[ 0 ];
+ i2 = ind1[ 1 ];
+ i3d = ind1[ 2 ];
+ i0t = ind2[ 1 ];
+ i2t = ind2[ 3 ];
+ }
+ else if (curNodes[ind1[ 1 ]] == curNodes[ind2[ 3 ]] &&
+ curNodes[ind1[ 3 ]] == curNodes[ind2[ 1 ]]) {
+ // stuck with 1-3 diagonal
+ i0 = ind1[ 0 ];
+ i1d = ind1[ 1 ];
+ i2 = ind1[ 2 ];
+ i3d = ind1[ 3 ];
+ i0t = ind2[ 0 ];
+ i2t = ind2[ 1 ];
+ }
+ else {
+ ASSERT(0);
+ }
+ // tetrahedron 1
+ uniqueNodes[ 0 ] = curNodes [ i0 ];
+ uniqueNodes[ 1 ] = curNodes [ i1d ];
+ uniqueNodes[ 2 ] = curNodes [ i3d ];
+ uniqueNodes[ 3 ] = curNodes [ i0t ];
+ nbUniqueNodes = 4;
+ // tetrahedron 2
+ SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ i1d ],
+ curNodes[ i2 ],
+ curNodes[ i3d ],
+ curNodes[ i2t ]);
+ myLastCreatedElems.Append(newElem);
+ if ( aShapeId )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ isOk = true;
+ }
+ else if (( nbTri == 2 && nbQuad == 3 ) || // merged (0-4)(1-5)
+ ( nbTri == 4 && nbQuad == 2 )) { // merged (7-4)(1-5)
+ // --------------------------------------------> prism
+ // find 2 opposite triangles
+ nbUniqueNodes = 6;
+ for ( iFace = 0; iFace + 1 < nbTri; iFace++ ) {
+ if ( hexa.GetOppFaceIndex( iTriFace[ iFace ] ) == iTriFace[ iFace + 1 ]) {
+ // find indices of kept and replaced nodes
+ // and fill unique nodes of 2 opposite triangles
+ const int *ind1 = hexa.GetFaceNodesIndices( iTriFace[ iFace ]);
+ const int *ind2 = hexa.GetFaceNodesIndices( iTriFace[ iFace + 1 ]);
+ const SMDS_MeshNode** hexanodes = hexa.GetNodes();
+ // fill unique nodes
+ iUnique = 0;
+ isOk = true;
+ for ( iCur = 0; iCur < 4 && isOk; iCur++ ) {
+ const SMDS_MeshNode* n = curNodes[ind1[ iCur ]];
+ const SMDS_MeshNode* nInit = hexanodes[ind1[ iCur ]];
+ if ( n == nInit ) {
+ // iCur of a linked node of the opposite face (make normals co-directed):
+ int iCurOpp = ( iCur == 1 || iCur == 3 ) ? 4 - iCur : iCur;
+ // check that correspondent corners of triangles are linked
+ if ( !hexa.IsLinked( ind1[ iCur ], ind2[ iCurOpp ] ))
+ isOk = false;
+ else {
+ uniqueNodes[ iUnique ] = n;
+ uniqueNodes[ iUnique + 3 ] = curNodes[ind2[ iCurOpp ]];
+ iUnique++;
+ }
+ }
+ }
+ break;
+ }
+ }
+ }
+ } // if ( nbUniqueNodes == 6 && nbRepl == 4 )
+ else
+ {
+ MESSAGE("MergeNodes() removes hexahedron "<< elem);
+ }
+ break;
+ } // HEXAHEDRON
+
+ default:
+ isOk = false;
+ } // switch ( nbNodes )
+
+ } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
+
+ if ( isOk ) // the non-poly elem remains valid after sticking nodes
+ {
+ if ( nbNodes != nbUniqueNodes ||
+ !aMesh->ChangeElementNodes( elem, & curNodes[0], nbNodes ))
+ {
+ elemType.Init( elem ).SetID( elem->GetID() );
+
+ SMESHDS_SubMesh * sm = aShapeId > 0 ? aMesh->MeshElements(aShapeId) : 0;
+ aMesh->RemoveFreeElement(elem, sm, /*fromGroups=*/false);
+
+ uniqueNodes.resize(nbUniqueNodes);
+ SMDS_MeshElement* newElem = this->AddElement( uniqueNodes, elemType );
+ if ( sm && newElem )
+ sm->AddElement( newElem );
+ if ( elem != newElem )
+ ReplaceElemInGroups( elem, newElem, aMesh );
+ }
+ }
+ else {
+ // Remove invalid regular element or invalid polygon
+ rmElemIds.push_back( elem->GetID() );
+ }
+
+ } // loop on elements
+
+ // Remove bad elements, then equal nodes (order important)
+
+ Remove( rmElemIds, false );
+ Remove( rmNodeIds, true );
+
+ return;
+}
+
+
+// ========================================================
+// class : SortableElement
+// purpose : allow sorting elements basing on their nodes
+// ========================================================
+class SortableElement : public set <const SMDS_MeshElement*>
+{
+public:
+
+ SortableElement( const SMDS_MeshElement* theElem )
+ {
+ myElem = theElem;
+ SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
+ while ( nodeIt->more() )
+ this->insert( nodeIt->next() );
+ }
+
+ const SMDS_MeshElement* Get() const
+ { return myElem; }
+
+private:
+ mutable const SMDS_MeshElement* myElem;
+};
+
+//=======================================================================
+//function : FindEqualElements
+//purpose : Return list of group of elements built on the same nodes.
+// Search among theElements or in the whole mesh if theElements is empty
+//=======================================================================
+
+void SMESH_MeshEditor::FindEqualElements(TIDSortedElemSet & theElements,
+ TListOfListOfElementsID & theGroupsOfElementsID)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ typedef map< SortableElement, int > TMapOfNodeSet;
+ typedef list<int> TGroupOfElems;
+
+ if ( theElements.empty() )
+ { // get all elements in the mesh
+ SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
+ while ( eIt->more() )
+ theElements.insert( theElements.end(), eIt->next() );
+ }
+
+ vector< TGroupOfElems > arrayOfGroups;
+ TGroupOfElems groupOfElems;
+ TMapOfNodeSet mapOfNodeSet;
+
+ TIDSortedElemSet::iterator elemIt = theElements.begin();
+ for ( int i = 0; elemIt != theElements.end(); ++elemIt )
+ {
+ const SMDS_MeshElement* curElem = *elemIt;
+ SortableElement SE(curElem);
+ // check uniqueness
+ pair< TMapOfNodeSet::iterator, bool> pp = mapOfNodeSet.insert(make_pair(SE, i));
+ if ( !pp.second ) { // one more coincident elem
+ TMapOfNodeSet::iterator& itSE = pp.first;
+ int ind = (*itSE).second;
+ arrayOfGroups[ind].push_back( curElem->GetID() );
+ }
+ else {
+ arrayOfGroups.push_back( groupOfElems );
+ arrayOfGroups.back().push_back( curElem->GetID() );
+ i++;
+ }
+ }
+
+ groupOfElems.clear();
+ vector< TGroupOfElems >::iterator groupIt = arrayOfGroups.begin();
+ for ( ; groupIt != arrayOfGroups.end(); ++groupIt )
+ {
+ if ( groupIt->size() > 1 ) {
+ //groupOfElems.sort(); -- theElements is sorted already
+ theGroupsOfElementsID.push_back( groupOfElems );
+ theGroupsOfElementsID.back().splice( theGroupsOfElementsID.back().end(), *groupIt );
+ }
+ }
+}
+
+//=======================================================================
+//function : MergeElements
+//purpose : In each given group, substitute all elements by the first one.
+//=======================================================================
+
+void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElementsID)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ typedef list<int> TListOfIDs;
+ TListOfIDs rmElemIds; // IDs of elems to remove
+
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+
+ TListOfListOfElementsID::iterator groupsIt = theGroupsOfElementsID.begin();
+ while ( groupsIt != theGroupsOfElementsID.end() ) {
+ TListOfIDs& aGroupOfElemID = *groupsIt;
+ aGroupOfElemID.sort();
+ int elemIDToKeep = aGroupOfElemID.front();
+ const SMDS_MeshElement* elemToKeep = aMesh->FindElement(elemIDToKeep);
+ aGroupOfElemID.pop_front();
+ TListOfIDs::iterator idIt = aGroupOfElemID.begin();
+ while ( idIt != aGroupOfElemID.end() ) {
+ int elemIDToRemove = *idIt;
+ const SMDS_MeshElement* elemToRemove = aMesh->FindElement(elemIDToRemove);
+ // add the kept element in groups of removed one (PAL15188)
+ AddToSameGroups( elemToKeep, elemToRemove, aMesh );
+ rmElemIds.push_back( elemIDToRemove );
+ ++idIt;
+ }
+ ++groupsIt;
+ }
+
+ Remove( rmElemIds, false );
+}
+
+//=======================================================================
+//function : MergeEqualElements
+//purpose : Remove all but one of elements built on the same nodes.
+//=======================================================================
+
+void SMESH_MeshEditor::MergeEqualElements()
+{
+ TIDSortedElemSet aMeshElements; /* empty input ==
+ to merge equal elements in the whole mesh */
+ TListOfListOfElementsID aGroupsOfElementsID;
+ FindEqualElements(aMeshElements, aGroupsOfElementsID);
+ MergeElements(aGroupsOfElementsID);
+}
+
+//=======================================================================
+//function : findAdjacentFace
+//purpose :
+//=======================================================================
+
+static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ const SMDS_MeshElement* elem)
+{
+ TIDSortedElemSet elemSet, avoidSet;
+ if ( elem )
+ avoidSet.insert ( elem );
+ return SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet );
+}
+
+//=======================================================================
+//function : findSegment
+//purpose : Return a mesh segment by two nodes one of which can be medium
+//=======================================================================
+
+static const SMDS_MeshElement* findSegment(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2)
+{
+ SMDS_ElemIteratorPtr it = n1->GetInverseElementIterator( SMDSAbs_Edge );
+ while ( it->more() )
+ {
+ const SMDS_MeshElement* seg = it->next();
+ if ( seg->GetNodeIndex( n2 ) >= 0 )
+ return seg;
+ }
+ return 0;
+}
+
+//=======================================================================
+//function : FindFreeBorder
+//purpose :
+//=======================================================================
+
+#define ControlFreeBorder SMESH::Controls::FreeEdges::IsFreeEdge
+
+bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode* theFirstNode,
+ const SMDS_MeshNode* theSecondNode,
+ const SMDS_MeshNode* theLastNode,
+ list< const SMDS_MeshNode* > & theNodes,
+ list< const SMDS_MeshElement* >& theFaces)
+{
+ if ( !theFirstNode || !theSecondNode )
+ return false;
+ // find border face between theFirstNode and theSecondNode
+ const SMDS_MeshElement* curElem = findAdjacentFace( theFirstNode, theSecondNode, 0 );
+ if ( !curElem )
+ return false;
+
+ theFaces.push_back( curElem );
+ theNodes.push_back( theFirstNode );
+ theNodes.push_back( theSecondNode );
+
+ const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
+ TIDSortedElemSet foundElems;
+ bool needTheLast = ( theLastNode != 0 );
+
+ 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();
+ vector<const SMDS_MeshNode*> nodes(nbNodes+1);
+
+ if ( e->IsQuadratic() ) {
+ const SMDS_VtkFace* F =
+ dynamic_cast<const SMDS_VtkFace*>(e);
+ if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
+ // use special nodes iterator
+ SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
+ while( anIter->more() ) {
+ nodes[ iNode++ ] = cast2Node(anIter->next());
+ }
+ }
+ else {
+ SMDS_ElemIteratorPtr nIt = e->nodesIterator();
+ while ( nIt->more() )
+ nodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
+ }
+ nodes[ iNode ] = nodes[ 0 ];
+ // check 2 links
+ for ( iNode = 0; iNode < nbNodes; iNode++ )
+ if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
+ (nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
+ ControlFreeBorder( &nodes[ iNode ], e->GetID() ))
+ {
+ nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart ? 1 : 0 )]);
+ curElemList.push_back( e );
+ }
+ }
+ }
+ // analyse the found
+
+ int nbNewBorders = curElemList.size();
+ if ( nbNewBorders == 0 ) {
+ // no free border furthermore
+ return !needTheLast;
+ }
+ else if ( nbNewBorders == 1 ) {
+ // one more element found
+ nIgnore = nStart;
+ nStart = nStartList.front();
+ curElem = curElemList.front();
+ theFaces.push_back( curElem );
+ theNodes.push_back( nStart );
+ }
+ else {
+ // several continuations found
+ list< const SMDS_MeshElement* >::iterator curElemIt;
+ list< const SMDS_MeshNode* >::iterator nStartIt;
+ // check if one of them reached the last node
+ if ( needTheLast ) {
+ for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
+ curElemIt!= curElemList.end();
+ curElemIt++, nStartIt++ )
+ if ( *nStartIt == theLastNode ) {
+ theFaces.push_back( *curElemIt );
+ theNodes.push_back( *nStartIt );
+ return true;
+ }
+ }
+ // find the best free border by the continuations
+ list<const SMDS_MeshNode*> contNodes[ 2 ], *cNL;
+ list<const SMDS_MeshElement*> contFaces[ 2 ], *cFL;
+ for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
+ curElemIt!= curElemList.end();
+ curElemIt++, nStartIt++ )
+ {
+ cNL = & contNodes[ contNodes[0].empty() ? 0 : 1 ];
+ cFL = & contFaces[ contFaces[0].empty() ? 0 : 1 ];
+ // find one more free border
+ if ( ! SMESH_MeshEditor::FindFreeBorder( nStart, *nStartIt, theLastNode, *cNL, *cFL )) {
+ cNL->clear();
+ cFL->clear();
+ }
+ else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
+ // choice: clear a worse one
+ int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
+ int iWorse = ( needTheLast ? 1 - iLongest : iLongest );
+ contNodes[ iWorse ].clear();
+ contFaces[ iWorse ].clear();
+ }
+ }
+ if ( contNodes[0].empty() && contNodes[1].empty() )
+ return false;
+
+ // append the best free border
+ cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
+ cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
+ theNodes.pop_back(); // remove nIgnore
+ theNodes.pop_back(); // remove nStart
+ theFaces.pop_back(); // remove curElem
+ list< const SMDS_MeshNode* >::iterator nIt = cNL->begin();
+ list< const SMDS_MeshElement* >::iterator fIt = cFL->begin();
+ for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt );
+ for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt );
+ return true;
+
+ } // several continuations found
+ } // while ( nStart != theLastNode )
+
+ return true;
+}
+
+//=======================================================================
+//function : CheckFreeBorderNodes
+//purpose : Return true if the tree nodes are on a free border
+//=======================================================================
+
+bool SMESH_MeshEditor::CheckFreeBorderNodes(const SMDS_MeshNode* theNode1,
+ const SMDS_MeshNode* theNode2,
+ const SMDS_MeshNode* theNode3)
+{
+ list< const SMDS_MeshNode* > nodes;
+ list< const SMDS_MeshElement* > faces;
+ return FindFreeBorder( theNode1, theNode2, theNode3, nodes, faces);
+}
+
+//=======================================================================
+//function : SewFreeBorder
+//purpose :
+//warning : for border-to-side sewing theSideSecondNode is considered as
+// the last side node and theSideThirdNode is not used
+//=======================================================================
+
+SMESH_MeshEditor::Sew_Error
+SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
+ const SMDS_MeshNode* theBordSecondNode,
+ const SMDS_MeshNode* theBordLastNode,
+ const SMDS_MeshNode* theSideFirstNode,
+ const SMDS_MeshNode* theSideSecondNode,
+ const SMDS_MeshNode* theSideThirdNode,
+ const bool theSideIsFreeBorder,
+ const bool toCreatePolygons,
+ const bool toCreatePolyedrs)
+{
+ 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[] = { (int)nSide[0].size(), (int)nSide[1].size() };
+ int maxNbNodes = Max( nbNodes[0], nbNodes[1] );
+
+ bool toMergeConformal = ( nbNodes[0] == nbNodes[1] );
+ if ( toMergeConformal && toCreatePolygons )
+ {
+ // do not merge quadrangles if polygons are OK (IPAL0052824)
+ eIt[0] = eSide[0].begin();
+ eIt[1] = eSide[1].begin();
+ bool allQuads[2] = { true, true };
+ for ( int iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
+ for ( ; allQuads[iBord] && eIt[iBord] != eSide[iBord].end(); ++eIt[iBord] )
+ allQuads[iBord] = ( (*eIt[iBord])->NbCornerNodes() == 4 );
+ }
+ toMergeConformal = ( !allQuads[0] && !allQuads[1] );
+ }
+
+ TListOfListOfNodes nodeGroupsToMerge;
+ if (( toMergeConformal ) ||
+ ( theSideIsFreeBorder && !theSideThirdNode )) {
+
+ // all nodes are to be merged
+
+ for (nIt[0] = nSide[0].begin(), nIt[1] = nSide[1].begin();
+ nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end();
+ nIt[0]++, nIt[1]++ )
+ {
+ nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
+ nodeGroupsToMerge.back().push_back( *nIt[1] ); // to keep
+ nodeGroupsToMerge.back().push_back( *nIt[0] ); // to remove
+ }
+ }
+ else {
+
+ // insert new nodes into the border and the side to get equal nb of segments
+
+ // get normalized parameters of nodes on the borders
+ vector< double > param[ 2 ];
+ param[0].resize( maxNbNodes );
+ param[1].resize( maxNbNodes );
+ int iNode, iBord;
+ for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
+ list< const SMDS_MeshNode* >& nodes = nSide[ iBord ];
+ list< const SMDS_MeshNode* >::iterator nIt = nodes.begin();
+ const SMDS_MeshNode* nPrev = *nIt;
+ double bordLength = 0;
+ for ( iNode = 0; nIt != nodes.end(); nIt++, iNode++ ) { // loop on border nodes
+ const SMDS_MeshNode* nCur = *nIt;
+ gp_XYZ segment (nCur->X() - nPrev->X(),
+ nCur->Y() - nPrev->Y(),
+ nCur->Z() - nPrev->Z());
+ double segmentLen = segment.Modulus();
+ bordLength += segmentLen;
+ param[ iBord ][ iNode ] = bordLength;
+ nPrev = nCur;
+ }
+ // normalize within [0,1]
+ for ( iNode = 0; iNode < nbNodes[ iBord ]; iNode++ ) {
+ param[ iBord ][ iNode ] /= bordLength;
+ }
+ }
+
+ // loop on border segments
+ const SMDS_MeshNode *nPrev[ 2 ] = { 0, 0 };
+ int i[ 2 ] = { 0, 0 };
+ nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
+ nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
+
+ 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 ( const SMDS_MeshElement* adjElem = findAdjacentFace( n12, n22, elem )) {
+ InsertNodesIntoLink( adjElem, n12, n22, nodeList, toCreatePolygons );
+ }
+ while ( const SMDS_MeshElement* seg = findSegment( n12, n22 )) {
+ InsertNodesIntoLink( seg, n12, n22, nodeList );
+ }
+ if (toCreatePolyedrs) {
+ // perform insertion into the links of adjacent volumes
+ UpdateVolumes(n12, n22, nodeList);
+ }
+ // 3. find an element appeared on n1 and n2 after the insertion
+ insertMap.erase( elem );
+ elem = findAdjacentFace( n1, n2, 0 );
+ }
+ if ( notFound || otherLink ) {
+ // add element and nodes of the side into the insertMap
+ insertMapIt = insertMap.insert( make_pair( elem, list<const SMDS_MeshNode*>() )).first;
+ (*insertMapIt).second.push_back( n1 );
+ (*insertMapIt).second.push_back( n2 );
+ }
+ // add node to be inserted into elem
+ (*insertMapIt).second.push_back( nIns );
+ next[ 1 - intoBord ] = true;
+ }
+
+ // go to the next segment
+ for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
+ if ( next[ iBord ] ) {
+ if ( i[ iBord ] != 0 && eIt[ iBord ] != eSide[ iBord ].end())
+ eIt[ iBord ]++;
+ nPrev[ iBord ] = *nIt[ iBord ];
+ nIt[ iBord ]++; i[ iBord ]++;
+ }
+ }
+ }
+ while ( nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end());
+
+ // perform insertion of nodes into elements
+
+ for (insertMapIt = insertMap.begin();
+ insertMapIt != insertMap.end();
+ insertMapIt++ )
+ {
+ const SMDS_MeshElement* elem = (*insertMapIt).first;
+ list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
+ const SMDS_MeshNode* n1 = nodeList.front(); nodeList.pop_front();
+ const SMDS_MeshNode* n2 = nodeList.front(); nodeList.pop_front();
+
+ InsertNodesIntoLink( elem, n1, n2, nodeList, toCreatePolygons );
+
+ while ( const SMDS_MeshElement* seg = findSegment( n1, n2 )) {
+ InsertNodesIntoLink( seg, n1, n2, nodeList );
+ }
+
+ if ( !theSideIsFreeBorder ) {
+ // look for and insert nodes into the faces adjacent to elem
+ while ( const SMDS_MeshElement* adjElem = findAdjacentFace( n1, n2, elem )) {
+ InsertNodesIntoLink( adjElem, n1, n2, nodeList, toCreatePolygons );
+ }
+ }
+ if (toCreatePolyedrs) {
+ // perform insertion into the links of adjacent volumes
+ UpdateVolumes(n1, n2, nodeList);
+ }
+ }
+ } // end: insert new nodes
+
+ MergeNodes ( nodeGroupsToMerge );
+
+
+ // Remove coincident segments
+
+ // get new segments
+ TIDSortedElemSet segments;
+ SMESH_SequenceOfElemPtr newFaces;
+ for ( int i = 1; i <= myLastCreatedElems.Length(); ++i )
+ {
+ if ( !myLastCreatedElems(i) ) continue;
+ if ( myLastCreatedElems(i)->GetType() == SMDSAbs_Edge )
+ segments.insert( segments.end(), myLastCreatedElems(i) );
+ else
+ newFaces.Append( myLastCreatedElems(i) );
+ }
+ // get segments adjacent to merged nodes
+ TListOfListOfNodes::iterator groupIt = nodeGroupsToMerge.begin();
+ for ( ; groupIt != nodeGroupsToMerge.end(); groupIt++ )
+ {
+ const list<const SMDS_MeshNode*>& nodes = *groupIt;
+ SMDS_ElemIteratorPtr segIt = nodes.front()->GetInverseElementIterator( SMDSAbs_Edge );
+ while ( segIt->more() )
+ segments.insert( segIt->next() );
+ }
+
+ // find coincident
+ TListOfListOfElementsID equalGroups;
+ if ( !segments.empty() )
+ FindEqualElements( segments, equalGroups );
+ if ( !equalGroups.empty() )
+ {
+ // remove from segments those that will be removed
+ TListOfListOfElementsID::iterator itGroups = equalGroups.begin();
+ for ( ; itGroups != equalGroups.end(); ++itGroups )
+ {
+ list< int >& group = *itGroups;
+ list< int >::iterator id = group.begin();
+ for ( ++id; id != group.end(); ++id )
+ if ( const SMDS_MeshElement* seg = GetMeshDS()->FindElement( *id ))
+ segments.erase( seg );
+ }
+ // remove equal segments
+ MergeElements( equalGroups );
+
+ // restore myLastCreatedElems
+ myLastCreatedElems = newFaces;
+ TIDSortedElemSet::iterator seg = segments.begin();
+ for ( ; seg != segments.end(); ++seg )
+ myLastCreatedElems.Append( *seg );
+ }
+
+ return aResult;
+}
+
+//=======================================================================
+//function : InsertNodesIntoLink
+//purpose : insert theNodesToInsert into theElement between theBetweenNode1
+// and theBetweenNode2 and split theElement
+//=======================================================================
+
+void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theElement,
+ const SMDS_MeshNode* theBetweenNode1,
+ const SMDS_MeshNode* theBetweenNode2,
+ list<const SMDS_MeshNode*>& theNodesToInsert,
+ const bool toCreatePoly)
+{
+ if ( !theElement ) return;
+
+ SMESHDS_Mesh *aMesh = GetMeshDS();
+ vector<const SMDS_MeshElement*> newElems;
+
+ if ( theElement->GetType() == SMDSAbs_Edge )
+ {
+ theNodesToInsert.push_front( theBetweenNode1 );
+ theNodesToInsert.push_back ( theBetweenNode2 );
+ list<const SMDS_MeshNode*>::iterator n = theNodesToInsert.begin();
+ const SMDS_MeshNode* n1 = *n;
+ for ( ++n; n != theNodesToInsert.end(); ++n )
+ {
+ const SMDS_MeshNode* n2 = *n;
+ if ( const SMDS_MeshElement* seg = aMesh->FindEdge( n1, n2 ))
+ AddToSameGroups( seg, theElement, aMesh );
+ else
+ newElems.push_back( aMesh->AddEdge ( n1, n2 ));
+ n1 = n2;
+ }
+ theNodesToInsert.pop_front();
+ theNodesToInsert.pop_back();
+
+ if ( theElement->IsQuadratic() ) // add a not split part
+ {
+ vector<const SMDS_MeshNode*> nodes( theElement->begin_nodes(),
+ theElement->end_nodes() );
+ int iOther = 0, nbN = nodes.size();
+ for ( ; iOther < nbN; ++iOther )
+ if ( nodes[iOther] != theBetweenNode1 &&
+ nodes[iOther] != theBetweenNode2 )
+ break;
+ if ( iOther == 0 )
+ {
+ if ( const SMDS_MeshElement* seg = aMesh->FindEdge( nodes[0], nodes[1] ))
+ AddToSameGroups( seg, theElement, aMesh );
+ else
+ newElems.push_back( aMesh->AddEdge ( nodes[0], nodes[1] ));
+ }
+ else if ( iOther == 2 )
+ {
+ if ( const SMDS_MeshElement* seg = aMesh->FindEdge( nodes[1], nodes[2] ))
+ AddToSameGroups( seg, theElement, aMesh );
+ else
+ newElems.push_back( aMesh->AddEdge ( nodes[1], nodes[2] ));
+ }
+ }
+ // treat new elements
+ for ( size_t i = 0; i < newElems.size(); ++i )
+ if ( newElems[i] )
+ {
+ aMesh->SetMeshElementOnShape( newElems[i], theElement->getshapeId() );
+ myLastCreatedElems.Append( newElems[i] );
+ }
+ ReplaceElemInGroups( theElement, newElems, aMesh );
+ aMesh->RemoveElement( theElement );
+ return;
+
+ } // if ( theElement->GetType() == SMDSAbs_Edge )
+
+ const SMDS_MeshElement* theFace = theElement;
+ if ( theFace->GetType() != SMDSAbs_Face ) return;
+
+ // find indices of 2 link nodes and of the rest nodes
+ int iNode = 0, il1, il2, i3, i4;
+ il1 = il2 = i3 = i4 = -1;
+ vector<const SMDS_MeshNode*> nodes( theFace->NbNodes() );
+
+ SMDS_NodeIteratorPtr nodeIt = theFace->interlacedNodesIterator();
+ while ( nodeIt->more() ) {
+ const SMDS_MeshNode* n = nodeIt->next();
+ if ( n == theBetweenNode1 )
+ il1 = iNode;
+ else if ( n == theBetweenNode2 )
+ il2 = iNode;
+ else if ( i3 < 0 )
+ i3 = iNode;
+ else
+ i4 = iNode;
+ nodes[ iNode++ ] = n;
+ }
+ if ( il1 < 0 || il2 < 0 || i3 < 0 )
+ return ;
+
+ // arrange link nodes to go one after another regarding the face orientation
+ bool reverse = ( Abs( il2 - il1 ) == 1 ? il2 < il1 : il1 < il2 );
+ list<const SMDS_MeshNode *> aNodesToInsert = theNodesToInsert;
+ if ( reverse ) {
+ iNode = il1;
+ il1 = il2;
+ il2 = iNode;
+ aNodesToInsert.reverse();
+ }
+ // check that not link nodes of a quadrangles are in good order
+ int nbFaceNodes = theFace->NbNodes();
+ if ( nbFaceNodes == 4 && i4 - i3 != 1 ) {
+ iNode = i3;
+ i3 = i4;
+ i4 = iNode;
+ }
+
+ if (toCreatePoly || theFace->IsPoly()) {
+
+ iNode = 0;
+ vector<const SMDS_MeshNode *> poly_nodes (nbFaceNodes + aNodesToInsert.size());
+
+ // add nodes of face up to first node of link
+ bool isFLN = false;
+
+ 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;
+ }
+ }
+
+ // make a new face
+ newElems.push_back( aMesh->AddPolygonalFace( poly_nodes ));
+ }
+
+ else if ( !theFace->IsQuadratic() )
+ {
+ // put aNodesToInsert between theBetweenNode1 and theBetweenNode2
+ int nbLinkNodes = 2 + aNodesToInsert.size();
+ //const SMDS_MeshNode* linkNodes[ nbLinkNodes ];
+ vector<const SMDS_MeshNode*> linkNodes( nbLinkNodes );
+ linkNodes[ 0 ] = nodes[ il1 ];
+ linkNodes[ nbLinkNodes - 1 ] = nodes[ il2 ];
+ list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
+ for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
+ linkNodes[ iNode++ ] = *nIt;
+ }
+ // decide how to split a quadrangle: compare possible variants
+ // and choose which of splits to be a quadrangle
+ int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad;
+ if ( nbFaceNodes == 3 ) {
+ iBestQuad = nbSplits;
+ i4 = i3;
+ }
+ else if ( nbFaceNodes == 4 ) {
+ SMESH::Controls::NumericalFunctorPtr aCrit( new SMESH::Controls::AspectRatio);
+ double aBestRate = DBL_MAX;
+ for ( int iQuad = 0; iQuad < nbSplits; iQuad++ ) {
+ i1 = 0; i2 = 1;
+ double aBadRate = 0;
+ // evaluate elements quality
+ for ( iSplit = 0; iSplit < nbSplits; iSplit++ ) {
+ if ( iSplit == iQuad ) {
+ SMDS_FaceOfNodes quad (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ i3 ],
+ nodes[ i4 ]);
+ aBadRate += getBadRate( &quad, aCrit );
+ }
+ else {
+ SMDS_FaceOfNodes tria (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ iSplit < iQuad ? i4 : i3 ]);
+ aBadRate += getBadRate( &tria, aCrit );
+ }
+ }
+ // choice
+ if ( aBadRate < aBestRate ) {
+ iBestQuad = iQuad;
+ aBestRate = aBadRate;
+ }
+ }
+ }
+
+ // create new elements
+ i1 = 0; i2 = 1;
+ for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ )
+ {
+ if ( iSplit == iBestQuad )
+ newElems.push_back( aMesh->AddFace (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ i3 ],
+ nodes[ i4 ]));
+ else
+ newElems.push_back( aMesh->AddFace (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ iSplit < iBestQuad ? i4 : i3 ]));
+ }
+
+ const SMDS_MeshNode* newNodes[ 4 ];
+ newNodes[ 0 ] = linkNodes[ i1 ];
+ newNodes[ 1 ] = linkNodes[ i2 ];
+ newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
+ newNodes[ 3 ] = nodes[ i4 ];
+ if (iSplit == iBestQuad)
+ newElems.push_back( aMesh->AddFace( newNodes[0], newNodes[1], newNodes[2], newNodes[3] ));
+ else
+ newElems.push_back( aMesh->AddFace( newNodes[0], newNodes[1], newNodes[2] ));
+
+ } // end if(!theFace->IsQuadratic())
+
+ else { // theFace is quadratic
+ // we have to split theFace on simple triangles and one simple quadrangle
+ int tmp = il1/2;
+ int nbshift = tmp*2;
+ // shift nodes in nodes[] by nbshift
+ int i,j;
+ for(i=0; i<nbshift; i++) {
+ const SMDS_MeshNode* n = nodes[0];
+ for(j=0; j<nbFaceNodes-1; j++) {
+ nodes[j] = nodes[j+1];
+ }
+ nodes[nbFaceNodes-1] = n;
+ }
+ il1 = il1 - nbshift;
+ // now have to insert nodes between n0 and n1 or n1 and n2 (see below)
+ // n0 n1 n2 n0 n1 n2
+ // +-----+-----+ +-----+-----+
+ // \ / | |
+ // \ / | |
+ // n5+ +n3 n7+ +n3
+ // \ / | |
+ // \ / | |
+ // + +-----+-----+
+ // n4 n6 n5 n4
+
+ // create new elements
+ int n1,n2,n3;
+ if ( nbFaceNodes == 6 ) { // quadratic triangle
+ newElems.push_back( aMesh->AddFace( nodes[3], nodes[4], nodes[5] ));
+ if ( theFace->IsMediumNode(nodes[il1]) ) {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[0], nodes[1], nodes[3], nodes[5] ));
+ n1 = 1;
+ n2 = 2;
+ n3 = 3;
+ }
+ else {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[1], nodes[2], nodes[3], nodes[5] ));
+ n1 = 0;
+ n2 = 1;
+ n3 = 5;
+ }
+ }
+ else { // nbFaceNodes==8 - quadratic quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[3], nodes[4], nodes[5] ));
+ newElems.push_back( aMesh->AddFace( nodes[5], nodes[6], nodes[7] ));
+ newElems.push_back( aMesh->AddFace( nodes[5], nodes[7], nodes[3] ));
+ if ( theFace->IsMediumNode( nodes[ il1 ])) {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[0], nodes[1], nodes[3], nodes[7] ));
+ n1 = 1;
+ n2 = 2;
+ n3 = 3;
+ }
+ else {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[1], nodes[2], nodes[3], nodes[7] ));
+ n1 = 0;
+ n2 = 1;
+ n3 = 7;
+ }
+ }
+ // create needed triangles using n1,n2,n3 and inserted nodes
+ int nbn = 2 + aNodesToInsert.size();
+ vector<const SMDS_MeshNode*> aNodes(nbn);
+ aNodes[0 ] = nodes[n1];
+ aNodes[nbn-1] = nodes[n2];
+ list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
+ for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
+ aNodes[iNode++] = *nIt;
+ }
+ for ( i = 1; i < nbn; i++ )
+ newElems.push_back( aMesh->AddFace( aNodes[i-1], aNodes[i], nodes[n3] ));
+ }
+
+ // remove the old face
+ for ( size_t i = 0; i < newElems.size(); ++i )
+ if ( newElems[i] )
+ {
+ aMesh->SetMeshElementOnShape( newElems[i], theFace->getshapeId() );
+ myLastCreatedElems.Append( newElems[i] );
+ }
+ ReplaceElemInGroups( theFace, newElems, aMesh );
+ aMesh->RemoveElement(theFace);
+
+} // InsertNodesIntoLink()
+
+//=======================================================================
+//function : UpdateVolumes
+//purpose :
+//=======================================================================
+
+void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode* theBetweenNode1,
+ const SMDS_MeshNode* theBetweenNode2,
+ list<const SMDS_MeshNode*>& theNodesToInsert)
+{
+ 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 the volume
+ SMESHDS_Mesh *aMesh = GetMeshDS();
+
+ if ( SMDS_MeshElement* newElem = aMesh->AddPolyhedralVolume( poly_nodes, quantities ))
+ {
+ aMesh->SetMeshElementOnShape( newElem, elem->getshapeId() );
+ myLastCreatedElems.Append( newElem );
+ ReplaceElemInGroups( elem, newElem, aMesh );
+ }
+ aMesh->RemoveElement( elem );
+ }
+}
+
+namespace
+{
+ //================================================================================
+ /*!
+ * \brief Transform any volume into data of SMDSEntity_Polyhedra
+ */
+ //================================================================================
+
+ void volumeToPolyhedron( const SMDS_MeshElement* elem,
+ vector<const SMDS_MeshNode *> & nodes,
+ vector<int> & nbNodeInFaces )
+ {
+ nodes.clear();
+ nbNodeInFaces.clear();
+ SMDS_VolumeTool vTool ( elem );
+ for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+ {
+ const SMDS_MeshNode** fNodes = vTool.GetFaceNodes( iF );
+ nodes.insert( nodes.end(), fNodes, fNodes + vTool.NbFaceNodes( iF ));
+ nbNodeInFaces.push_back( vTool.NbFaceNodes( iF ));
+ }
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Convert elements contained in a sub-mesh to quadratic
+ * \return int - nb of checked elements
+ */
+//=======================================================================
+
+int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm,
+ SMESH_MesherHelper& theHelper,
+ const bool theForce3d)
+{
+ int nbElem = 0;
+ if( !theSm ) return nbElem;
+
+ vector<int> nbNodeInFaces;
+ vector<const SMDS_MeshNode *> nodes;
+ SMDS_ElemIteratorPtr ElemItr = theSm->GetElements();
+ while(ElemItr->more())
+ {
+ nbElem++;
+ const SMDS_MeshElement* elem = ElemItr->next();
+ if( !elem ) continue;
+
+ // analyse a necessity of conversion
+ const SMDSAbs_ElementType aType = elem->GetType();
+ if ( aType < SMDSAbs_Edge || aType > SMDSAbs_Volume )
+ continue;
+ const SMDSAbs_EntityType aGeomType = elem->GetEntityType();
+ bool hasCentralNodes = false;
+ if ( elem->IsQuadratic() )
+ {
+ bool alreadyOK;
+ switch ( aGeomType ) {
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa:
+ alreadyOK = !theHelper.GetIsBiQuadratic(); break;
+
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ case SMDSEntity_TriQuad_Hexa:
+ alreadyOK = theHelper.GetIsBiQuadratic();
+ hasCentralNodes = true;
+ break;
+ default:
+ alreadyOK = true;
+ }
+ // take into account already present modium nodes
+ switch ( aType ) {
+ case SMDSAbs_Volume:
+ theHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( elem )); break;
+ case SMDSAbs_Face:
+ theHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( elem )); break;
+ case SMDSAbs_Edge:
+ theHelper.AddTLinks( static_cast< const SMDS_MeshEdge* >( elem )); break;
+ default:;
+ }
+ if ( alreadyOK )
+ continue;
+ }
+ // get elem data needed to re-create it
+ //
+ const int id = elem->GetID();
+ const int nbNodes = elem->NbCornerNodes();
+ nodes.assign(elem->begin_nodes(), elem->end_nodes());
+ if ( aGeomType == SMDSEntity_Polyhedra )
+ nbNodeInFaces = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
+ else if ( aGeomType == SMDSEntity_Hexagonal_Prism )
+ volumeToPolyhedron( elem, nodes, nbNodeInFaces );
+
+ // remove a linear element
+ GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
+
+ // remove central nodes of biquadratic elements (biquad->quad convertion)
+ if ( hasCentralNodes )
+ for ( size_t i = nbNodes * 2; i < nodes.size(); ++i )
+ if ( nodes[i]->NbInverseElements() == 0 )
+ GetMeshDS()->RemoveFreeNode( nodes[i], theSm, /*fromGroups=*/true );
+
+ const SMDS_MeshElement* NewElem = 0;
+
+ switch( aType )
+ {
+ case SMDSAbs_Edge :
+ {
+ NewElem = theHelper.AddEdge(nodes[0], nodes[1], id, theForce3d);
+ break;
+ }
+ case SMDSAbs_Face :
+ {
+ switch(nbNodes)
+ {
+ case 3:
+ NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
+ break;
+ case 4:
+ NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ default:
+ NewElem = theHelper.AddPolygonalFace(nodes, id, theForce3d);
+ }
+ break;
+ }
+ case SMDSAbs_Volume :
+ {
+ switch( aGeomType )
+ {
+ case SMDSEntity_Tetra:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ case SMDSEntity_Pyramid:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], id, theForce3d);
+ break;
+ case SMDSEntity_Penta:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], id, theForce3d);
+ break;
+ case SMDSEntity_Hexa:
+ case SMDSEntity_Quad_Hexa:
+ case SMDSEntity_TriQuad_Hexa:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+ break;
+ case SMDSEntity_Hexagonal_Prism:
+ default:
+ NewElem = theHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
+ }
+ break;
+ }
+ default :
+ continue;
+ }
+ ReplaceElemInGroups( elem, NewElem, GetMeshDS());
+ if( NewElem && NewElem->getshapeId() < 1 )
+ theSm->AddElement( NewElem );
+ }
+ return nbElem;
+}
+//=======================================================================
+//function : ConvertToQuadratic
+//purpose :
+//=======================================================================
+
+void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theToBiQuad)
+{
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+
+ SMESH_MesherHelper aHelper(*myMesh);
+
+ aHelper.SetIsQuadratic( true );
+ aHelper.SetIsBiQuadratic( theToBiQuad );
+ aHelper.SetElementsOnShape(true);
+ aHelper.ToFixNodeParameters( true );
+
+ // convert elements assigned to sub-meshes
+ int nbCheckedElems = 0;
+ if ( myMesh->HasShapeToMesh() )
+ {
+ if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
+ {
+ SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
+ while ( smIt->more() ) {
+ SMESH_subMesh* sm = smIt->next();
+ if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() ) {
+ aHelper.SetSubShape( sm->GetSubShape() );
+ nbCheckedElems += convertElemToQuadratic(smDS, aHelper, theForce3d);
+ }
+ }
+ }
+ }
+
+ // convert elements NOT assigned to sub-meshes
+ int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
+ if ( nbCheckedElems < totalNbElems ) // not all elements are in sub-meshes
+ {
+ aHelper.SetElementsOnShape(false);
+ SMESHDS_SubMesh *smDS = 0;
+
+ // convert edges
+ SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator();
+ while( aEdgeItr->more() )
+ {
+ const SMDS_MeshEdge* edge = aEdgeItr->next();
+ if ( !edge->IsQuadratic() )
+ {
+ int id = edge->GetID();
+ const SMDS_MeshNode* n1 = edge->GetNode(0);
+ const SMDS_MeshNode* n2 = edge->GetNode(1);
+
+ meshDS->RemoveFreeElement(edge, smDS, /*fromGroups=*/false);
+
+ const SMDS_MeshEdge* NewEdge = aHelper.AddEdge(n1, n2, id, theForce3d);
+ ReplaceElemInGroups( edge, NewEdge, GetMeshDS());
+ }
+ else
+ {
+ aHelper.AddTLinks( static_cast< const SMDS_MeshEdge* >( edge ));
+ }
+ }
+
+ // convert faces
+ SMDS_FaceIteratorPtr aFaceItr = meshDS->facesIterator();
+ while( aFaceItr->more() )
+ {
+ const SMDS_MeshFace* face = aFaceItr->next();
+ if ( !face ) continue;
+
+ const SMDSAbs_EntityType type = face->GetEntityType();
+ bool alreadyOK;
+ switch( type )
+ {
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ alreadyOK = !theToBiQuad;
+ aHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( face ));
+ break;
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ alreadyOK = theToBiQuad;
+ aHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( face ));
+ break;
+ default: alreadyOK = false;
+ }
+ if ( alreadyOK )
+ continue;
+
+ const int id = face->GetID();
+ vector<const SMDS_MeshNode *> nodes ( face->begin_nodes(), face->end_nodes());
+
+ meshDS->RemoveFreeElement(face, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshFace * NewFace = 0;
+ switch( type )
+ {
+ case SMDSEntity_Triangle:
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_BiQuad_Triangle:
+ NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
+ if ( nodes.size() == 7 && nodes[6]->NbInverseElements() == 0 ) // rm a central node
+ GetMeshDS()->RemoveFreeNode( nodes[6], /*sm=*/0, /*fromGroups=*/true );
+ break;
+
+ case SMDSEntity_Quadrangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ if ( nodes.size() == 9 && nodes[8]->NbInverseElements() == 0 ) // rm a central node
+ GetMeshDS()->RemoveFreeNode( nodes[8], /*sm=*/0, /*fromGroups=*/true );
+ break;
+
+ default:;
+ NewFace = aHelper.AddPolygonalFace(nodes, id, theForce3d);
+ }
+ ReplaceElemInGroups( face, NewFace, GetMeshDS());
+ }
+
+ // convert volumes
+ vector<int> nbNodeInFaces;
+ SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator();
+ while(aVolumeItr->more())
+ {
+ const SMDS_MeshVolume* volume = aVolumeItr->next();
+ if ( !volume ) continue;
+
+ const SMDSAbs_EntityType type = volume->GetEntityType();
+ if ( volume->IsQuadratic() )
+ {
+ bool alreadyOK;
+ switch ( type )
+ {
+ case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break;
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break;
+ default: alreadyOK = true;
+ }
+ if ( alreadyOK )
+ {
+ aHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( volume ));
+ continue;
+ }
+ }
+ const int id = volume->GetID();
+ vector<const SMDS_MeshNode *> nodes (volume->begin_nodes(), volume->end_nodes());
+ if ( type == SMDSEntity_Polyhedra )
+ nbNodeInFaces = static_cast<const SMDS_VtkVolume* >(volume)->GetQuantities();
+ else if ( type == SMDSEntity_Hexagonal_Prism )
+ volumeToPolyhedron( volume, nodes, nbNodeInFaces );
+
+ meshDS->RemoveFreeElement(volume, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshVolume * NewVolume = 0;
+ switch ( type )
+ {
+ case SMDSEntity_Tetra:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d );
+ break;
+ case SMDSEntity_Hexa:
+ case SMDSEntity_Quad_Hexa:
+ case SMDSEntity_TriQuad_Hexa:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+ for ( size_t i = 20; i < nodes.size(); ++i ) // rm central nodes
+ if ( nodes[i]->NbInverseElements() == 0 )
+ GetMeshDS()->RemoveFreeNode( nodes[i], /*sm=*/0, /*fromGroups=*/true );
+ break;
+ case SMDSEntity_Pyramid:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+ nodes[3], nodes[4], id, theForce3d);
+ break;
+ case SMDSEntity_Penta:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+ nodes[3], nodes[4], nodes[5], id, theForce3d);
+ break;
+ case SMDSEntity_Hexagonal_Prism:
+ default:
+ NewVolume = aHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
+ }
+ ReplaceElemInGroups(volume, NewVolume, meshDS);
+ }
+ }
+
+ if ( !theForce3d )
+ { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
+ // aHelper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
+ // aHelper.FixQuadraticElements(myError);
+ SMESH_MesherHelper( *myMesh ).FixQuadraticElements(myError);
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Makes given elements quadratic
+ * \param theForce3d - if true, the medium nodes will be placed in the middle of link
+ * \param theElements - elements to make quadratic
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d,
+ TIDSortedElemSet& theElements,
+ const bool theToBiQuad)
+{
+ if ( theElements.empty() ) return;
+
+ // we believe that all theElements are of the same type
+ const SMDSAbs_ElementType elemType = (*theElements.begin())->GetType();
+
+ // get all nodes shared by theElements
+ TIDSortedNodeSet allNodes;
+ TIDSortedElemSet::iterator eIt = theElements.begin();
+ for ( ; eIt != theElements.end(); ++eIt )
+ allNodes.insert( (*eIt)->begin_nodes(), (*eIt)->end_nodes() );
+
+ // complete theElements with elements of lower dim whose all nodes are in allNodes
+
+ TIDSortedElemSet quadAdjacentElems [ SMDSAbs_NbElementTypes ]; // quadratic adjacent elements
+ TIDSortedElemSet checkedAdjacentElems [ SMDSAbs_NbElementTypes ];
+ TIDSortedNodeSet::iterator nIt = allNodes.begin();
+ for ( ; nIt != allNodes.end(); ++nIt )
+ {
+ const SMDS_MeshNode* n = *nIt;
+ SMDS_ElemIteratorPtr invIt = n->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* e = invIt->next();
+ const SMDSAbs_ElementType type = e->GetType();
+ if ( e->IsQuadratic() )
+ {
+ quadAdjacentElems[ type ].insert( e );
+
+ bool alreadyOK;
+ switch ( e->GetEntityType() ) {
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break;
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break;
+ default: alreadyOK = true;
+ }
+ if ( alreadyOK )
+ continue;
+ }
+ if ( type >= elemType )
+ continue; // same type or more complex linear element
+
+ if ( !checkedAdjacentElems[ type ].insert( e ).second )
+ continue; // e is already checked
+
+ // check nodes
+ bool allIn = true;
+ SMDS_NodeIteratorPtr nodeIt = e->nodeIterator();
+ while ( nodeIt->more() && allIn )
+ allIn = allNodes.count( nodeIt->next() );
+ if ( allIn )
+ theElements.insert(e );
+ }
+ }
+
+ SMESH_MesherHelper helper(*myMesh);
+ helper.SetIsQuadratic( true );
+ helper.SetIsBiQuadratic( theToBiQuad );
+
+ // add links of quadratic adjacent elements to the helper
+
+ if ( !quadAdjacentElems[SMDSAbs_Edge].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Edge].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Edge].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshEdge*> (*eIt) );
+ }
+ if ( !quadAdjacentElems[SMDSAbs_Face].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Face].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Face].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshFace*> (*eIt) );
+ }
+ if ( !quadAdjacentElems[SMDSAbs_Volume].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Volume].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Volume].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshVolume*> (*eIt) );
+ }
+
+ // make quadratic (or bi-tri-quadratic) elements instead of linear ones
+
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+ SMESHDS_SubMesh* smDS = 0;
+ for ( eIt = theElements.begin(); eIt != theElements.end(); ++eIt )
+ {
+ const SMDS_MeshElement* elem = *eIt;
+
+ bool alreadyOK;
+ int nbCentralNodes = 0;
+ switch ( elem->GetEntityType() ) {
+ // linear convertible
+ case SMDSEntity_Edge:
+ case SMDSEntity_Triangle:
+ case SMDSEntity_Quadrangle:
+ case SMDSEntity_Tetra:
+ case SMDSEntity_Pyramid:
+ case SMDSEntity_Hexa:
+ case SMDSEntity_Penta: alreadyOK = false; nbCentralNodes = 0; break;
+ // quadratic that can become bi-quadratic
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa: alreadyOK =!theToBiQuad; nbCentralNodes = 0; break;
+ // bi-quadratic
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle: alreadyOK = theToBiQuad; nbCentralNodes = 1; break;
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; nbCentralNodes = 7; break;
+ // the rest
+ default: alreadyOK = true;
+ }
+ if ( alreadyOK ) continue;
+
+ const SMDSAbs_ElementType type = elem->GetType();
+ const int id = elem->GetID();
+ const int nbNodes = elem->NbCornerNodes();
+ vector<const SMDS_MeshNode *> nodes ( elem->begin_nodes(), elem->end_nodes());
+
+ helper.SetSubShape( elem->getshapeId() );
+
+ if ( !smDS || !smDS->Contains( elem ))
+ smDS = meshDS->MeshElements( elem->getshapeId() );
+ meshDS->RemoveFreeElement(elem, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshElement * newElem = 0;
+ switch( nbNodes )
+ {
+ case 4: // cases for most frequently used element types go first (for optimization)
+ if ( type == SMDSAbs_Volume )
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ else
+ newElem = helper.AddFace (nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ case 8:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+ break;
+ case 3:
+ newElem = helper.AddFace (nodes[0], nodes[1], nodes[2], id, theForce3d);
+ break;
+ case 2:
+ newElem = helper.AddEdge(nodes[0], nodes[1], id, theForce3d);
+ break;
+ case 5:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], id, theForce3d);
+ break;
+ case 6:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], id, theForce3d);
+ break;
+ default:;
+ }
+ ReplaceElemInGroups( elem, newElem, meshDS);
+ if( newElem && smDS )
+ smDS->AddElement( newElem );
+
+ // remove central nodes
+ for ( size_t i = nodes.size() - nbCentralNodes; i < nodes.size(); ++i )
+ if ( nodes[i]->NbInverseElements() == 0 )
+ meshDS->RemoveFreeNode( nodes[i], smDS, /*fromGroups=*/true );
+
+ } // loop on theElements
+
+ if ( !theForce3d )
+ { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
+ // helper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
+ // helper.FixQuadraticElements( myError );
+ SMESH_MesherHelper( *myMesh ).FixQuadraticElements(myError);
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Convert quadratic elements to linear ones and remove quadratic nodes
+ * \return int - nb of checked elements
+ */
+//=======================================================================
+
+int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh * theSm,
+ SMDS_ElemIteratorPtr theItr,
+ const int theShapeID)
+{
+ int nbElem = 0;
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+ ElemFeatures elemType;
+ vector<const SMDS_MeshNode *> nodes;
+
+ while( theItr->more() )
+ {
+ const SMDS_MeshElement* elem = theItr->next();
+ nbElem++;
+ if( elem && elem->IsQuadratic())
+ {
+ // get elem data
+ int nbCornerNodes = elem->NbCornerNodes();
+ nodes.assign( elem->begin_nodes(), elem->end_nodes() );
+
+ elemType.Init( elem, /*basicOnly=*/false ).SetID( elem->GetID() ).SetQuad( false );
+
+ //remove a quadratic element
+ if ( !theSm || !theSm->Contains( elem ))
+ theSm = meshDS->MeshElements( elem->getshapeId() );
+ meshDS->RemoveFreeElement( elem, theSm, /*fromGroups=*/false );
+
+ // remove medium nodes
+ for ( size_t i = nbCornerNodes; i < nodes.size(); ++i )
+ if ( nodes[i]->NbInverseElements() == 0 )
+ meshDS->RemoveFreeNode( nodes[i], theSm );
+
+ // add a linear element
+ nodes.resize( nbCornerNodes );
+ SMDS_MeshElement * newElem = AddElement( nodes, elemType );
+ ReplaceElemInGroups(elem, newElem, meshDS);
+ if( theSm && newElem )
+ theSm->AddElement( newElem );
+ }
+ }
+ return nbElem;
+}
+
+//=======================================================================
+//function : ConvertFromQuadratic
+//purpose :
+//=======================================================================
+
+bool SMESH_MeshEditor::ConvertFromQuadratic()
+{
+ int nbCheckedElems = 0;
+ if ( myMesh->HasShapeToMesh() )
+ {
+ if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
+ {
+ SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
+ while ( smIt->more() ) {
+ SMESH_subMesh* sm = smIt->next();
+ if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() )
+ nbCheckedElems += removeQuadElem( smDS, smDS->GetElements(), sm->GetId() );
+ }
+ }
+ }
+
+ int totalNbElems =
+ GetMeshDS()->NbEdges() + GetMeshDS()->NbFaces() + GetMeshDS()->NbVolumes();
+ if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
+ {
+ SMESHDS_SubMesh *aSM = 0;
+ removeQuadElem( aSM, GetMeshDS()->elementsIterator(), 0 );
+ }
+
+ return true;
+}
+
+namespace
+{
+ //================================================================================
+ /*!
+ * \brief Return true if all medium nodes of the element are in the node set
+ */
+ //================================================================================
+
+ bool allMediumNodesIn(const SMDS_MeshElement* elem, TIDSortedNodeSet& nodeSet )
+ {
+ for ( int i = elem->NbCornerNodes(); i < elem->NbNodes(); ++i )
+ if ( !nodeSet.count( elem->GetNode(i) ))
+ return false;
+ return true;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Makes given elements linear
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
+{
+ if ( theElements.empty() ) return;
+
+ // collect IDs of medium nodes of theElements; some of these nodes will be removed
+ set<int> mediumNodeIDs;
+ TIDSortedElemSet::iterator eIt = theElements.begin();
+ for ( ; eIt != theElements.end(); ++eIt )
+ {
+ const SMDS_MeshElement* e = *eIt;
+ for ( int i = e->NbCornerNodes(); i < e->NbNodes(); ++i )
+ mediumNodeIDs.insert( e->GetNode(i)->GetID() );
+ }
+
+ // replace given elements by linear ones
+ SMDS_ElemIteratorPtr elemIt = elemSetIterator( theElements );
+ removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
+
+ // we need to convert remaining elements whose all medium nodes are in mediumNodeIDs
+ // except those elements sharing medium nodes of quadratic element whose medium nodes
+ // are not all in mediumNodeIDs
+
+ // get remaining medium nodes
+ TIDSortedNodeSet mediumNodes;
+ set<int>::iterator nIdsIt = mediumNodeIDs.begin();
+ for ( ; nIdsIt != mediumNodeIDs.end(); ++nIdsIt )
+ if ( const SMDS_MeshNode* n = GetMeshDS()->FindNode( *nIdsIt ))
+ mediumNodes.insert( mediumNodes.end(), n );
+
+ // find more quadratic elements to convert
+ TIDSortedElemSet moreElemsToConvert;
+ TIDSortedNodeSet::iterator nIt = mediumNodes.begin();
+ for ( ; nIt != mediumNodes.end(); ++nIt )
+ {
+ SMDS_ElemIteratorPtr invIt = (*nIt)->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* e = invIt->next();
+ if ( e->IsQuadratic() && allMediumNodesIn( e, mediumNodes ))
+ {
+ // find a more complex element including e and
+ // whose medium nodes are not in mediumNodes
+ bool complexFound = false;
+ for ( int type = e->GetType() + 1; type < SMDSAbs_0DElement; ++type )
+ {
+ SMDS_ElemIteratorPtr invIt2 =
+ (*nIt)->GetInverseElementIterator( SMDSAbs_ElementType( type ));
+ while ( invIt2->more() )
+ {
+ const SMDS_MeshElement* eComplex = invIt2->next();
+ if ( eComplex->IsQuadratic() && !allMediumNodesIn( eComplex, mediumNodes))
+ {
+ int nbCommonNodes = SMESH_MeshAlgos::GetCommonNodes( e, eComplex ).size();
+ if ( nbCommonNodes == e->NbNodes())
+ {
+ complexFound = true;
+ type = SMDSAbs_NbElementTypes; // to quit from the outer loop
+ break;
+ }
+ }
+ }
+ }
+ if ( !complexFound )
+ moreElemsToConvert.insert( e );
+ }
+ }
+ }
+ elemIt = elemSetIterator( moreElemsToConvert );
+ removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
+}
+
+//=======================================================================
+//function : SewSideElements
+//purpose :
+//=======================================================================
+
+SMESH_MeshEditor::Sew_Error
+SMESH_MeshEditor::SewSideElements (TIDSortedElemSet& theSide1,
+ TIDSortedElemSet& theSide2,
+ const SMDS_MeshNode* theFirstNode1,
+ const SMDS_MeshNode* theFirstNode2,
+ const SMDS_MeshNode* theSecondNode1,
+ const SMDS_MeshNode* theSecondNode2)
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ MESSAGE ("::::SewSideElements()");
+ if ( theSide1.size() != theSide2.size() )
+ return SEW_DIFF_NB_OF_ELEMENTS;
+
+ Sew_Error aResult = SEW_OK;
+ // Algo:
+ // 1. Build set of faces representing each side
+ // 2. Find which nodes of the side 1 to merge with ones on the side 2
+ // 3. Replace nodes in elements of the side 1 and remove replaced nodes
+
+ // =======================================================================
+ // 1. Build set of faces representing each side:
+ // =======================================================================
+ // a. build set of nodes belonging to faces
+ // b. complete set of faces: find missing faces whose nodes are in set of nodes
+ // c. create temporary faces representing side of volumes if correspondent
+ // face does not exist
+
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+ // TODO algoritm not OK with vtkUnstructuredGrid: 2 meshes can't share nodes
+ //SMDS_Mesh aTmpFacesMesh; // try to use the same mesh
+ TIDSortedElemSet faceSet1, faceSet2;
+ set<const SMDS_MeshElement*> volSet1, volSet2;
+ set<const SMDS_MeshNode*> nodeSet1, nodeSet2;
+ TIDSortedElemSet * faceSetPtr[] = { &faceSet1, &faceSet2 };
+ set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
+ set<const SMDS_MeshNode*> * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
+ TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
+ int iSide, iFace, iNode;
+
+ list<const SMDS_MeshElement* > tempFaceList;
+ for ( iSide = 0; iSide < 2; iSide++ ) {
+ set<const SMDS_MeshNode*> * nodeSet = nodeSetPtr[ iSide ];
+ TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
+ TIDSortedElemSet * faceSet = faceSetPtr[ iSide ];
+ set<const SMDS_MeshElement*> * volSet = volSetPtr [ iSide ];
+ set<const SMDS_MeshElement*>::iterator vIt;
+ TIDSortedElemSet::iterator eIt;
+ set<const SMDS_MeshNode*>::iterator nIt;
+
+ // check that given nodes belong to given elements
+ const SMDS_MeshNode* n1 = ( iSide == 0 ) ? theFirstNode1 : theFirstNode2;
+ const SMDS_MeshNode* n2 = ( iSide == 0 ) ? theSecondNode1 : theSecondNode2;
+ int firstIndex = -1, secondIndex = -1;
+ for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
+ const SMDS_MeshElement* elem = *eIt;
+ if ( firstIndex < 0 ) firstIndex = elem->GetNodeIndex( n1 );
+ if ( secondIndex < 0 ) secondIndex = elem->GetNodeIndex( n2 );
+ if ( firstIndex > -1 && secondIndex > -1 ) break;
+ }
+ if ( firstIndex < 0 || secondIndex < 0 ) {
+ // we can simply return until temporary faces created
+ return (iSide == 0 ) ? SEW_BAD_SIDE1_NODES : SEW_BAD_SIDE2_NODES;
+ }
+
+ // -----------------------------------------------------------
+ // 1a. Collect nodes of existing faces
+ // and build set of face nodes in order to detect missing
+ // faces corresponding to sides of volumes
+ // -----------------------------------------------------------
+
+ set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
+
+ // loop on the given element of a side
+ for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
+ //const SMDS_MeshElement* elem = *eIt;
+ const SMDS_MeshElement* elem = *eIt;
+ if ( elem->GetType() == SMDSAbs_Face ) {
+ faceSet->insert( elem );
+ set <const SMDS_MeshNode*> faceNodeSet;
+ SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+ while ( nodeIt->more() ) {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ nodeSet->insert( n );
+ faceNodeSet.insert( n );
+ }
+ setOfFaceNodeSet.insert( faceNodeSet );
+ }
+ else if ( elem->GetType() == SMDSAbs_Volume )
+ volSet->insert( elem );
+ }
+ // ------------------------------------------------------------------------------
+ // 1b. Complete set of faces: find missing faces whose nodes are in set of nodes
+ // ------------------------------------------------------------------------------
+
+ for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
+ SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() ) { // loop on faces sharing a node
+ const SMDS_MeshElement* f = fIt->next();
+ if ( faceSet->find( f ) == faceSet->end() ) {
+ // check if all nodes are in nodeSet and
+ // complete setOfFaceNodeSet if they are
+ set <const SMDS_MeshNode*> faceNodeSet;
+ SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
+ bool allInSet = true;
+ while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ if ( nodeSet->find( n ) == nodeSet->end() )
+ allInSet = false;
+ else
+ faceNodeSet.insert( n );
+ }
+ if ( allInSet ) {
+ faceSet->insert( f );
+ setOfFaceNodeSet.insert( faceNodeSet );