+ // 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 ( size_t i = 0; i < myLastCreatedElems.size(); ++i )
+ {
+ if ( !myLastCreatedElems[i] ) continue;
+ if ( myLastCreatedElems[i]->GetType() == SMDSAbs_Edge )
+ segments.insert( segments.end(), myLastCreatedElems[i] );
+ else
+ newFaces.push_back( myLastCreatedElems[i] );
+ }
+ // get segments adjacent to merged nodes
+ TListOfListOfNodes::iterator groupIt = nodeGroupsToMerge.begin();
+ for ( ; groupIt != nodeGroupsToMerge.end(); groupIt++ )
+ {
+ const list<const SMDS_MeshNode*>& nodes = *groupIt;
+ if ( nodes.front()->IsNull() ) continue;
+ SMDS_ElemIteratorPtr segIt = nodes.front()->GetInverseElementIterator( SMDSAbs_Edge );
+ while ( segIt->more() )
+ segments.insert( segIt->next() );
+ }
+
+ // find coincident
+ TListOfListOfElementsID equalGroups;
+ if ( !segments.empty() )
+ FindEqualElements( segments, equalGroups );
+ if ( !equalGroups.empty() )
+ {
+ // remove from segments those that will be removed
+ TListOfListOfElementsID::iterator itGroups = equalGroups.begin();
+ for ( ; itGroups != equalGroups.end(); ++itGroups )
+ {
+ list< int >& group = *itGroups;
+ list< int >::iterator id = group.begin();
+ for ( ++id; id != group.end(); ++id )
+ if ( const SMDS_MeshElement* seg = GetMeshDS()->FindElement( *id ))
+ segments.erase( seg );
+ }
+ // remove equal segments
+ MergeElements( equalGroups );
+
+ // restore myLastCreatedElems
+ myLastCreatedElems = newFaces;
+ TIDSortedElemSet::iterator seg = segments.begin();
+ for ( ; seg != segments.end(); ++seg )
+ myLastCreatedElems.push_back( *seg );
+ }
+
+ return aResult;
+}
+
+//=======================================================================
+//function : InsertNodesIntoLink
+//purpose : insert theNodesToInsert into theElement between theBetweenNode1
+// and theBetweenNode2 and split theElement
+//=======================================================================
+
+void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theElement,
+ const SMDS_MeshNode* theBetweenNode1,
+ const SMDS_MeshNode* theBetweenNode2,
+ list<const SMDS_MeshNode*>& theNodesToInsert,
+ const bool toCreatePoly)
+{
+ if ( !theElement ) return;
+
+ SMESHDS_Mesh *aMesh = GetMeshDS();
+ vector<const SMDS_MeshElement*> newElems;
+
+ if ( theElement->GetType() == SMDSAbs_Edge )
+ {
+ theNodesToInsert.push_front( theBetweenNode1 );
+ theNodesToInsert.push_back ( theBetweenNode2 );
+ list<const SMDS_MeshNode*>::iterator n = theNodesToInsert.begin();
+ const SMDS_MeshNode* n1 = *n;
+ for ( ++n; n != theNodesToInsert.end(); ++n )
+ {
+ const SMDS_MeshNode* n2 = *n;
+ if ( const SMDS_MeshElement* seg = aMesh->FindEdge( n1, n2 ))
+ AddToSameGroups( seg, theElement, aMesh );
+ else
+ newElems.push_back( aMesh->AddEdge ( n1, n2 ));
+ n1 = n2;
+ }
+ theNodesToInsert.pop_front();
+ theNodesToInsert.pop_back();
+
+ if ( theElement->IsQuadratic() ) // add a not split part
+ {
+ vector<const SMDS_MeshNode*> nodes( theElement->begin_nodes(),
+ theElement->end_nodes() );
+ int iOther = 0, nbN = nodes.size();
+ for ( ; iOther < nbN; ++iOther )
+ if ( nodes[iOther] != theBetweenNode1 &&
+ nodes[iOther] != theBetweenNode2 )
+ break;
+ if ( iOther == 0 )
+ {
+ if ( const SMDS_MeshElement* seg = aMesh->FindEdge( nodes[0], nodes[1] ))
+ AddToSameGroups( seg, theElement, aMesh );
+ else
+ newElems.push_back( aMesh->AddEdge ( nodes[0], nodes[1] ));
+ }
+ else if ( iOther == 2 )
+ {
+ if ( const SMDS_MeshElement* seg = aMesh->FindEdge( nodes[1], nodes[2] ))
+ AddToSameGroups( seg, theElement, aMesh );
+ else
+ newElems.push_back( aMesh->AddEdge ( nodes[1], nodes[2] ));
+ }
+ }
+ // treat new elements
+ for ( size_t i = 0; i < newElems.size(); ++i )
+ if ( newElems[i] )
+ {
+ aMesh->SetMeshElementOnShape( newElems[i], theElement->getshapeId() );
+ myLastCreatedElems.push_back( newElems[i] );
+ }
+ ReplaceElemInGroups( theElement, newElems, aMesh );
+ aMesh->RemoveElement( theElement );
+ return;
+
+ } // if ( theElement->GetType() == SMDSAbs_Edge )
+
+ const SMDS_MeshElement* theFace = theElement;
+ if ( theFace->GetType() != SMDSAbs_Face ) return;
+
+ // find indices of 2 link nodes and of the rest nodes
+ int iNode = 0, il1, il2, i3, i4;
+ il1 = il2 = i3 = i4 = -1;
+ vector<const SMDS_MeshNode*> nodes( theFace->NbNodes() );
+
+ SMDS_NodeIteratorPtr nodeIt = theFace->interlacedNodesIterator();
+ while ( nodeIt->more() ) {
+ const SMDS_MeshNode* n = nodeIt->next();
+ if ( n == theBetweenNode1 )
+ il1 = iNode;
+ else if ( n == theBetweenNode2 )
+ il2 = iNode;
+ else if ( i3 < 0 )
+ i3 = iNode;
+ else
+ i4 = iNode;
+ nodes[ iNode++ ] = n;
+ }
+ if ( il1 < 0 || il2 < 0 || i3 < 0 )
+ return ;
+
+ // arrange link nodes to go one after another regarding the face orientation
+ bool reverse = ( Abs( il2 - il1 ) == 1 ? il2 < il1 : il1 < il2 );
+ list<const SMDS_MeshNode *> aNodesToInsert = theNodesToInsert;
+ if ( reverse ) {
+ iNode = il1;
+ il1 = il2;
+ il2 = iNode;
+ aNodesToInsert.reverse();
+ }
+ // check that not link nodes of a quadrangles are in good order
+ int nbFaceNodes = theFace->NbNodes();
+ if ( nbFaceNodes == 4 && i4 - i3 != 1 ) {
+ iNode = i3;
+ i3 = i4;
+ i4 = iNode;
+ }
+
+ if (toCreatePoly || theFace->IsPoly()) {
+
+ iNode = 0;
+ vector<const SMDS_MeshNode *> poly_nodes (nbFaceNodes + aNodesToInsert.size());
+
+ // add nodes of face up to first node of link
+ bool isFLN = false;
+ SMDS_NodeIteratorPtr nodeIt = theFace->interlacedNodesIterator();
+ while ( nodeIt->more() && !isFLN ) {
+ const SMDS_MeshNode* n = nodeIt->next();
+ poly_nodes[iNode++] = n;
+ isFLN = ( n == nodes[il1] );
+ }
+ // add nodes to insert
+ list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
+ for (; nIt != aNodesToInsert.end(); nIt++) {
+ poly_nodes[iNode++] = *nIt;
+ }
+ // add nodes of face starting from last node of link
+ while ( nodeIt->more() ) {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ poly_nodes[iNode++] = n;
+ }
+
+ // make a new face
+ newElems.push_back( aMesh->AddPolygonalFace( poly_nodes ));
+ }
+
+ else if ( !theFace->IsQuadratic() )
+ {
+ // put aNodesToInsert between theBetweenNode1 and theBetweenNode2
+ int nbLinkNodes = 2 + aNodesToInsert.size();
+ //const SMDS_MeshNode* linkNodes[ nbLinkNodes ];
+ vector<const SMDS_MeshNode*> linkNodes( nbLinkNodes );
+ linkNodes[ 0 ] = nodes[ il1 ];
+ linkNodes[ nbLinkNodes - 1 ] = nodes[ il2 ];
+ list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
+ for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
+ linkNodes[ iNode++ ] = *nIt;
+ }
+ // decide how to split a quadrangle: compare possible variants
+ // and choose which of splits to be a quadrangle
+ int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad = 0;
+ if ( nbFaceNodes == 3 ) {
+ iBestQuad = nbSplits;
+ i4 = i3;
+ }
+ else if ( nbFaceNodes == 4 ) {
+ SMESH::Controls::NumericalFunctorPtr aCrit( new SMESH::Controls::AspectRatio);
+ double aBestRate = DBL_MAX;
+ for ( int iQuad = 0; iQuad < nbSplits; iQuad++ ) {
+ i1 = 0; i2 = 1;
+ double aBadRate = 0;
+ // evaluate elements quality
+ for ( iSplit = 0; iSplit < nbSplits; iSplit++ ) {
+ if ( iSplit == iQuad ) {
+ SMDS_FaceOfNodes quad (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ i3 ],
+ nodes[ i4 ]);
+ aBadRate += getBadRate( &quad, aCrit );
+ }
+ else {
+ SMDS_FaceOfNodes tria (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ iSplit < iQuad ? i4 : i3 ]);
+ aBadRate += getBadRate( &tria, aCrit );
+ }
+ }
+ // choice
+ if ( aBadRate < aBestRate ) {
+ iBestQuad = iQuad;
+ aBestRate = aBadRate;
+ }
+ }
+ }
+
+ // create new elements
+ i1 = 0; i2 = 1;
+ for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ )
+ {
+ if ( iSplit == iBestQuad )
+ newElems.push_back( aMesh->AddFace (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ i3 ],
+ nodes[ i4 ]));
+ else
+ newElems.push_back( aMesh->AddFace (linkNodes[ i1++ ],
+ linkNodes[ i2++ ],
+ nodes[ iSplit < iBestQuad ? i4 : i3 ]));
+ }
+
+ const SMDS_MeshNode* newNodes[ 4 ];
+ newNodes[ 0 ] = linkNodes[ i1 ];
+ newNodes[ 1 ] = linkNodes[ i2 ];
+ newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
+ newNodes[ 3 ] = nodes[ i4 ];
+ if (iSplit == iBestQuad)
+ newElems.push_back( aMesh->AddFace( newNodes[0], newNodes[1], newNodes[2], newNodes[3] ));
+ else
+ newElems.push_back( aMesh->AddFace( newNodes[0], newNodes[1], newNodes[2] ));
+
+ } // end if(!theFace->IsQuadratic())
+
+ else { // theFace is quadratic
+ // we have to split theFace on simple triangles and one simple quadrangle
+ int tmp = il1/2;
+ int nbshift = tmp*2;
+ // shift nodes in nodes[] by nbshift
+ int i,j;
+ for(i=0; i<nbshift; i++) {
+ const SMDS_MeshNode* n = nodes[0];
+ for(j=0; j<nbFaceNodes-1; j++) {
+ nodes[j] = nodes[j+1];
+ }
+ nodes[nbFaceNodes-1] = n;
+ }
+ il1 = il1 - nbshift;
+ // now have to insert nodes between n0 and n1 or n1 and n2 (see below)
+ // n0 n1 n2 n0 n1 n2
+ // +-----+-----+ +-----+-----+
+ // \ / | |
+ // \ / | |
+ // n5+ +n3 n7+ +n3
+ // \ / | |
+ // \ / | |
+ // + +-----+-----+
+ // n4 n6 n5 n4
+
+ // create new elements
+ int n1,n2,n3;
+ if ( nbFaceNodes == 6 ) { // quadratic triangle
+ newElems.push_back( aMesh->AddFace( nodes[3], nodes[4], nodes[5] ));
+ if ( theFace->IsMediumNode(nodes[il1]) ) {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[0], nodes[1], nodes[3], nodes[5] ));
+ n1 = 1;
+ n2 = 2;
+ n3 = 3;
+ }
+ else {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[1], nodes[2], nodes[3], nodes[5] ));
+ n1 = 0;
+ n2 = 1;
+ n3 = 5;
+ }
+ }
+ else { // nbFaceNodes==8 - quadratic quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[3], nodes[4], nodes[5] ));
+ newElems.push_back( aMesh->AddFace( nodes[5], nodes[6], nodes[7] ));
+ newElems.push_back( aMesh->AddFace( nodes[5], nodes[7], nodes[3] ));
+ if ( theFace->IsMediumNode( nodes[ il1 ])) {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[0], nodes[1], nodes[3], nodes[7] ));
+ n1 = 1;
+ n2 = 2;
+ n3 = 3;
+ }
+ else {
+ // create quadrangle
+ newElems.push_back( aMesh->AddFace( nodes[1], nodes[2], nodes[3], nodes[7] ));
+ n1 = 0;
+ n2 = 1;
+ n3 = 7;
+ }
+ }
+ // create needed triangles using n1,n2,n3 and inserted nodes
+ int nbn = 2 + aNodesToInsert.size();
+ vector<const SMDS_MeshNode*> aNodes(nbn);
+ aNodes[0 ] = nodes[n1];
+ aNodes[nbn-1] = nodes[n2];
+ list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
+ for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
+ aNodes[iNode++] = *nIt;
+ }
+ for ( i = 1; i < nbn; i++ )
+ newElems.push_back( aMesh->AddFace( aNodes[i-1], aNodes[i], nodes[n3] ));
+ }
+
+ // remove the old face
+ for ( size_t i = 0; i < newElems.size(); ++i )
+ if ( newElems[i] )
+ {
+ aMesh->SetMeshElementOnShape( newElems[i], theFace->getshapeId() );
+ myLastCreatedElems.push_back( newElems[i] );
+ }
+ ReplaceElemInGroups( theFace, newElems, aMesh );
+ aMesh->RemoveElement(theFace);
+
+} // InsertNodesIntoLink()
+
+//=======================================================================
+//function : UpdateVolumes
+//purpose :
+//=======================================================================
+
+void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode* theBetweenNode1,
+ const SMDS_MeshNode* theBetweenNode2,
+ list<const SMDS_MeshNode*>& theNodesToInsert)
+{
+ ClearLastCreated();
+
+ SMDS_ElemIteratorPtr invElemIt = theBetweenNode1->GetInverseElementIterator(SMDSAbs_Volume);
+ while (invElemIt->more()) { // loop on inverse elements of theBetweenNode1
+ const SMDS_MeshElement* elem = invElemIt->next();
+
+ // check, if current volume has link theBetweenNode1 - theBetweenNode2
+ SMDS_VolumeTool aVolume (elem);
+ if (!aVolume.IsLinked(theBetweenNode1, theBetweenNode2))
+ continue;
+
+ // insert new nodes in all faces of the volume, sharing link theBetweenNode1 - theBetweenNode2
+ int iface, nbFaces = aVolume.NbFaces();
+ vector<const SMDS_MeshNode *> poly_nodes;
+ vector<int> quantities (nbFaces);
+
+ for (iface = 0; iface < nbFaces; iface++) {
+ int nbFaceNodes = aVolume.NbFaceNodes(iface), nbInserted = 0;
+ // faceNodes will contain (nbFaceNodes + 1) nodes, last = first
+ const SMDS_MeshNode** faceNodes = aVolume.GetFaceNodes(iface);
+
+ for (int inode = 0; inode < nbFaceNodes; inode++) {
+ poly_nodes.push_back(faceNodes[inode]);
+
+ if (nbInserted == 0) {
+ if (faceNodes[inode] == theBetweenNode1) {
+ if (faceNodes[inode + 1] == theBetweenNode2) {
+ nbInserted = theNodesToInsert.size();
+
+ // add nodes to insert
+ list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.begin();
+ for (; nIt != theNodesToInsert.end(); nIt++) {
+ poly_nodes.push_back(*nIt);
+ }
+ }
+ }
+ else if (faceNodes[inode] == theBetweenNode2) {
+ if (faceNodes[inode + 1] == theBetweenNode1) {
+ nbInserted = theNodesToInsert.size();
+
+ // add nodes to insert in reversed order
+ list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.end();
+ nIt--;
+ for (; nIt != theNodesToInsert.begin(); nIt--) {
+ poly_nodes.push_back(*nIt);
+ }
+ poly_nodes.push_back(*nIt);
+ }
+ }
+ else {
+ }
+ }
+ }
+ quantities[iface] = nbFaceNodes + nbInserted;
+ }
+
+ // Replace the volume
+ SMESHDS_Mesh *aMesh = GetMeshDS();
+
+ if ( SMDS_MeshElement* newElem = aMesh->AddPolyhedralVolume( poly_nodes, quantities ))
+ {
+ aMesh->SetMeshElementOnShape( newElem, elem->getshapeId() );
+ myLastCreatedElems.push_back( newElem );
+ ReplaceElemInGroups( elem, newElem, aMesh );
+ }
+ aMesh->RemoveElement( elem );
+ }
+}
+
+namespace
+{
+ //================================================================================
+ /*!
+ * \brief Transform any volume into data of SMDSEntity_Polyhedra
+ */
+ //================================================================================
+
+ void volumeToPolyhedron( const SMDS_MeshElement* elem,
+ vector<const SMDS_MeshNode *> & nodes,
+ vector<int> & nbNodeInFaces )
+ {
+ nodes.clear();
+ nbNodeInFaces.clear();
+ SMDS_VolumeTool vTool ( elem );
+ for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+ {
+ const SMDS_MeshNode** fNodes = vTool.GetFaceNodes( iF );
+ nodes.insert( nodes.end(), fNodes, fNodes + vTool.NbFaceNodes( iF ));
+ nbNodeInFaces.push_back( vTool.NbFaceNodes( iF ));
+ }
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Convert elements contained in a sub-mesh to quadratic
+ * \return int - nb of checked elements
+ */
+//=======================================================================
+
+int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm,
+ SMESH_MesherHelper& theHelper,
+ const bool theForce3d)
+{
+ //MESSAGE("convertElemToQuadratic");
+ int nbElem = 0;
+ if( !theSm ) return nbElem;
+
+ vector<int> nbNodeInFaces;
+ vector<const SMDS_MeshNode *> nodes;
+ SMDS_ElemIteratorPtr ElemItr = theSm->GetElements();
+ while(ElemItr->more())
+ {
+ nbElem++;
+ const SMDS_MeshElement* elem = ElemItr->next();
+ if( !elem ) continue;
+
+ // analyse a necessity of conversion
+ const SMDSAbs_ElementType aType = elem->GetType();
+ if ( aType < SMDSAbs_Edge || aType > SMDSAbs_Volume )
+ continue;
+ const SMDSAbs_EntityType aGeomType = elem->GetEntityType();
+ bool hasCentralNodes = false;
+ if ( elem->IsQuadratic() )
+ {
+ bool alreadyOK;
+ switch ( aGeomType ) {
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa:
+ case SMDSEntity_Quad_Penta:
+ alreadyOK = !theHelper.GetIsBiQuadratic(); break;
+
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ case SMDSEntity_TriQuad_Hexa:
+ case SMDSEntity_BiQuad_Penta:
+ alreadyOK = theHelper.GetIsBiQuadratic();
+ hasCentralNodes = true;
+ break;
+ default:
+ alreadyOK = true;
+ }
+ // take into account already present medium nodes
+ switch ( aType ) {
+ case SMDSAbs_Volume:
+ theHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( elem )); break;
+ case SMDSAbs_Face:
+ theHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( elem )); break;
+ case SMDSAbs_Edge:
+ theHelper.AddTLinks( static_cast< const SMDS_MeshEdge* >( elem )); break;
+ default:;
+ }
+ if ( alreadyOK )
+ continue;
+ }
+ // get elem data needed to re-create it
+ //
+ const int id = elem->GetID();
+ const int nbNodes = elem->NbCornerNodes();
+ nodes.assign(elem->begin_nodes(), elem->end_nodes());
+ if ( aGeomType == SMDSEntity_Polyhedra )
+ nbNodeInFaces = static_cast<const SMDS_MeshVolume* >( elem )->GetQuantities();
+ else if ( aGeomType == SMDSEntity_Hexagonal_Prism )
+ volumeToPolyhedron( elem, nodes, nbNodeInFaces );
+
+ // remove a linear element
+ GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
+
+ // remove central nodes of biquadratic elements (biquad->quad conversion)
+ if ( hasCentralNodes )
+ for ( size_t i = nbNodes * 2; i < nodes.size(); ++i )
+ if ( nodes[i]->NbInverseElements() == 0 )
+ GetMeshDS()->RemoveFreeNode( nodes[i], theSm, /*fromGroups=*/true );
+
+ const SMDS_MeshElement* NewElem = 0;
+
+ switch( aType )
+ {
+ case SMDSAbs_Edge :
+ {
+ NewElem = theHelper.AddEdge(nodes[0], nodes[1], id, theForce3d);
+ break;
+ }
+ case SMDSAbs_Face :
+ {
+ switch(nbNodes)
+ {
+ case 3:
+ NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
+ break;
+ case 4:
+ NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ default:
+ NewElem = theHelper.AddPolygonalFace(nodes, id, theForce3d);
+ }
+ break;
+ }
+ case SMDSAbs_Volume :
+ {
+ switch( aGeomType )
+ {
+ case SMDSEntity_Tetra:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ case SMDSEntity_Pyramid:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], id, theForce3d);
+ break;
+ case SMDSEntity_Penta:
+ case SMDSEntity_Quad_Penta:
+ case SMDSEntity_BiQuad_Penta:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], id, theForce3d);
+ break;
+ case SMDSEntity_Hexa:
+ case SMDSEntity_Quad_Hexa:
+ case SMDSEntity_TriQuad_Hexa:
+ NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+ break;
+ case SMDSEntity_Hexagonal_Prism:
+ default:
+ NewElem = theHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
+ }
+ break;
+ }
+ default :
+ continue;
+ }
+ ReplaceElemInGroups( elem, NewElem, GetMeshDS());
+ if( NewElem && NewElem->getshapeId() < 1 )
+ theSm->AddElement( NewElem );
+ }
+ return nbElem;
+}
+//=======================================================================
+//function : ConvertToQuadratic
+//purpose :
+//=======================================================================
+
+void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theToBiQuad)
+{
+ //MESSAGE("ConvertToQuadratic "<< theForce3d << " " << theToBiQuad);
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+
+ SMESH_MesherHelper aHelper(*myMesh);
+
+ aHelper.SetIsQuadratic( true );
+ aHelper.SetIsBiQuadratic( theToBiQuad );
+ aHelper.SetElementsOnShape(true);
+ aHelper.ToFixNodeParameters( true );
+
+ // convert elements assigned to sub-meshes
+ int nbCheckedElems = 0;
+ if ( myMesh->HasShapeToMesh() )
+ {
+ if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
+ {
+ SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
+ while ( smIt->more() ) {
+ SMESH_subMesh* sm = smIt->next();
+ if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() ) {
+ aHelper.SetSubShape( sm->GetSubShape() );
+ nbCheckedElems += convertElemToQuadratic(smDS, aHelper, theForce3d);
+ }
+ }
+ }
+ }
+
+ // convert elements NOT assigned to sub-meshes
+ int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
+ if ( nbCheckedElems < totalNbElems ) // not all elements are in sub-meshes
+ {
+ aHelper.SetElementsOnShape(false);
+ SMESHDS_SubMesh *smDS = 0;
+
+ // convert edges
+ SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator();
+ while( aEdgeItr->more() )
+ {
+ const SMDS_MeshEdge* edge = aEdgeItr->next();
+ if ( !edge->IsQuadratic() )
+ {
+ int id = edge->GetID();
+ const SMDS_MeshNode* n1 = edge->GetNode(0);
+ const SMDS_MeshNode* n2 = edge->GetNode(1);
+
+ meshDS->RemoveFreeElement(edge, smDS, /*fromGroups=*/false);
+
+ const SMDS_MeshEdge* NewEdge = aHelper.AddEdge(n1, n2, id, theForce3d);
+ ReplaceElemInGroups( edge, NewEdge, GetMeshDS());
+ }
+ else
+ {
+ aHelper.AddTLinks( static_cast< const SMDS_MeshEdge* >( edge ));
+ }
+ }
+
+ // convert faces
+ SMDS_FaceIteratorPtr aFaceItr = meshDS->facesIterator();
+ while( aFaceItr->more() )
+ {
+ const SMDS_MeshFace* face = aFaceItr->next();
+ if ( !face ) continue;
+
+ const SMDSAbs_EntityType type = face->GetEntityType();
+ bool alreadyOK;
+ switch( type )
+ {
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ alreadyOK = !theToBiQuad;
+ aHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( face ));
+ break;
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ alreadyOK = theToBiQuad;
+ aHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( face ));
+ break;
+ default: alreadyOK = false;
+ }
+ if ( alreadyOK )
+ continue;
+
+ const int id = face->GetID();
+ vector<const SMDS_MeshNode *> nodes ( face->begin_nodes(), face->end_nodes());
+
+ meshDS->RemoveFreeElement(face, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshFace * NewFace = 0;
+ switch( type )
+ {
+ case SMDSEntity_Triangle:
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_BiQuad_Triangle:
+ NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
+ if ( nodes.size() == 7 && nodes[6]->NbInverseElements() == 0 ) // rm a central node
+ GetMeshDS()->RemoveFreeNode( nodes[6], /*sm=*/0, /*fromGroups=*/true );
+ break;
+
+ case SMDSEntity_Quadrangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ if ( nodes.size() == 9 && nodes[8]->NbInverseElements() == 0 ) // rm a central node
+ GetMeshDS()->RemoveFreeNode( nodes[8], /*sm=*/0, /*fromGroups=*/true );
+ break;
+
+ default:;
+ NewFace = aHelper.AddPolygonalFace(nodes, id, theForce3d);
+ }
+ ReplaceElemInGroups( face, NewFace, GetMeshDS());
+ }
+
+ // convert volumes
+ vector<int> nbNodeInFaces;
+ SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator();
+ while(aVolumeItr->more())
+ {
+ const SMDS_MeshVolume* volume = aVolumeItr->next();
+ if ( !volume ) continue;
+
+ const SMDSAbs_EntityType type = volume->GetEntityType();
+ if ( volume->IsQuadratic() )
+ {
+ bool alreadyOK;
+ switch ( type )
+ {
+ case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break;
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break;
+ case SMDSEntity_Quad_Penta: alreadyOK = !theToBiQuad; break;
+ case SMDSEntity_BiQuad_Penta: alreadyOK = theToBiQuad; break;
+ default: alreadyOK = true;
+ }
+ if ( alreadyOK )
+ {
+ aHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( volume ));
+ continue;
+ }
+ }
+ const int id = volume->GetID();
+ vector<const SMDS_MeshNode *> nodes (volume->begin_nodes(), volume->end_nodes());
+ if ( type == SMDSEntity_Polyhedra )
+ nbNodeInFaces = static_cast<const SMDS_MeshVolume* >(volume)->GetQuantities();
+ else if ( type == SMDSEntity_Hexagonal_Prism )
+ volumeToPolyhedron( volume, nodes, nbNodeInFaces );
+
+ meshDS->RemoveFreeElement(volume, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshVolume * NewVolume = 0;
+ switch ( type )
+ {
+ case SMDSEntity_Tetra:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d );
+ break;
+ case SMDSEntity_Hexa:
+ case SMDSEntity_Quad_Hexa:
+ case SMDSEntity_TriQuad_Hexa:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+ for ( size_t i = 20; i < nodes.size(); ++i ) // rm central nodes
+ if ( nodes[i]->NbInverseElements() == 0 )
+ GetMeshDS()->RemoveFreeNode( nodes[i], /*sm=*/0, /*fromGroups=*/true );
+ break;
+ case SMDSEntity_Pyramid:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+ nodes[3], nodes[4], id, theForce3d);
+ break;
+ case SMDSEntity_Penta:
+ case SMDSEntity_Quad_Penta:
+ case SMDSEntity_BiQuad_Penta:
+ NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+ nodes[3], nodes[4], nodes[5], id, theForce3d);
+ for ( size_t i = 15; i < nodes.size(); ++i ) // rm central nodes
+ if ( nodes[i]->NbInverseElements() == 0 )
+ GetMeshDS()->RemoveFreeNode( nodes[i], /*sm=*/0, /*fromGroups=*/true );
+ break;
+ case SMDSEntity_Hexagonal_Prism:
+ default:
+ NewVolume = aHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
+ }
+ ReplaceElemInGroups(volume, NewVolume, meshDS);
+ }
+ }
+
+ if ( !theForce3d )
+ { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
+ // aHelper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
+ // aHelper.FixQuadraticElements(myError);
+ SMESH_MesherHelper( *myMesh ).FixQuadraticElements(myError);
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Makes given elements quadratic
+ * \param theForce3d - if true, the medium nodes will be placed in the middle of link
+ * \param theElements - elements to make quadratic
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d,
+ TIDSortedElemSet& theElements,
+ const bool theToBiQuad)
+{
+ if ( theElements.empty() ) return;
+
+ // we believe that all theElements are of the same type
+ const SMDSAbs_ElementType elemType = (*theElements.begin())->GetType();
+
+ // get all nodes shared by theElements
+ TIDSortedNodeSet allNodes;
+ TIDSortedElemSet::iterator eIt = theElements.begin();
+ for ( ; eIt != theElements.end(); ++eIt )
+ allNodes.insert( (*eIt)->begin_nodes(), (*eIt)->end_nodes() );
+
+ // complete theElements with elements of lower dim whose all nodes are in allNodes
+
+ TIDSortedElemSet quadAdjacentElems [ SMDSAbs_NbElementTypes ]; // quadratic adjacent elements
+ TIDSortedElemSet checkedAdjacentElems [ SMDSAbs_NbElementTypes ];
+ TIDSortedNodeSet::iterator nIt = allNodes.begin();
+ for ( ; nIt != allNodes.end(); ++nIt )
+ {
+ const SMDS_MeshNode* n = *nIt;
+ SMDS_ElemIteratorPtr invIt = n->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* e = invIt->next();
+ const SMDSAbs_ElementType type = e->GetType();
+ if ( e->IsQuadratic() )
+ {
+ quadAdjacentElems[ type ].insert( e );
+
+ bool alreadyOK;
+ switch ( e->GetEntityType() ) {
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break;
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break;
+ default: alreadyOK = true;
+ }
+ if ( alreadyOK )
+ continue;
+ }
+ if ( type >= elemType )
+ continue; // same type or more complex linear element
+
+ if ( !checkedAdjacentElems[ type ].insert( e ).second )
+ continue; // e is already checked
+
+ // check nodes
+ bool allIn = true;
+ SMDS_NodeIteratorPtr nodeIt = e->nodeIterator();
+ while ( nodeIt->more() && allIn )
+ allIn = allNodes.count( nodeIt->next() );
+ if ( allIn )
+ theElements.insert(e );
+ }
+ }
+
+ SMESH_MesherHelper helper(*myMesh);
+ helper.SetIsQuadratic( true );
+ helper.SetIsBiQuadratic( theToBiQuad );
+
+ // add links of quadratic adjacent elements to the helper
+
+ if ( !quadAdjacentElems[SMDSAbs_Edge].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Edge].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Edge].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshEdge*> (*eIt) );
+ }
+ if ( !quadAdjacentElems[SMDSAbs_Face].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Face].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Face].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshFace*> (*eIt) );
+ }
+ if ( !quadAdjacentElems[SMDSAbs_Volume].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Volume].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Volume].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshVolume*> (*eIt) );
+ }
+
+ // make quadratic (or bi-tri-quadratic) elements instead of linear ones
+
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+ SMESHDS_SubMesh* smDS = 0;
+ for ( eIt = theElements.begin(); eIt != theElements.end(); ++eIt )
+ {
+ const SMDS_MeshElement* elem = *eIt;
+
+ bool alreadyOK;
+ int nbCentralNodes = 0;
+ switch ( elem->GetEntityType() ) {
+ // linear convertible
+ case SMDSEntity_Edge:
+ case SMDSEntity_Triangle:
+ case SMDSEntity_Quadrangle:
+ case SMDSEntity_Tetra:
+ case SMDSEntity_Pyramid:
+ case SMDSEntity_Hexa:
+ case SMDSEntity_Penta: alreadyOK = false; nbCentralNodes = 0; break;
+ // quadratic that can become bi-quadratic
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_Quad_Hexa: alreadyOK =!theToBiQuad; nbCentralNodes = 0; break;
+ // bi-quadratic
+ case SMDSEntity_BiQuad_Triangle:
+ case SMDSEntity_BiQuad_Quadrangle: alreadyOK = theToBiQuad; nbCentralNodes = 1; break;
+ case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; nbCentralNodes = 7; break;
+ // the rest
+ default: alreadyOK = true;
+ }
+ if ( alreadyOK ) continue;
+
+ const SMDSAbs_ElementType type = elem->GetType();
+ const int id = elem->GetID();
+ const int nbNodes = elem->NbCornerNodes();
+ vector<const SMDS_MeshNode *> nodes ( elem->begin_nodes(), elem->end_nodes());
+
+ helper.SetSubShape( elem->getshapeId() );
+
+ if ( !smDS || !smDS->Contains( elem ))
+ smDS = meshDS->MeshElements( elem->getshapeId() );
+ meshDS->RemoveFreeElement(elem, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshElement * newElem = 0;
+ switch( nbNodes )
+ {
+ case 4: // cases for most frequently used element types go first (for optimization)
+ if ( type == SMDSAbs_Volume )
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ else
+ newElem = helper.AddFace (nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ case 8:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+ break;
+ case 3:
+ newElem = helper.AddFace (nodes[0], nodes[1], nodes[2], id, theForce3d);
+ break;
+ case 2:
+ newElem = helper.AddEdge(nodes[0], nodes[1], id, theForce3d);
+ break;
+ case 5:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], id, theForce3d);
+ break;
+ case 6:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], id, theForce3d);
+ break;
+ default:;
+ }
+ ReplaceElemInGroups( elem, newElem, meshDS);
+ if( newElem && smDS )
+ smDS->AddElement( newElem );
+
+ // remove central nodes
+ for ( size_t i = nodes.size() - nbCentralNodes; i < nodes.size(); ++i )
+ if ( nodes[i]->NbInverseElements() == 0 )
+ meshDS->RemoveFreeNode( nodes[i], smDS, /*fromGroups=*/true );
+
+ } // loop on theElements
+
+ if ( !theForce3d )
+ { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
+ // helper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
+ // helper.FixQuadraticElements( myError );
+ SMESH_MesherHelper( *myMesh ).FixQuadraticElements(myError);
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Convert quadratic elements to linear ones and remove quadratic nodes
+ * \return int - nb of checked elements
+ */
+//=======================================================================
+
+int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh * theSm,
+ SMDS_ElemIteratorPtr theItr,
+ const int theShapeID)
+{
+ int nbElem = 0;
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+ ElemFeatures elemType;
+ vector<const SMDS_MeshNode *> nodes;
+
+ while( theItr->more() )
+ {
+ const SMDS_MeshElement* elem = theItr->next();
+ nbElem++;
+ if( elem && elem->IsQuadratic())
+ {
+ // get elem data
+ int nbCornerNodes = elem->NbCornerNodes();
+ nodes.assign( elem->begin_nodes(), elem->end_nodes() );
+
+ elemType.Init( elem, /*basicOnly=*/false ).SetID( elem->GetID() ).SetQuad( false );
+
+ //remove a quadratic element
+ if ( !theSm || !theSm->Contains( elem ))
+ theSm = meshDS->MeshElements( elem->getshapeId() );
+ meshDS->RemoveFreeElement( elem, theSm, /*fromGroups=*/false );
+
+ // remove medium nodes
+ for ( size_t i = nbCornerNodes; i < nodes.size(); ++i )
+ if ( nodes[i]->NbInverseElements() == 0 )
+ meshDS->RemoveFreeNode( nodes[i], theSm );
+
+ // add a linear element
+ nodes.resize( nbCornerNodes );
+ SMDS_MeshElement * newElem = AddElement( nodes, elemType );
+ ReplaceElemInGroups(elem, newElem, meshDS);
+ if( theSm && newElem )
+ theSm->AddElement( newElem );
+ }
+ }
+ return nbElem;
+}
+
+//=======================================================================
+//function : ConvertFromQuadratic
+//purpose :
+//=======================================================================
+
+bool SMESH_MeshEditor::ConvertFromQuadratic()
+{
+ int nbCheckedElems = 0;
+ if ( myMesh->HasShapeToMesh() )
+ {
+ if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
+ {
+ SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
+ while ( smIt->more() ) {
+ SMESH_subMesh* sm = smIt->next();
+ if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() )
+ nbCheckedElems += removeQuadElem( smDS, smDS->GetElements(), sm->GetId() );
+ }
+ }
+ }
+
+ int totalNbElems =
+ GetMeshDS()->NbEdges() + GetMeshDS()->NbFaces() + GetMeshDS()->NbVolumes();
+ if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
+ {
+ SMESHDS_SubMesh *aSM = 0;
+ removeQuadElem( aSM, GetMeshDS()->elementsIterator(), 0 );
+ }
+
+ return true;
+}
+
+namespace
+{
+ //================================================================================
+ /*!
+ * \brief Return true if all medium nodes of the element are in the node set
+ */
+ //================================================================================
+
+ bool allMediumNodesIn(const SMDS_MeshElement* elem, TIDSortedNodeSet& nodeSet )
+ {
+ for ( int i = elem->NbCornerNodes(); i < elem->NbNodes(); ++i )
+ if ( !nodeSet.count( elem->GetNode(i) ))
+ return false;
+ return true;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Makes given elements linear
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
+{
+ if ( theElements.empty() ) return;
+
+ // collect IDs of medium nodes of theElements; some of these nodes will be removed
+ set<int> mediumNodeIDs;
+ TIDSortedElemSet::iterator eIt = theElements.begin();
+ for ( ; eIt != theElements.end(); ++eIt )
+ {
+ const SMDS_MeshElement* e = *eIt;
+ for ( int i = e->NbCornerNodes(); i < e->NbNodes(); ++i )
+ mediumNodeIDs.insert( e->GetNode(i)->GetID() );
+ }
+
+ // replace given elements by linear ones
+ SMDS_ElemIteratorPtr elemIt = SMESHUtils::elemSetIterator( theElements );
+ removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
+
+ // we need to convert remaining elements whose all medium nodes are in mediumNodeIDs
+ // except those elements sharing medium nodes of quadratic element whose medium nodes
+ // are not all in mediumNodeIDs
+
+ // get remaining medium nodes
+ TIDSortedNodeSet mediumNodes;
+ set<int>::iterator nIdsIt = mediumNodeIDs.begin();
+ for ( ; nIdsIt != mediumNodeIDs.end(); ++nIdsIt )
+ if ( const SMDS_MeshNode* n = GetMeshDS()->FindNode( *nIdsIt ))
+ mediumNodes.insert( mediumNodes.end(), n );
+
+ // find more quadratic elements to convert
+ TIDSortedElemSet moreElemsToConvert;
+ TIDSortedNodeSet::iterator nIt = mediumNodes.begin();
+ for ( ; nIt != mediumNodes.end(); ++nIt )
+ {
+ SMDS_ElemIteratorPtr invIt = (*nIt)->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* e = invIt->next();
+ if ( e->IsQuadratic() && allMediumNodesIn( e, mediumNodes ))
+ {
+ // find a more complex element including e and
+ // whose medium nodes are not in mediumNodes
+ bool complexFound = false;
+ for ( int type = e->GetType() + 1; type < SMDSAbs_0DElement; ++type )
+ {
+ SMDS_ElemIteratorPtr invIt2 =
+ (*nIt)->GetInverseElementIterator( SMDSAbs_ElementType( type ));
+ while ( invIt2->more() )
+ {
+ const SMDS_MeshElement* eComplex = invIt2->next();
+ if ( eComplex->IsQuadratic() && !allMediumNodesIn( eComplex, mediumNodes))
+ {
+ int nbCommonNodes = SMESH_MeshAlgos::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 = SMESHUtils::elemSetIterator( moreElemsToConvert );
+ removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
+}
+
+//=======================================================================
+//function : SewSideElements
+//purpose :
+//=======================================================================
+
+SMESH_MeshEditor::Sew_Error
+SMESH_MeshEditor::SewSideElements (TIDSortedElemSet& theSide1,
+ TIDSortedElemSet& theSide2,
+ const SMDS_MeshNode* theFirstNode1,
+ const SMDS_MeshNode* theFirstNode2,
+ const SMDS_MeshNode* theSecondNode1,
+ const SMDS_MeshNode* theSecondNode2)
+{
+ ClearLastCreated();
+
+ if ( theSide1.size() != theSide2.size() )
+ return SEW_DIFF_NB_OF_ELEMENTS;
+
+ Sew_Error aResult = SEW_OK;
+ // Algo:
+ // 1. Build set of faces representing each side
+ // 2. Find which nodes of the side 1 to merge with ones on the side 2
+ // 3. Replace nodes in elements of the side 1 and remove replaced nodes
+
+ // =======================================================================
+ // 1. Build set of faces representing each side:
+ // =======================================================================
+ // a. build set of nodes belonging to faces
+ // b. complete set of faces: find missing faces whose nodes are in set of nodes
+ // c. create temporary faces representing side of volumes if correspondent
+ // face does not exist
+
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+ // TODO algorithm not OK with vtkUnstructuredGrid: 2 meshes can't share nodes
+ //SMDS_Mesh aTmpFacesMesh; // try to use the same mesh
+ TIDSortedElemSet faceSet1, faceSet2;
+ set<const SMDS_MeshElement*> volSet1, volSet2;
+ set<const SMDS_MeshNode*> nodeSet1, nodeSet2;
+ TIDSortedElemSet * faceSetPtr[] = { &faceSet1, &faceSet2 };
+ set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
+ set<const SMDS_MeshNode*> * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
+ TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
+ int iSide, iFace, iNode;
+
+ list<const SMDS_MeshElement* > tempFaceList;
+ for ( iSide = 0; iSide < 2; iSide++ ) {
+ set<const SMDS_MeshNode*> * nodeSet = nodeSetPtr[ iSide ];
+ TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
+ TIDSortedElemSet * faceSet = faceSetPtr[ iSide ];
+ set<const SMDS_MeshElement*> * volSet = volSetPtr [ iSide ];
+ set<const SMDS_MeshElement*>::iterator vIt;
+ TIDSortedElemSet::iterator eIt;
+ set<const SMDS_MeshNode*>::iterator nIt;
+
+ // check that given nodes belong to given elements
+ const SMDS_MeshNode* n1 = ( iSide == 0 ) ? theFirstNode1 : theFirstNode2;
+ const SMDS_MeshNode* n2 = ( iSide == 0 ) ? theSecondNode1 : theSecondNode2;
+ int firstIndex = -1, secondIndex = -1;
+ for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
+ const SMDS_MeshElement* elem = *eIt;
+ if ( firstIndex < 0 ) firstIndex = elem->GetNodeIndex( n1 );
+ if ( secondIndex < 0 ) secondIndex = elem->GetNodeIndex( n2 );
+ if ( firstIndex > -1 && secondIndex > -1 ) break;
+ }
+ if ( firstIndex < 0 || secondIndex < 0 ) {
+ // we can simply return until temporary faces created
+ return (iSide == 0 ) ? SEW_BAD_SIDE1_NODES : SEW_BAD_SIDE2_NODES;
+ }
+
+ // -----------------------------------------------------------
+ // 1a. Collect nodes of existing faces
+ // and build set of face nodes in order to detect missing
+ // faces corresponding to sides of volumes
+ // -----------------------------------------------------------
+
+ set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
+
+ // loop on the given element of a side
+ for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
+ //const SMDS_MeshElement* elem = *eIt;
+ const SMDS_MeshElement* elem = *eIt;
+ if ( elem->GetType() == SMDSAbs_Face ) {
+ faceSet->insert( elem );
+ set <const SMDS_MeshNode*> faceNodeSet;
+ SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+ while ( nodeIt->more() ) {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ nodeSet->insert( n );
+ faceNodeSet.insert( n );
+ }
+ setOfFaceNodeSet.insert( faceNodeSet );
+ }
+ else if ( elem->GetType() == SMDSAbs_Volume )
+ volSet->insert( elem );
+ }
+ // ------------------------------------------------------------------------------
+ // 1b. Complete set of faces: find missing faces whose nodes are in set of nodes
+ // ------------------------------------------------------------------------------
+
+ for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
+ SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() ) { // loop on faces sharing a node
+ const SMDS_MeshElement* f = fIt->next();
+ if ( faceSet->find( f ) == faceSet->end() ) {
+ // check if all nodes are in nodeSet and
+ // complete setOfFaceNodeSet if they are
+ set <const SMDS_MeshNode*> faceNodeSet;
+ SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
+ bool allInSet = true;
+ while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ if ( nodeSet->find( n ) == nodeSet->end() )
+ allInSet = false;
+ else
+ faceNodeSet.insert( n );
+ }
+ if ( allInSet ) {
+ faceSet->insert( f );
+ setOfFaceNodeSet.insert( faceNodeSet );
+ }
+ }
+ }
+ }
+
+ // -------------------------------------------------------------------------
+ // 1c. Create temporary faces representing sides of volumes if correspondent
+ // face does not exist
+ // -------------------------------------------------------------------------
+
+ if ( !volSet->empty() ) {
+ //int nodeSetSize = nodeSet->size();
+
+ // loop on given volumes
+ for ( vIt = volSet->begin(); vIt != volSet->end(); vIt++ ) {
+ SMDS_VolumeTool vol (*vIt);
+ // loop on volume faces: find free faces
+ // --------------------------------------
+ list<const SMDS_MeshElement* > freeFaceList;
+ for ( iFace = 0; iFace < vol.NbFaces(); iFace++ ) {
+ if ( !vol.IsFreeFace( iFace ))
+ continue;
+ // check if there is already a face with same nodes in a face set
+ const SMDS_MeshElement* aFreeFace = 0;
+ const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iFace );
+ int nbNodes = vol.NbFaceNodes( iFace );
+ set <const SMDS_MeshNode*> faceNodeSet;
+ vol.GetFaceNodes( iFace, faceNodeSet );
+ bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
+ if ( isNewFace ) {
+ // no such a face is given but it still can exist, check it
+ vector<const SMDS_MeshNode *> nodes ( fNodes, fNodes + nbNodes);
+ aFreeFace = aMesh->FindElement( nodes, SMDSAbs_Face, /*noMedium=*/false );
+ }
+ if ( !aFreeFace ) {
+ // create a temporary face
+ if ( nbNodes == 3 ) {
+ //aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
+ aFreeFace = aMesh->AddFace( fNodes[0],fNodes[1],fNodes[2] );
+ }
+ else if ( nbNodes == 4 ) {
+ //aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
+ aFreeFace = aMesh->AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
+ }
+ else {
+ vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
+ //aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
+ aFreeFace = aMesh->AddPolygonalFace(poly_nodes);
+ }
+ if ( aFreeFace )
+ tempFaceList.push_back( aFreeFace );
+ }
+
+ if ( aFreeFace )
+ freeFaceList.push_back( aFreeFace );
+
+ } // loop on faces of a volume
+
+ // choose one of several free faces of a volume
+ // --------------------------------------------
+ if ( freeFaceList.size() > 1 ) {
+ // choose a face having max nb of nodes shared by other elems of a side
+ int maxNbNodes = -1;
+ list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
+ while ( fIt != freeFaceList.end() ) { // loop on free faces
+ int nbSharedNodes = 0;
+ SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
+ while ( nodeIt->more() ) { // loop on free face nodes
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
+ while ( invElemIt->more() ) {
+ const SMDS_MeshElement* e = invElemIt->next();
+ nbSharedNodes += faceSet->count( e );
+ nbSharedNodes += elemSet->count( e );
+ }
+ }
+ if ( nbSharedNodes > maxNbNodes ) {
+ maxNbNodes = nbSharedNodes;
+ freeFaceList.erase( freeFaceList.begin(), fIt++ );
+ }
+ else if ( nbSharedNodes == maxNbNodes ) {
+ fIt++;
+ }
+ else {
+ freeFaceList.erase( fIt++ ); // here fIt++ occurs before erase
+ }
+ }
+ if ( freeFaceList.size() > 1 )
+ {
+ // could not choose one face, use another way
+ // choose a face most close to the bary center of the opposite side
+ gp_XYZ aBC( 0., 0., 0. );
+ set <const SMDS_MeshNode*> addedNodes;
+ TIDSortedElemSet * elemSet2 = elemSetPtr[ 1 - iSide ];
+ eIt = elemSet2->begin();
+ for ( eIt = elemSet2->begin(); eIt != elemSet2->end(); eIt++ ) {
+ SMDS_ElemIteratorPtr nodeIt = (*eIt)->nodesIterator();
+ while ( nodeIt->more() ) { // loop on free face nodes
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ if ( addedNodes.insert( n ).second )
+ aBC += gp_XYZ( n->X(),n->Y(),n->Z() );
+ }
+ }
+ aBC /= addedNodes.size();
+ double minDist = DBL_MAX;
+ fIt = freeFaceList.begin();
+ while ( fIt != freeFaceList.end() ) { // loop on free faces
+ double dist = 0;
+ SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
+ while ( nodeIt->more() ) { // loop on free face nodes
+ const SMDS_MeshNode* n =
+ static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ gp_XYZ p( n->X(),n->Y(),n->Z() );
+ dist += ( aBC - p ).SquareModulus();
+ }
+ if ( dist < minDist ) {
+ minDist = dist;
+ freeFaceList.erase( freeFaceList.begin(), fIt++ );
+ }
+ else
+ fIt = freeFaceList.erase( fIt++ );
+ }
+ }
+ } // choose one of several free faces of a volume
+
+ if ( freeFaceList.size() == 1 ) {
+ const SMDS_MeshElement* aFreeFace = freeFaceList.front();
+ faceSet->insert( aFreeFace );
+ // complete a node set with nodes of a found free face
+ // for ( iNode = 0; iNode < ; iNode++ )
+ // nodeSet->insert( fNodes[ iNode ] );
+ }
+
+ } // loop on volumes of a side
+
+ // // complete a set of faces if new nodes in a nodeSet appeared
+ // // ----------------------------------------------------------
+ // if ( nodeSetSize != nodeSet->size() ) {
+ // for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
+ // SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
+ // while ( fIt->more() ) { // loop on faces sharing a node
+ // const SMDS_MeshElement* f = fIt->next();
+ // if ( faceSet->find( f ) == faceSet->end() ) {
+ // // check if all nodes are in nodeSet and
+ // // complete setOfFaceNodeSet if they are
+ // set <const SMDS_MeshNode*> faceNodeSet;
+ // SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
+ // bool allInSet = true;
+ // while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
+ // const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ // if ( nodeSet->find( n ) == nodeSet->end() )
+ // allInSet = false;
+ // else
+ // faceNodeSet.insert( n );
+ // }
+ // if ( allInSet ) {
+ // faceSet->insert( f );
+ // setOfFaceNodeSet.insert( faceNodeSet );
+ // }
+ // }
+ // }
+ // }
+ // }
+ } // Create temporary faces, if there are volumes given
+ } // loop on sides
+
+ if ( faceSet1.size() != faceSet2.size() ) {
+ // delete temporary faces: they are in reverseElements of actual nodes
+ // SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
+ // while ( tmpFaceIt->more() )
+ // aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
+ // list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
+ // for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
+ // aMesh->RemoveElement(*tmpFaceIt);
+ MESSAGE("Diff nb of faces");
+ return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ }
+
+ // ============================================================
+ // 2. Find nodes to merge:
+ // bind a node to remove to a node to put instead
+ // ============================================================
+
+ TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
+ if ( theFirstNode1 != theFirstNode2 )
+ nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
+ if ( theSecondNode1 != theSecondNode2 )
+ nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
+
+ LinkID_Gen aLinkID_Gen( GetMeshDS() );
+ set< long > linkIdSet; // links to process
+ linkIdSet.insert( aLinkID_Gen.GetLinkID( theFirstNode1, theSecondNode1 ));
+
+ typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > NLink;
+ list< NLink > linkList[2];
+ linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
+ linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
+ // loop on links in linkList; find faces by links and append links
+ // of the found faces to linkList
+ list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
+ for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ )
+ {
+ NLink link[] = { *linkIt[0], *linkIt[1] };
+ long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
+ if ( !linkIdSet.count( linkID ) )
+ continue;
+
+ // by links, find faces in the face sets,
+ // and find indices of link nodes in the found faces;
+ // in a face set, there is only one or no face sharing a link
+ // ---------------------------------------------------------------
+
+ const SMDS_MeshElement* face[] = { 0, 0 };
+ vector<const SMDS_MeshNode*> fnodes[2];
+ int iLinkNode[2][2];
+ TIDSortedElemSet avoidSet;
+ for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
+ const SMDS_MeshNode* n1 = link[iSide].first;
+ const SMDS_MeshNode* n2 = link[iSide].second;
+ //cout << "Side " << iSide << " ";
+ //cout << "L( " << n1->GetID() << ", " << n2->GetID() << " ) " << endl;
+ // find a face by two link nodes
+ face[ iSide ] = SMESH_MeshAlgos::FindFaceInSet( n1, n2,
+ *faceSetPtr[ iSide ], avoidSet,
+ &iLinkNode[iSide][0],
+ &iLinkNode[iSide][1] );
+ if ( face[ iSide ])
+ {
+ //cout << " F " << face[ iSide]->GetID() <<endl;
+ faceSetPtr[ iSide ]->erase( face[ iSide ]);
+ // put face nodes to fnodes
+ SMDS_MeshElement::iterator nIt( face[ iSide ]->interlacedNodesIterator() ), nEnd;
+ fnodes[ iSide ].assign( nIt, nEnd );
+ fnodes[ iSide ].push_back( fnodes[ iSide ].front());
+ }
+ }
+
+ // check similarity of elements of the sides
+ if (aResult == SEW_OK && (( face[0] && !face[1] ) || ( !face[0] && face[1] ))) {
+ MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
+ if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
+ aResult = ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
+ }
+ else {
+ aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ }
+ break; // do not return because it's necessary to remove tmp faces
+ }
+
+ // set nodes to merge
+ // -------------------
+
+ if ( face[0] && face[1] ) {
+ const int nbNodes = face[0]->NbNodes();
+ if ( nbNodes != face[1]->NbNodes() ) {
+ MESSAGE("Diff nb of face nodes");
+ aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ break; // do not return because it s necessary to remove tmp faces
+ }
+ bool reverse[] = { false, false }; // order of nodes in the link
+ for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
+ // analyse link orientation in faces
+ int i1 = iLinkNode[ iSide ][ 0 ];
+ int i2 = iLinkNode[ iSide ][ 1 ];
+ reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
+ }
+ int di1 = reverse[0] ? -1 : +1, i1 = iLinkNode[0][1] + di1;
+ int di2 = reverse[1] ? -1 : +1, i2 = iLinkNode[1][1] + di2;
+ for ( int i = nbNodes - 2; i > 0; --i, i1 += di1, i2 += di2 )
+ {
+ nReplaceMap.insert ( make_pair ( fnodes[0][ ( i1 + nbNodes ) % nbNodes ],
+ fnodes[1][ ( i2 + nbNodes ) % nbNodes ]));
+ }
+
+ // add other links of the faces to linkList
+ // -----------------------------------------
+
+ for ( iNode = 0; iNode < nbNodes; iNode++ ) {
+ linkID = aLinkID_Gen.GetLinkID( fnodes[0][iNode], fnodes[0][iNode+1] );
+ pair< set<long>::iterator, bool > iter_isnew = linkIdSet.insert( linkID );
+ if ( !iter_isnew.second ) { // already in a set: no need to process
+ linkIdSet.erase( iter_isnew.first );
+ }
+ else // new in set == encountered for the first time: add
+ {
+ const SMDS_MeshNode* n1 = fnodes[0][ iNode ];
+ const SMDS_MeshNode* n2 = fnodes[0][ iNode + 1];
+ linkList[0].push_back ( NLink( n1, n2 ));
+ linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
+ }
+ }
+ } // 2 faces found
+
+ if ( faceSetPtr[0]->empty() || faceSetPtr[1]->empty() )
+ break;
+
+ } // loop on link lists
+
+ if ( aResult == SEW_OK &&
+ ( //linkIt[0] != linkList[0].end() ||
+ !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
+ MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
+ " " << (faceSetPtr[1]->empty()));
+ aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ }
+
+ // ====================================================================
+ // 3. Replace nodes in elements of the side 1 and remove replaced nodes
+ // ====================================================================
+
+ // delete temporary faces
+ // SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
+ // while ( tmpFaceIt->more() )
+ // aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
+ list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
+ for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
+ aMesh->RemoveElement(*tmpFaceIt);
+
+ if ( aResult != SEW_OK)
+ return aResult;
+
+ list< int > nodeIDsToRemove;
+ vector< const SMDS_MeshNode*> nodes;
+ ElemFeatures elemType;
+
+ // loop on nodes replacement map
+ TNodeNodeMap::iterator nReplaceMapIt = nReplaceMap.begin(), nnIt;
+ for ( ; nReplaceMapIt != nReplaceMap.end(); nReplaceMapIt++ )
+ if ( (*nReplaceMapIt).first != (*nReplaceMapIt).second )
+ {
+ const SMDS_MeshNode* nToRemove = (*nReplaceMapIt).first;
+ nodeIDsToRemove.push_back( nToRemove->GetID() );
+ // loop on elements sharing nToRemove
+ SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
+ while ( invElemIt->more() ) {
+ const SMDS_MeshElement* e = invElemIt->next();
+ // get a new suite of nodes: make replacement
+ int nbReplaced = 0, i = 0, nbNodes = e->NbNodes();
+ nodes.resize( nbNodes );
+ SMDS_ElemIteratorPtr nIt = e->nodesIterator();
+ while ( nIt->more() ) {
+ const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
+ nnIt = nReplaceMap.find( n );
+ if ( nnIt != nReplaceMap.end() ) {
+ nbReplaced++;
+ n = (*nnIt).second;
+ }
+ nodes[ i++ ] = n;
+ }
+ // if ( nbReplaced == nbNodes && e->GetType() == SMDSAbs_Face )
+ // elemIDsToRemove.push_back( e->GetID() );
+ // else
+ if ( nbReplaced )
+ {
+ elemType.Init( e, /*basicOnly=*/false ).SetID( e->GetID() );
+ aMesh->RemoveElement( e );
+
+ if ( SMDS_MeshElement* newElem = this->AddElement( nodes, elemType ))
+ {
+ AddToSameGroups( newElem, e, aMesh );
+ if ( int aShapeId = e->getshapeId() )
+ aMesh->SetMeshElementOnShape( newElem, aShapeId );
+ }
+ }
+ }
+ }
+
+ Remove( nodeIDsToRemove, true );
+
+ return aResult;
+}
+
+//================================================================================
+/*!
+ * \brief Find corresponding nodes in two sets of faces
+ * \param theSide1 - first face set
+ * \param theSide2 - second first face
+ * \param theFirstNode1 - a boundary node of set 1
+ * \param theFirstNode2 - a node of set 2 corresponding to theFirstNode1
+ * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
+ * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
+ * \param nReplaceMap - output map of corresponding nodes
+ * \return bool - is a success or not
+ */
+//================================================================================
+
+#ifdef _DEBUG_
+//#define DEBUG_MATCHING_NODES
+#endif
+
+SMESH_MeshEditor::Sew_Error
+SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
+ set<const SMDS_MeshElement*>& theSide2,
+ const SMDS_MeshNode* theFirstNode1,
+ const SMDS_MeshNode* theFirstNode2,
+ const SMDS_MeshNode* theSecondNode1,
+ const SMDS_MeshNode* theSecondNode2,
+ TNodeNodeMap & nReplaceMap)
+{
+ set<const SMDS_MeshElement*> * faceSetPtr[] = { &theSide1, &theSide2 };
+
+ nReplaceMap.clear();
+ if ( theFirstNode1 != theFirstNode2 )
+ nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
+ if ( theSecondNode1 != theSecondNode2 )
+ nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
+
+ set< SMESH_TLink > linkSet; // set of nodes where order of nodes is ignored
+ linkSet.insert( SMESH_TLink( theFirstNode1, theSecondNode1 ));
+
+ list< NLink > linkList[2];
+ linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
+ linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
+
+ // loop on links in linkList; find faces by links and append links
+ // of the found faces to linkList
+ list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
+ for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
+ NLink link[] = { *linkIt[0], *linkIt[1] };
+ if ( linkSet.find( link[0] ) == linkSet.end() )
+ continue;
+
+ // by links, find faces in the face sets,
+ // and find indices of link nodes in the found faces;
+ // in a face set, there is only one or no face sharing a link
+ // ---------------------------------------------------------------
+
+ const SMDS_MeshElement* face[] = { 0, 0 };
+ list<const SMDS_MeshNode*> notLinkNodes[2];
+ //bool reverse[] = { false, false }; // order of notLinkNodes
+ int nbNodes[2];
+ for ( int iSide = 0; iSide < 2; iSide++ ) // loop on 2 sides
+ {
+ const SMDS_MeshNode* n1 = link[iSide].first;
+ const SMDS_MeshNode* n2 = link[iSide].second;
+ set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
+ set< const SMDS_MeshElement* > facesOfNode1;
+ for ( int iNode = 0; iNode < 2; iNode++ ) // loop on 2 nodes of a link
+ {
+ // during a loop of the first node, we find all faces around n1,
+ // during a loop of the second node, we find one face sharing both n1 and n2
+ const SMDS_MeshNode* n = iNode ? n1 : n2; // a node of a link
+ SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() ) { // loop on faces sharing a node
+ const SMDS_MeshElement* f = fIt->next();
+ if (faceSet->find( f ) != faceSet->end() && // f is in face set
+ ! facesOfNode1.insert( f ).second ) // f encounters twice
+ {
+ if ( face[ iSide ] ) {
+ MESSAGE( "2 faces per link " );
+ return ( iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
+ }
+ face[ iSide ] = f;
+ faceSet->erase( f );
+
+ // get not link nodes
+ int nbN = f->NbNodes();
+ if ( f->IsQuadratic() )
+ nbN /= 2;
+ nbNodes[ iSide ] = nbN;
+ list< const SMDS_MeshNode* > & nodes = notLinkNodes[ iSide ];
+ int i1 = f->GetNodeIndex( n1 );
+ int i2 = f->GetNodeIndex( n2 );
+ int iEnd = nbN, iBeg = -1, iDelta = 1;
+ bool reverse = ( Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1 );
+ if ( reverse ) {
+ std::swap( iEnd, iBeg ); iDelta = -1;
+ }
+ int i = i2;
+ while ( true ) {
+ i += iDelta;
+ if ( i == iEnd ) i = iBeg + iDelta;
+ if ( i == i1 ) break;
+ nodes.push_back ( f->GetNode( i ) );
+ }
+ }
+ }
+ }
+ }
+ // check similarity of elements of the sides
+ if (( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
+ MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
+ if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
+ return ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
+ }
+ else {
+ return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ }
+ }
+
+ // set nodes to merge
+ // -------------------
+
+ if ( face[0] && face[1] ) {
+ if ( nbNodes[0] != nbNodes[1] ) {
+ MESSAGE("Diff nb of face nodes");
+ return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
+ }
+#ifdef DEBUG_MATCHING_NODES
+ MESSAGE ( " Link 1: " << link[0].first->GetID() <<" "<< link[0].second->GetID()
+ << " F 1: " << face[0] << "| Link 2: " << link[1].first->GetID() <<" "
+ << link[1].second->GetID() << " F 2: " << face[1] << " | Bind: " ) ;
+#endif
+ int nbN = nbNodes[0];
+ {
+ list<const SMDS_MeshNode*>::iterator n1 = notLinkNodes[0].begin();
+ list<const SMDS_MeshNode*>::iterator n2 = notLinkNodes[1].begin();
+ for ( int i = 0 ; i < nbN - 2; ++i ) {
+#ifdef DEBUG_MATCHING_NODES
+ MESSAGE ( (*n1)->GetID() << " to " << (*n2)->GetID() );
+#endif
+ nReplaceMap.insert( make_pair( *(n1++), *(n2++) ));
+ }
+ }
+
+ // add other links of the face 1 to linkList
+ // -----------------------------------------
+
+ const SMDS_MeshElement* f0 = face[0];
+ const SMDS_MeshNode* n1 = f0->GetNode( nbN - 1 );
+ for ( int i = 0; i < nbN; i++ )
+ {
+ const SMDS_MeshNode* n2 = f0->GetNode( i );
+ pair< set< SMESH_TLink >::iterator, bool > iter_isnew =
+ linkSet.insert( SMESH_TLink( n1, n2 ));
+ if ( !iter_isnew.second ) { // already in a set: no need to process
+ linkSet.erase( iter_isnew.first );
+ }
+ else // new in set == encountered for the first time: add
+ {
+#ifdef DEBUG_MATCHING_NODES
+ MESSAGE ( "Add link 1: " << n1->GetID() << " " << n2->GetID() << " "
+ << " | link 2: " << nReplaceMap[n1]->GetID() << " " << nReplaceMap[n2]->GetID() << " " );
+#endif
+ linkList[0].push_back ( NLink( n1, n2 ));
+ linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
+ }
+ n1 = n2;
+ }
+ } // 2 faces found
+ } // loop on link lists
+
+ return SEW_OK;
+}
+
+namespace // automatically find theAffectedElems for DoubleNodes()
+{
+ bool isOut( const SMDS_MeshNode* n, const gp_XYZ& norm, const SMDS_MeshElement* elem );
+
+ //--------------------------------------------------------------------------------
+ // Nodes shared by adjacent FissureBorder's.
+ // 1 node if FissureBorder separates faces
+ // 2 nodes if FissureBorder separates volumes
+ struct SubBorder
+ {
+ const SMDS_MeshNode* _nodes[2];
+ int _nbNodes;
+
+ SubBorder( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2 = 0 )
+ {
+ _nodes[0] = n1;
+ _nodes[1] = n2;
+ _nbNodes = bool( n1 ) + bool( n2 );
+ if ( _nbNodes == 2 && n1 > n2 )
+ std::swap( _nodes[0], _nodes[1] );
+ }
+ bool operator<( const SubBorder& other ) const
+ {
+ for ( int i = 0; i < _nbNodes; ++i )
+ {
+ if ( _nodes[i] < other._nodes[i] ) return true;
+ if ( _nodes[i] > other._nodes[i] ) return false;
+ }
+ return false;
+ }
+ };
+
+ //--------------------------------------------------------------------------------
+ // Map a SubBorder to all FissureBorder it bounds
+ struct FissureBorder;
+ typedef std::map< SubBorder, std::vector< FissureBorder* > > TBorderLinks;
+ typedef TBorderLinks::iterator TMappedSub;
+
+ //--------------------------------------------------------------------------------
+ /*!
+ * \brief Element border (volume facet or face edge) at a fissure
+ */
+ struct FissureBorder
+ {
+ std::vector< const SMDS_MeshNode* > _nodes; // border nodes
+ const SMDS_MeshElement* _elems[2]; // volume or face adjacent to fissure
+
+ std::vector< TMappedSub > _mappedSubs; // Sub() in TBorderLinks map
+ std::vector< const SMDS_MeshNode* > _sortedNodes; // to compare FissureBorder's
+
+ FissureBorder( FissureBorder && from ) // move constructor
+ {
+ std::swap( _nodes, from._nodes );
+ std::swap( _sortedNodes, from._sortedNodes );
+ _elems[0] = from._elems[0];
+ _elems[1] = from._elems[1];
+ }
+
+ FissureBorder( const SMDS_MeshElement* elemToDuplicate,
+ std::vector< const SMDS_MeshElement* > & adjElems)
+ : _nodes( elemToDuplicate->NbCornerNodes() )
+ {
+ for ( size_t i = 0; i < _nodes.size(); ++i )
+ _nodes[i] = elemToDuplicate->GetNode( i );
+
+ SMDSAbs_ElementType type = SMDSAbs_ElementType( elemToDuplicate->GetType() + 1 );
+ findAdjacent( type, adjElems );
+ }
+
+ FissureBorder( const SMDS_MeshNode** nodes,
+ const size_t nbNodes,
+ const SMDSAbs_ElementType adjElemsType,
+ std::vector< const SMDS_MeshElement* > & adjElems)
+ : _nodes( nodes, nodes + nbNodes )
+ {
+ findAdjacent( adjElemsType, adjElems );
+ }
+
+ void findAdjacent( const SMDSAbs_ElementType adjElemsType,
+ std::vector< const SMDS_MeshElement* > & adjElems)
+ {
+ _elems[0] = _elems[1] = 0;
+ adjElems.clear();
+ if ( SMDS_Mesh::GetElementsByNodes( _nodes, adjElems, adjElemsType ))
+ for ( size_t i = 0; i < adjElems.size() && i < 2; ++i )
+ _elems[i] = adjElems[i];
+ }
+
+ bool operator<( const FissureBorder& other ) const
+ {
+ return GetSortedNodes() < other.GetSortedNodes();
+ }
+
+ const std::vector< const SMDS_MeshNode* >& GetSortedNodes() const
+ {
+ if ( _sortedNodes.empty() && !_nodes.empty() )
+ {
+ FissureBorder* me = const_cast<FissureBorder*>( this );
+ me->_sortedNodes = me->_nodes;
+ std::sort( me->_sortedNodes.begin(), me->_sortedNodes.end() );
+ }
+ return _sortedNodes;
+ }
+
+ size_t NbSub() const
+ {
+ return _nodes.size();
+ }
+
+ SubBorder Sub(size_t i) const
+ {
+ return SubBorder( _nodes[i], NbSub() > 2 ? _nodes[ (i+1)%NbSub() ] : 0 );
+ }
+
+ void AddSelfTo( TBorderLinks& borderLinks )
+ {
+ _mappedSubs.resize( NbSub() );
+ for ( size_t i = 0; i < NbSub(); ++i )
+ {
+ TBorderLinks::iterator s2b =
+ borderLinks.insert( std::make_pair( Sub(i), TBorderLinks::mapped_type() )).first;
+ s2b->second.push_back( this );
+ _mappedSubs[ i ] = s2b;
+ }
+ }
+
+ void Clear()
+ {
+ _nodes.clear();
+ }
+
+ const SMDS_MeshElement* GetMarkedElem() const
+ {
+ if ( _nodes.empty() ) return 0; // cleared
+ if ( _elems[0] && _elems[0]->isMarked() ) return _elems[0];
+ if ( _elems[1] && _elems[1]->isMarked() ) return _elems[1];
+ return 0;
+ }
+
+ gp_XYZ GetNorm() const // normal to the border
+ {
+ gp_XYZ norm;
+ if ( _nodes.size() == 2 )
+ {
+ gp_XYZ avgNorm( 0,0,0 ); // sum of normals of adjacent faces
+ if ( SMESH_MeshAlgos::FaceNormal( _elems[0], norm ))
+ avgNorm += norm;
+ if ( SMESH_MeshAlgos::FaceNormal( _elems[1], norm ))
+ avgNorm += norm;
+
+ gp_XYZ bordDir( SMESH_NodeXYZ( _nodes[0] ) - SMESH_NodeXYZ( _nodes[1] ));
+ norm = bordDir ^ avgNorm;
+ }
+ else
+ {
+ SMESH_NodeXYZ p0( _nodes[0] );
+ SMESH_NodeXYZ p1( _nodes[1] );
+ SMESH_NodeXYZ p2( _nodes[2] );
+ norm = ( p0 - p1 ) ^ ( p2 - p1 );
+ }
+ if ( isOut( _nodes[0], norm, GetMarkedElem() ))
+ norm.Reverse();
+
+ return norm;
+ }
+
+ void ChooseSide() // mark an _elem located at positive side of fissure
+ {
+ _elems[0]->setIsMarked( true );
+ gp_XYZ norm = GetNorm();
+ double maxX = norm.Coord(1);
+ if ( Abs( maxX ) < Abs( norm.Coord(2)) ) maxX = norm.Coord(2);
+ if ( Abs( maxX ) < Abs( norm.Coord(3)) ) maxX = norm.Coord(3);
+ if ( maxX < 0 )
+ {
+ _elems[0]->setIsMarked( false );
+ _elems[1]->setIsMarked( true );
+ }
+ }
+
+ }; // struct FissureBorder
+
+ //--------------------------------------------------------------------------------
+ /*!
+ * \brief Classifier of elements at fissure edge
+ */
+ class FissureNormal
+ {
+ std::vector< gp_XYZ > _normals;
+ bool _bothIn;
+
+ public:
+ void Add( const SMDS_MeshNode* n, const FissureBorder& bord )
+ {
+ _bothIn = false;
+ _normals.reserve(2);
+ _normals.push_back( bord.GetNorm() );
+ if ( _normals.size() == 2 )
+ _bothIn = !isOut( n, _normals[0], bord.GetMarkedElem() );
+ }
+
+ bool IsIn( const SMDS_MeshNode* n, const SMDS_MeshElement* elem ) const
+ {
+ bool isIn = false;
+ switch ( _normals.size() ) {
+ case 1:
+ {
+ isIn = !isOut( n, _normals[0], elem );
+ break;
+ }
+ case 2:
+ {
+ bool in1 = !isOut( n, _normals[0], elem );
+ bool in2 = !isOut( n, _normals[1], elem );
+ isIn = _bothIn ? ( in1 && in2 ) : ( in1 || in2 );
+ }
+ }
+ return isIn;
+ }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Classify an element by a plane passing through a node
+ */
+ //================================================================================
+
+ bool isOut( const SMDS_MeshNode* n, const gp_XYZ& norm, const SMDS_MeshElement* elem )
+ {
+ SMESH_NodeXYZ p = n;
+ double sumDot = 0;
+ for ( int i = 0, nb = elem->NbCornerNodes(); i < nb; ++i )
+ {
+ SMESH_NodeXYZ pi = elem->GetNode( i );
+ sumDot += norm * ( pi - p );
+ }
+ return sumDot < -1e-100;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Find FissureBorder's by nodes to duplicate
+ */
+ //================================================================================
+
+ void findFissureBorders( const TIDSortedElemSet& theNodes,
+ std::vector< FissureBorder > & theFissureBorders )
+ {
+ TIDSortedElemSet::const_iterator nIt = theNodes.begin();
+ const SMDS_MeshNode* n = dynamic_cast< const SMDS_MeshNode*>( *nIt );
+ if ( !n ) return;
+ SMDSAbs_ElementType elemType = SMDSAbs_Volume;
+ if ( n->NbInverseElements( elemType ) == 0 )
+ {
+ elemType = SMDSAbs_Face;
+ if ( n->NbInverseElements( elemType ) == 0 )
+ return;
+ }
+ // unmark elements touching the fissure
+ for ( ; nIt != theNodes.end(); ++nIt )
+ SMESH_MeshAlgos::MarkElems( cast2Node(*nIt)->GetInverseElementIterator(), false );
+
+ // loop on elements touching the fissure to get their borders belonging to the fissure
+ std::set< FissureBorder > fissureBorders;
+ std::vector< const SMDS_MeshElement* > adjElems;
+ std::vector< const SMDS_MeshNode* > nodes;
+ SMDS_VolumeTool volTool;
+ for ( nIt = theNodes.begin(); nIt != theNodes.end(); ++nIt )
+ {
+ SMDS_ElemIteratorPtr invIt = cast2Node(*nIt)->GetInverseElementIterator( elemType );
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* eInv = invIt->next();
+ if ( eInv->isMarked() ) continue;
+ eInv->setIsMarked( true );
+
+ if ( elemType == SMDSAbs_Volume )
+ {
+ volTool.Set( eInv );
+ int iQuad = eInv->IsQuadratic() ? 2 : 1;
+ for ( int iF = 0, nbF = volTool.NbFaces(); iF < nbF; ++iF )
+ {
+ const SMDS_MeshNode** nn = volTool.GetFaceNodes( iF );
+ int nbN = volTool.NbFaceNodes( iF ) / iQuad;
+ nodes.clear();
+ bool allOnFissure = true;
+ for ( int iN = 0; iN < nbN && allOnFissure; iN += iQuad )
+ if (( allOnFissure = theNodes.count( nn[ iN ])))
+ nodes.push_back( nn[ iN ]);
+ if ( allOnFissure )
+ fissureBorders.insert( std::move( FissureBorder( &nodes[0], nodes.size(),
+ elemType, adjElems )));
+ }
+ }
+ else // elemType == SMDSAbs_Face
+ {
+ const SMDS_MeshNode* nn[2] = { eInv->GetNode( eInv->NbCornerNodes()-1 ), 0 };
+ bool onFissure0 = theNodes.count( nn[0] ), onFissure1;
+ for ( int iN = 0, nbN = eInv->NbCornerNodes(); iN < nbN; ++iN )
+ {
+ nn[1] = eInv->GetNode( iN );
+ onFissure1 = theNodes.count( nn[1] );
+ if ( onFissure0 && onFissure1 )
+ fissureBorders.insert( std::move( FissureBorder( nn, 2, elemType, adjElems )));
+ nn[0] = nn[1];
+ onFissure0 = onFissure1;
+ }
+ }
+ }
+ }
+
+ theFissureBorders.reserve( theFissureBorders.size() + fissureBorders.size());
+ std::set< FissureBorder >::iterator bord = fissureBorders.begin();
+ for ( ; bord != fissureBorders.end(); ++bord )
+ {
+ theFissureBorders.push_back( std::move( const_cast<FissureBorder&>( *bord ) ));
+ }
+ return;
+ } // findFissureBorders()
+
+ //================================================================================
+ /*!
+ * \brief Find elements on one side of a fissure defined by elements or nodes to duplicate
+ * \param [in] theElemsOrNodes - elements or nodes to duplicate
+ * \param [in] theNodesNot - nodes not to duplicate
+ * \param [out] theAffectedElems - the found elements
+ */
+ //================================================================================
+
+ void findAffectedElems( const TIDSortedElemSet& theElemsOrNodes,
+ TIDSortedElemSet& theAffectedElems)
+ {
+ if ( theElemsOrNodes.empty() ) return;
+
+ // find FissureBorder's
+
+ std::vector< FissureBorder > fissure;
+ std::vector< const SMDS_MeshElement* > elemsByFacet;
+
+ TIDSortedElemSet::const_iterator elIt = theElemsOrNodes.begin();
+ if ( (*elIt)->GetType() == SMDSAbs_Node )
+ {
+ findFissureBorders( theElemsOrNodes, fissure );
+ }
+ else
+ {
+ fissure.reserve( theElemsOrNodes.size() );
+ for ( ; elIt != theElemsOrNodes.end(); ++elIt )
+ fissure.push_back( std::move( FissureBorder( *elIt, elemsByFacet )));
+ }
+ if ( fissure.empty() )
+ return;
+
+ // fill borderLinks
+
+ TBorderLinks borderLinks;
+
+ for ( size_t i = 0; i < fissure.size(); ++i )
+ {
+ fissure[i].AddSelfTo( borderLinks );
+ }
+
+ // get theAffectedElems
+
+ // unmark elements having nodes on the fissure, theAffectedElems elements will be marked
+ for ( size_t i = 0; i < fissure.size(); ++i )
+ for ( size_t j = 0; j < fissure[i]._nodes.size(); ++j )
+ {
+ SMESH_MeshAlgos::MarkElemNodes( fissure[i]._nodes[j]->GetInverseElementIterator(),
+ false, /*markElem=*/true );
+ }
+
+ std::vector<const SMDS_MeshNode *> facetNodes;
+ std::map< const SMDS_MeshNode*, FissureNormal > fissEdgeNodes2Norm;
+ boost::container::flat_set< const SMDS_MeshNode* > fissureNodes;
+
+ // choose a side of fissure
+ fissure[0].ChooseSide();
+ theAffectedElems.insert( fissure[0].GetMarkedElem() );
+
+ size_t nbCheckedBorders = 0;
+ while ( nbCheckedBorders < fissure.size() )
+ {
+ // find a FissureBorder to treat
+ FissureBorder* bord = 0;
+ for ( size_t i = 0; i < fissure.size() && !bord; ++i )
+ if ( fissure[i].GetMarkedElem() )
+ bord = & fissure[i];
+ for ( size_t i = 0; i < fissure.size() && !bord; ++i )
+ if ( fissure[i].NbSub() > 0 && fissure[i]._elems[0] )
+ {
+ bord = & fissure[i];
+ bord->ChooseSide();
+ theAffectedElems.insert( bord->GetMarkedElem() );
+ }
+ if ( !bord ) return;
+ ++nbCheckedBorders;
+
+ // treat FissureBorder's linked to bord
+ fissureNodes.clear();
+ fissureNodes.insert( bord->_nodes.begin(), bord->_nodes.end() );
+ for ( size_t i = 0; i < bord->NbSub(); ++i )
+ {
+ TBorderLinks::iterator l2b = bord->_mappedSubs[ i ];
+ if ( l2b == borderLinks.end() || l2b->second.empty() ) continue;
+ std::vector< FissureBorder* >& linkedBorders = l2b->second;
+ const SubBorder& sb = l2b->first;
+ const SMDS_MeshElement* bordElem = bord->GetMarkedElem();
+
+ if ( linkedBorders.size() == 1 ) // fissure edge reached, fill fissEdgeNodes2Norm
+ {
+ for ( int j = 0; j < sb._nbNodes; ++j )
+ fissEdgeNodes2Norm[ sb._nodes[j] ].Add( sb._nodes[j], *bord );
+ continue;
+ }
+
+ // add to theAffectedElems elems sharing nodes of a SubBorder and a node of bordElem
+ // until an elem adjacent to a neighbour FissureBorder is found
+ facetNodes.clear();
+ facetNodes.insert( facetNodes.end(), sb._nodes, sb._nodes + sb._nbNodes );
+ facetNodes.resize( sb._nbNodes + 1 );
+
+ while ( bordElem )
+ {
+ // check if bordElem is adjacent to a neighbour FissureBorder
+ for ( size_t j = 0; j < linkedBorders.size(); ++j )
+ {
+ FissureBorder* bord2 = linkedBorders[j];
+ if ( bord2 == bord ) continue;
+ if ( bordElem == bord2->_elems[0] || bordElem == bord2->_elems[1] )
+ bordElem = 0;
+ else
+ fissureNodes.insert( bord2->_nodes.begin(), bord2->_nodes.end() );
+ }
+ if ( !bordElem )
+ break;
+
+ // find the next bordElem
+ const SMDS_MeshElement* nextBordElem = 0;
+ for ( int iN = 0, nbN = bordElem->NbCornerNodes(); iN < nbN && !nextBordElem; ++iN )
+ {
+ const SMDS_MeshNode* n = bordElem->GetNode( iN );
+ if ( fissureNodes.count( n )) continue;
+
+ facetNodes[ sb._nbNodes ] = n;
+ elemsByFacet.clear();
+ if ( SMDS_Mesh::GetElementsByNodes( facetNodes, elemsByFacet ) > 1 )
+ {
+ for ( size_t iE = 0; iE < elemsByFacet.size(); ++iE )
+ if ( elemsByFacet[ iE ] != bordElem &&
+ !elemsByFacet[ iE ]->isMarked() )
+ {
+ theAffectedElems.insert( elemsByFacet[ iE ]);
+ elemsByFacet[ iE ]->setIsMarked( true );
+ if ( elemsByFacet[ iE ]->GetType() == bordElem->GetType() )
+ nextBordElem = elemsByFacet[ iE ];
+ }
+ }
+ }
+ bordElem = nextBordElem;
+
+ } // while ( bordElem )
+
+ linkedBorders.clear(); // not to treat this link any more
+
+ } // loop on SubBorder's of a FissureBorder
+
+ bord->Clear();
+
+ } // loop on FissureBorder's
+
+
+ // add elements sharing only one node of the fissure, except those sharing fissure edge nodes
+
+ // mark nodes of theAffectedElems
+ SMESH_MeshAlgos::MarkElemNodes( theAffectedElems.begin(), theAffectedElems.end(), true );
+
+ // unmark nodes of the fissure
+ elIt = theElemsOrNodes.begin();
+ if ( (*elIt)->GetType() == SMDSAbs_Node )
+ SMESH_MeshAlgos::MarkElems( elIt, theElemsOrNodes.end(), false );
+ else
+ SMESH_MeshAlgos::MarkElemNodes( elIt, theElemsOrNodes.end(), false );
+
+ std::vector< gp_XYZ > normVec;
+
+ // loop on nodes of the fissure, add elements having marked nodes
+ for ( elIt = theElemsOrNodes.begin(); elIt != theElemsOrNodes.end(); ++elIt )
+ {
+ const SMDS_MeshElement* e = (*elIt);
+ if ( e->GetType() != SMDSAbs_Node )
+ e->setIsMarked( true ); // avoid adding a fissure element
+
+ for ( int iN = 0, nbN = e->NbCornerNodes(); iN < nbN; ++iN )
+ {
+ const SMDS_MeshNode* n = e->GetNode( iN );
+ if ( fissEdgeNodes2Norm.count( n ))
+ continue;
+
+ SMDS_ElemIteratorPtr invIt = n->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* eInv = invIt->next();
+ if ( eInv->isMarked() ) continue;
+ eInv->setIsMarked( true );
+
+ SMDS_ElemIteratorPtr nIt = eInv->nodesIterator();
+ while( nIt->more() )
+ if ( nIt->next()->isMarked())
+ {
+ theAffectedElems.insert( eInv );
+ SMESH_MeshAlgos::MarkElems( eInv->nodesIterator(), true );
+ n->setIsMarked( false );
+ break;
+ }
+ }
+ }
+ }
+
+ // add elements on the fissure edge
+ std::map< const SMDS_MeshNode*, FissureNormal >::iterator n2N;
+ for ( n2N = fissEdgeNodes2Norm.begin(); n2N != fissEdgeNodes2Norm.end(); ++n2N )
+ {
+ const SMDS_MeshNode* edgeNode = n2N->first;
+ const FissureNormal & normals = n2N->second;
+
+ SMDS_ElemIteratorPtr invIt = edgeNode->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* eInv = invIt->next();
+ if ( eInv->isMarked() ) continue;
+ eInv->setIsMarked( true );
+
+ // classify eInv using normals
+ bool toAdd = normals.IsIn( edgeNode, eInv );
+ if ( toAdd ) // check if all nodes lie on the fissure edge
+ {
+ bool notOnEdge = false;
+ for ( int iN = 0, nbN = eInv->NbCornerNodes(); iN < nbN && !notOnEdge; ++iN )
+ notOnEdge = !fissEdgeNodes2Norm.count( eInv->GetNode( iN ));
+ toAdd = notOnEdge;
+ }
+ if ( toAdd )
+ {
+ theAffectedElems.insert( eInv );
+ }
+ }
+ }
+
+ return;
+ } // findAffectedElems()
+} // namespace
+
+//================================================================================
+/*!
+ * \brief Create elements equal (on same nodes) to given ones
+ * \param [in] theElements - a set of elems to duplicate. If it is empty, all
+ * elements of the uppest dimension are duplicated.
+ */
+//================================================================================
+
+void SMESH_MeshEditor::DoubleElements( const TIDSortedElemSet& theElements )
+{
+ ClearLastCreated();
+ SMESHDS_Mesh* mesh = GetMeshDS();
+
+ // get an element type and an iterator over elements
+
+ SMDSAbs_ElementType type = SMDSAbs_All;
+ SMDS_ElemIteratorPtr elemIt;
+ if ( theElements.empty() )
+ {
+ if ( mesh->NbNodes() == 0 )
+ return;
+ // get most complex type
+ SMDSAbs_ElementType types[SMDSAbs_NbElementTypes] = {
+ SMDSAbs_Volume, SMDSAbs_Face, SMDSAbs_Edge,
+ SMDSAbs_0DElement, SMDSAbs_Ball, SMDSAbs_Node
+ };
+ for ( int i = 0; i < SMDSAbs_NbElementTypes; ++i )
+ if ( mesh->GetMeshInfo().NbElements( types[i] ))
+ {
+ type = types[i];
+ elemIt = mesh->elementsIterator( type );
+ break;
+ }
+ }
+ else
+ {
+ type = (*theElements.begin())->GetType();
+ elemIt = SMESHUtils::elemSetIterator( theElements );
+ }
+
+ // un-mark all elements to avoid duplicating just created elements
+ SMESH_MeshAlgos::MarkElems( mesh->elementsIterator( type ), false );
+
+ // duplicate elements
+
+ ElemFeatures elemType;
+
+ vector< const SMDS_MeshNode* > nodes;
+ while ( elemIt->more() )
+ {
+ const SMDS_MeshElement* elem = elemIt->next();
+ if ( elem->GetType() != type || elem->isMarked() )
+ continue;
+
+ elemType.Init( elem, /*basicOnly=*/false );
+ nodes.assign( elem->begin_nodes(), elem->end_nodes() );
+
+ if ( const SMDS_MeshElement* newElem = AddElement( nodes, elemType ))
+ newElem->setIsMarked( true );
+ }
+}
+
+//================================================================================
+/*!
+ \brief Creates a hole in a mesh by doubling the nodes of some particular elements
+ \param theElems - the list of elements (edges or faces) to be replicated
+ The nodes for duplication could be found from these elements
+ \param theNodesNot - list of nodes to NOT replicate
+ \param theAffectedElems - the list of elements (cells and edges) to which the
+ replicated nodes should be associated to.
+ \return TRUE if operation has been completed successfully, FALSE otherwise
+*/
+//================================================================================
+
+bool SMESH_MeshEditor::DoubleNodes( const TIDSortedElemSet& theElems,
+ const TIDSortedElemSet& theNodesNot,
+ const TIDSortedElemSet& theAffectedElems )
+{
+ ClearLastCreated();
+
+ if ( theElems.size() == 0 )
+ return false;
+
+ SMESHDS_Mesh* aMeshDS = GetMeshDS();
+ if ( !aMeshDS )
+ return false;
+
+ bool res = false;
+ TNodeNodeMap anOldNodeToNewNode;
+ // duplicate elements and nodes
+ res = doubleNodes( aMeshDS, theElems, theNodesNot, anOldNodeToNewNode, true );
+ // replce nodes by duplications
+ res = doubleNodes( aMeshDS, theAffectedElems, theNodesNot, anOldNodeToNewNode, false );
+ return res;
+}
+
+//================================================================================
+/*!
+ \brief Creates a hole in a mesh by doubling the nodes of some particular elements
+ \param theMeshDS - mesh instance
+ \param theElems - the elements replicated or modified (nodes should be changed)
+ \param theNodesNot - nodes to NOT replicate
+ \param theNodeNodeMap - relation of old node to new created node
+ \param theIsDoubleElem - flag os to replicate element or modify
+ \return TRUE if operation has been completed successfully, FALSE otherwise
+*/
+//================================================================================
+
+bool SMESH_MeshEditor::doubleNodes(SMESHDS_Mesh* theMeshDS,
+ const TIDSortedElemSet& theElems,
+ const TIDSortedElemSet& theNodesNot,
+ TNodeNodeMap& theNodeNodeMap,
+ const bool theIsDoubleElem )
+{
+ // iterate through element and duplicate them (by nodes duplication)
+ bool res = false;
+ std::vector<const SMDS_MeshNode*> newNodes;
+ ElemFeatures elemType;
+
+ TIDSortedElemSet::const_iterator elemItr = theElems.begin();
+ for ( ; elemItr != theElems.end(); ++elemItr )
+ {
+ const SMDS_MeshElement* anElem = *elemItr;
+ // if (!anElem)
+ // continue;
+
+ // duplicate nodes to duplicate element
+ bool isDuplicate = false;
+ newNodes.resize( anElem->NbNodes() );
+ SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
+ int ind = 0;
+ while ( anIter->more() )
+ {
+ const SMDS_MeshNode* aCurrNode = static_cast<const SMDS_MeshNode*>( anIter->next() );
+ const SMDS_MeshNode* aNewNode = aCurrNode;
+ TNodeNodeMap::iterator n2n = theNodeNodeMap.find( aCurrNode );
+ if ( n2n != theNodeNodeMap.end() )
+ {
+ aNewNode = n2n->second;
+ }
+ else if ( theIsDoubleElem && !theNodesNot.count( aCurrNode ))
+ {
+ // duplicate node
+ aNewNode = theMeshDS->AddNode( aCurrNode->X(), aCurrNode->Y(), aCurrNode->Z() );
+ copyPosition( aCurrNode, aNewNode );
+ theNodeNodeMap[ aCurrNode ] = aNewNode;
+ myLastCreatedNodes.push_back( aNewNode );
+ }
+ isDuplicate |= (aCurrNode != aNewNode);
+ newNodes[ ind++ ] = aNewNode;
+ }
+ if ( !isDuplicate )
+ continue;
+
+ if ( theIsDoubleElem )
+ AddElement( newNodes, elemType.Init( anElem, /*basicOnly=*/false ));
+ else
+ theMeshDS->ChangeElementNodes( anElem, &newNodes[ 0 ], newNodes.size() );
+
+ res = true;
+ }
+ return res;
+}
+
+//================================================================================
+/*!
+ \brief Creates a hole in a mesh by doubling the nodes of some particular elements
+ \param theNodes - identifiers of nodes to be doubled
+ \param theModifiedElems - identifiers of elements to be updated by the new (doubled)
+ nodes. If list of element identifiers is empty then nodes are doubled but
+ they not assigned to elements
+ \return TRUE if operation has been completed successfully, FALSE otherwise
+*/
+//================================================================================
+
+bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
+ const std::list< int >& theListOfModifiedElems )
+{
+ ClearLastCreated();
+
+ if ( theListOfNodes.size() == 0 )
+ return false;
+
+ SMESHDS_Mesh* aMeshDS = GetMeshDS();
+ if ( !aMeshDS )
+ return false;
+
+ // iterate through nodes and duplicate them
+
+ std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > anOldNodeToNewNode;
+
+ std::list< int >::const_iterator aNodeIter;
+ for ( aNodeIter = theListOfNodes.begin(); aNodeIter != theListOfNodes.end(); ++aNodeIter )
+ {
+ const SMDS_MeshNode* aNode = aMeshDS->FindNode( *aNodeIter );
+ if ( !aNode )
+ continue;
+
+ // duplicate node
+
+ const SMDS_MeshNode* aNewNode = aMeshDS->AddNode( aNode->X(), aNode->Y(), aNode->Z() );
+ if ( aNewNode )
+ {
+ copyPosition( aNode, aNewNode );
+ anOldNodeToNewNode[ aNode ] = aNewNode;
+ myLastCreatedNodes.push_back( aNewNode );
+ }
+ }
+
+ // Change nodes of elements
+
+ std::vector<const SMDS_MeshNode*> aNodeArr;
+
+ std::list< int >::const_iterator anElemIter;
+ for ( anElemIter = theListOfModifiedElems.begin();
+ anElemIter != theListOfModifiedElems.end();
+ anElemIter++ )
+ {
+ const SMDS_MeshElement* anElem = aMeshDS->FindElement( *anElemIter );
+ if ( !anElem )
+ continue;
+
+ aNodeArr.assign( anElem->begin_nodes(), anElem->end_nodes() );
+ for( size_t i = 0; i < aNodeArr.size(); ++i )
+ {
+ std::map< const SMDS_MeshNode*, const SMDS_MeshNode* >::iterator n2n =
+ anOldNodeToNewNode.find( aNodeArr[ i ]);
+ if ( n2n != anOldNodeToNewNode.end() )
+ aNodeArr[ i ] = n2n->second;
+ }
+ aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], aNodeArr.size() );
+ }
+
+ return true;
+}
+
+namespace {
+
+ //================================================================================
+ /*!
+ \brief Check if element located inside shape
+ \return TRUE if IN or ON shape, FALSE otherwise
+ */
+ //================================================================================
+
+ template<class Classifier>
+ bool isInside(const SMDS_MeshElement* theElem,
+ Classifier& theClassifier,
+ const double theTol)
+ {
+ gp_XYZ centerXYZ (0, 0, 0);
+ for ( SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator(); aNodeItr->more(); )
+ centerXYZ += SMESH_NodeXYZ( aNodeItr->next() );
+
+ gp_Pnt aPnt = centerXYZ / theElem->NbNodes();
+ theClassifier.Perform(aPnt, theTol);
+ TopAbs_State aState = theClassifier.State();
+ return (aState == TopAbs_IN || aState == TopAbs_ON );
+ }
+
+ //================================================================================
+ /*!
+ * \brief Classifier of the 3D point on the TopoDS_Face
+ * with interaface suitable for isInside()
+ */
+ //================================================================================
+
+ struct _FaceClassifier
+ {
+ Extrema_ExtPS _extremum;
+ BRepAdaptor_Surface _surface;
+ TopAbs_State _state;
+
+ _FaceClassifier(const TopoDS_Face& face):_extremum(),_surface(face),_state(TopAbs_OUT)
+ {
+ _extremum.Initialize( _surface,
+ _surface.FirstUParameter(), _surface.LastUParameter(),
+ _surface.FirstVParameter(), _surface.LastVParameter(),
+ _surface.Tolerance(), _surface.Tolerance() );
+ }
+ void Perform(const gp_Pnt& aPnt, double theTol)
+ {
+ theTol *= theTol;
+ _state = TopAbs_OUT;
+ _extremum.Perform(aPnt);
+ if ( _extremum.IsDone() )
+ for ( int iSol = 1; iSol <= _extremum.NbExt() && _state == TopAbs_OUT; ++iSol)
+ _state = ( _extremum.SquareDistance(iSol) <= theTol ? TopAbs_IN : TopAbs_OUT );
+ }
+ TopAbs_State State() const
+ {
+ return _state;
+ }
+ };
+}
+
+//================================================================================
+/*!
+ \brief Identify the elements that will be affected by node duplication (actual duplication is not performed).
+ This method is the first step of DoubleNodeElemGroupsInRegion.
+ \param theElems - list of groups of elements (edges or faces) to be replicated
+ \param theNodesNot - list of groups of nodes not to replicated
+ \param theShape - shape to detect affected elements (element which geometric center
+ located on or inside shape). If the shape is null, detection is done on faces orientations
+ (select elements with a gravity center on the side given by faces normals).
+ This mode (null shape) is faster, but works only when theElems are faces, with coherents orientations.
+ The replicated nodes should be associated to affected elements.
+ \return true
+ \sa DoubleNodeElemGroupsInRegion()
+*/
+//================================================================================
+
+bool SMESH_MeshEditor::AffectedElemGroupsInRegion( const TIDSortedElemSet& theElems,
+ const TIDSortedElemSet& theNodesNot,
+ const TopoDS_Shape& theShape,
+ TIDSortedElemSet& theAffectedElems)
+{
+ if ( theShape.IsNull() )
+ {
+ findAffectedElems( theElems, theAffectedElems );
+ }
+ else
+ {
+ const double aTol = Precision::Confusion();
+ std::unique_ptr< BRepClass3d_SolidClassifier> bsc3d;
+ std::unique_ptr<_FaceClassifier> aFaceClassifier;
+ if ( theShape.ShapeType() == TopAbs_SOLID )
+ {
+ bsc3d.reset( new BRepClass3d_SolidClassifier(theShape));;
+ bsc3d->PerformInfinitePoint(aTol);
+ }
+ else if (theShape.ShapeType() == TopAbs_FACE )
+ {
+ aFaceClassifier.reset( new _FaceClassifier(TopoDS::Face(theShape)));
+ }
+
+ // iterates on indicated elements and get elements by back references from their nodes
+ TIDSortedElemSet::const_iterator elemItr = theElems.begin();
+ for ( ; elemItr != theElems.end(); ++elemItr )
+ {
+ SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
+ SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
+ while ( nodeItr->more() )
+ {
+ const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
+ if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
+ continue;
+ SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
+ while ( backElemItr->more() )
+ {
+ const SMDS_MeshElement* curElem = backElemItr->next();
+ if ( curElem && theElems.find(curElem) == theElems.end() &&
+ ( bsc3d.get() ?
+ isInside( curElem, *bsc3d, aTol ) :
+ isInside( curElem, *aFaceClassifier, aTol )))
+ theAffectedElems.insert( curElem );
+ }
+ }
+ }
+ }
+ return true;
+}
+
+//================================================================================
+/*!
+ \brief Creates a hole in a mesh by doubling the nodes of some particular elements
+ \param theElems - group of of elements (edges or faces) to be replicated
+ \param theNodesNot - group of nodes not to replicate
+ \param theShape - shape to detect affected elements (element which geometric center
+ located on or inside shape).
+ The replicated nodes should be associated to affected elements.
+ \return TRUE if operation has been completed successfully, FALSE otherwise
+*/
+//================================================================================
+
+bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems,
+ const TIDSortedElemSet& theNodesNot,
+ const TopoDS_Shape& theShape )
+{
+ if ( theShape.IsNull() )
+ return false;
+
+ const double aTol = Precision::Confusion();
+ SMESHUtils::Deleter< BRepClass3d_SolidClassifier> bsc3d;
+ SMESHUtils::Deleter<_FaceClassifier> aFaceClassifier;
+ if ( theShape.ShapeType() == TopAbs_SOLID )
+ {
+ bsc3d._obj = new BRepClass3d_SolidClassifier( theShape );
+ bsc3d->PerformInfinitePoint(aTol);
+ }
+ else if (theShape.ShapeType() == TopAbs_FACE )
+ {
+ aFaceClassifier._obj = new _FaceClassifier( TopoDS::Face( theShape ));
+ }
+
+ // iterates on indicated elements and get elements by back references from their nodes
+ TIDSortedElemSet anAffected;
+ TIDSortedElemSet::const_iterator elemItr = theElems.begin();
+ for ( ; elemItr != theElems.end(); ++elemItr )
+ {
+ SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
+ if (!anElem)
+ continue;
+
+ SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
+ while ( nodeItr->more() )
+ {
+ const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
+ if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
+ continue;
+ SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
+ while ( backElemItr->more() )
+ {
+ const SMDS_MeshElement* curElem = backElemItr->next();
+ if ( curElem && theElems.find(curElem) == theElems.end() &&
+ ( bsc3d ?
+ isInside( curElem, *bsc3d, aTol ) :
+ isInside( curElem, *aFaceClassifier, aTol )))
+ anAffected.insert( curElem );
+ }
+ }
+ }
+ return DoubleNodes( theElems, theNodesNot, anAffected );
+}
+
+/*!
+ * \brief compute an oriented angle between two planes defined by four points.
+ * The vector (p0,p1) defines the intersection of the 2 planes (p0,p1,g1) and (p0,p1,g2)
+ * @param p0 base of the rotation axe
+ * @param p1 extremity of the rotation axe
+ * @param g1 belongs to the first plane
+ * @param g2 belongs to the second plane
+ */
+double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const gp_Pnt& g1, const gp_Pnt& g2)
+{
+ gp_Vec vref(p0, p1);
+ gp_Vec v1(p0, g1);
+ gp_Vec v2(p0, g2);
+ gp_Vec n1 = vref.Crossed(v1);
+ gp_Vec n2 = vref.Crossed(v2);
+ try {
+ return n2.AngleWithRef(n1, vref);
+ }
+ catch ( Standard_Failure ) {
+ }
+ return Max( v1.Magnitude(), v2.Magnitude() );
+}
+
+/*!
+ * \brief Double nodes on shared faces between groups of volumes and create flat elements on demand.
+ * The list of groups must contain at least two groups. The groups have to be disjoint: no common element into two different groups.
+ * The nodes of the internal faces at the boundaries of the groups are doubled. Optionally, the internal faces are replaced by flat elements.
+ * Triangles are transformed into prisms, and quadrangles into hexahedrons.
+ * The flat elements are stored in groups of volumes. These groups are named according to the position of the group in the list:
+ * the group j_n_p is the group of the flat elements that are built between the group #n and the group #p in the list.
+ * If there is no shared faces between the group #n and the group #p in the list, the group j_n_p is not created.
+ * All the flat elements are gathered into the group named "joints3D" (or "joints2D" in 2D situation).
+ * The flat element of the multiple junctions between the simple junction are stored in a group named "jointsMultiples".
+ * \param theElems - list of groups of volumes, where a group of volume is a set of
+ * SMDS_MeshElements sorted by Id.
+ * \param createJointElems - if TRUE, create the elements
+ * \param onAllBoundaries - if TRUE, the nodes and elements are also created on
+ * the boundary between \a theDomains and the rest mesh
+ * \return TRUE if operation has been completed successfully, FALSE otherwise
+ */
+bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSortedElemSet>& theElems,
+ bool createJointElems,
+ bool onAllBoundaries)
+{
+ // MESSAGE("----------------------------------------------");
+ // MESSAGE("SMESH_MeshEditor::doubleNodesOnGroupBoundaries");
+ // MESSAGE("----------------------------------------------");
+
+ SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
+ meshDS->BuildDownWardConnectivity(true);
+ CHRONO(50);
+ SMDS_UnstructuredGrid *grid = meshDS->GetGrid();
+
+ // --- build the list of faces shared by 2 domains (group of elements), with their domain and volume indexes
+ // build the list of cells with only a node or an edge on the border, with their domain and volume indexes
+ // build the list of nodes shared by 2 or more domains, with their domain indexes
+
+ std::map<DownIdType, std::map<int,int>, DownIdCompare> faceDomains; // face --> (id domain --> id volume)
+ std::map<int,int>celldom; // cell vtkId --> domain
+ std::map<DownIdType, std::map<int,int>, DownIdCompare> cellDomains; // oldNode --> (id domain --> id cell)
+ std::map<int, std::map<int,int> > nodeDomains; // oldId --> (domainId --> newId)
+ faceDomains.clear();
+ celldom.clear();
+ cellDomains.clear();
+ nodeDomains.clear();
+ std::map<int,int> emptyMap;
+ std::set<int> emptySet;
+ emptyMap.clear();
+
+ //MESSAGE(".. Number of domains :"<<theElems.size());
+
+ TIDSortedElemSet theRestDomElems;
+ const int iRestDom = -1;
+ const int idom0 = onAllBoundaries ? iRestDom : 0;
+ const int nbDomains = theElems.size();
+
+ // Check if the domains do not share an element
+ for (int idom = 0; idom < nbDomains-1; idom++)
+ {
+ // MESSAGE("... Check of domain #" << idom);
+ const TIDSortedElemSet& domain = theElems[idom];
+ TIDSortedElemSet::const_iterator elemItr = domain.begin();
+ for (; elemItr != domain.end(); ++elemItr)
+ {
+ const SMDS_MeshElement* anElem = *elemItr;
+ int idombisdeb = idom + 1 ;
+ // check if the element belongs to a domain further in the list
+ for ( size_t idombis = idombisdeb; idombis < theElems.size(); idombis++ )
+ {
+ const TIDSortedElemSet& domainbis = theElems[idombis];
+ if ( domainbis.count( anElem ))
+ {
+ MESSAGE(".... Domain #" << idom);
+ MESSAGE(".... Domain #" << idombis);
+ throw SALOME_Exception("The domains are not disjoint.");
+ return false ;
+ }
+ }
+ }
+ }
+
+ for (int idom = 0; idom < nbDomains; idom++)
+ {
+
+ // --- build a map (face to duplicate --> volume to modify)
+ // with all the faces shared by 2 domains (group of elements)
+ // and corresponding volume of this domain, for each shared face.
+ // a volume has a face shared by 2 domains if it has a neighbor which is not in his domain.
+
+ //MESSAGE("... Neighbors of domain #" << idom);
+ const TIDSortedElemSet& domain = theElems[idom];
+ TIDSortedElemSet::const_iterator elemItr = domain.begin();
+ for (; elemItr != domain.end(); ++elemItr)
+ {
+ const SMDS_MeshElement* anElem = *elemItr;
+ if (!anElem)
+ continue;
+ int vtkId = anElem->GetVtkID();
+ //MESSAGE(" vtkId " << vtkId << " smdsId " << anElem->GetID());
+ int neighborsVtkIds[NBMAXNEIGHBORS];
+ int downIds[NBMAXNEIGHBORS];
+ unsigned char downTypes[NBMAXNEIGHBORS];
+ int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId);
+ for (int n = 0; n < nbNeighbors; n++)
+ {
+ int smdsId = meshDS->FromVtkToSmds(neighborsVtkIds[n]);
+ const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
+ if (elem && ! domain.count(elem)) // neighbor is in another domain : face is shared
+ {
+ bool ok = false;
+ for ( size_t idombis = 0; idombis < theElems.size() && !ok; idombis++) // check if the neighbor belongs to another domain of the list
+ {
+ // MESSAGE("Domain " << idombis);
+ const TIDSortedElemSet& domainbis = theElems[idombis];
+ if ( domainbis.count(elem)) ok = true ; // neighbor is in a correct domain : face is kept
+ }
+ if ( ok || onAllBoundaries ) // the characteristics of the face is stored
+ {
+ DownIdType face(downIds[n], downTypes[n]);
+ if (!faceDomains[face].count(idom))
+ {
+ faceDomains[face][idom] = vtkId; // volume associated to face in this domain
+ celldom[vtkId] = idom;
+ //MESSAGE(" cell with a border " << vtkId << " domain " << idom);
+ }
+ if ( !ok )
+ {
+ theRestDomElems.insert( elem );
+ faceDomains[face][iRestDom] = neighborsVtkIds[n];
+ celldom[neighborsVtkIds[n]] = iRestDom;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ //MESSAGE("Number of shared faces " << faceDomains.size());
+ std::map<DownIdType, std::map<int, int>, DownIdCompare>::iterator itface;
+
+ // --- explore the shared faces domain by domain,
+ // explore the nodes of the face and see if they belong to a cell in the domain,
+ // which has only a node or an edge on the border (not a shared face)
+
+ for (int idomain = idom0; idomain < nbDomains; idomain++)
+ {
+ //MESSAGE("Domain " << idomain);
+ const TIDSortedElemSet& domain = (idomain == iRestDom) ? theRestDomElems : theElems[idomain];
+ itface = faceDomains.begin();
+ for (; itface != faceDomains.end(); ++itface)
+ {
+ const std::map<int, int>& domvol = itface->second;
+ if (!domvol.count(idomain))
+ continue;
+ DownIdType face = itface->first;
+ //MESSAGE(" --- face " << face.cellId);
+ std::set<int> oldNodes;
+ oldNodes.clear();
+ grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+ std::set<int>::iterator itn = oldNodes.begin();
+ for (; itn != oldNodes.end(); ++itn)
+ {
+ int oldId = *itn;
+ //MESSAGE(" node " << oldId);
+ vtkCellLinks::Link l = grid->GetCellLinks()->GetLink(oldId);
+ for (int i=0; i<l.ncells; i++)
+ {
+ int vtkId = l.cells[i];
+ const SMDS_MeshElement* anElem = GetMeshDS()->FindElement(GetMeshDS()->FromVtkToSmds(vtkId));
+ if (!domain.count(anElem))
+ continue;
+ int vtkType = grid->GetCellType(vtkId);
+ int downId = grid->CellIdToDownId(vtkId);
+ if (downId < 0)
+ {
+ MESSAGE("doubleNodesOnGroupBoundaries: internal algorithm problem");
+ continue; // not OK at this stage of the algorithm:
+ //no cells created after BuildDownWardConnectivity
+ }
+ DownIdType aCell(downId, vtkType);
+ cellDomains[aCell][idomain] = vtkId;
+ celldom[vtkId] = idomain;
+ //MESSAGE(" cell " << vtkId << " domain " << idomain);
+ }
+ }
+ }
+ }
+
+ // --- explore the shared faces domain by domain, to duplicate the nodes in a coherent way
+ // for each shared face, get the nodes
+ // for each node, for each domain of the face, create a clone of the node
+
+ // --- edges at the intersection of 3 or 4 domains, with the order of domains to build
+ // junction elements of type prism or hexa. the key is the pair of nodesId (lower first)
+ // the value is the ordered domain ids. (more than 4 domains not taken into account)
+
+ std::map<std::vector<int>, std::vector<int> > edgesMultiDomains; // nodes of edge --> ordered domains
+ std::map<int, std::vector<int> > mutipleNodes; // nodes multi domains with domain order
+ std::map<int, std::vector<int> > mutipleNodesToFace; // nodes multi domains with domain order to transform in Face (junction between 3 or more 2D domains)
+
+ //MESSAGE(".. Duplication of the nodes");
+ for (int idomain = idom0; idomain < nbDomains; idomain++)
+ {
+ itface = faceDomains.begin();
+ for (; itface != faceDomains.end(); ++itface)
+ {
+ const std::map<int, int>& domvol = itface->second;
+ if (!domvol.count(idomain))
+ continue;
+ DownIdType face = itface->first;
+ //MESSAGE(" --- face " << face.cellId);
+ std::set<int> oldNodes;
+ oldNodes.clear();
+ grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+ std::set<int>::iterator itn = oldNodes.begin();
+ for (; itn != oldNodes.end(); ++itn)
+ {
+ int oldId = *itn;
+ if (nodeDomains[oldId].empty())
+ {
+ nodeDomains[oldId][idomain] = oldId; // keep the old node in the first domain
+ //MESSAGE("-+-+-b oldNode " << oldId << " domain " << idomain);
+ }
+ std::map<int, int>::const_iterator itdom = domvol.begin();
+ for (; itdom != domvol.end(); ++itdom)
+ {
+ int idom = itdom->first;
+ //MESSAGE(" domain " << idom);
+ if (!nodeDomains[oldId].count(idom)) // --- node to clone
+ {
+ if (nodeDomains[oldId].size() >= 2) // a multiple node
+ {
+ vector<int> orderedDoms;
+ //MESSAGE("multiple node " << oldId);
+ if (mutipleNodes.count(oldId))
+ orderedDoms = mutipleNodes[oldId];
+ else
+ {
+ map<int,int>::iterator it = nodeDomains[oldId].begin();
+ for (; it != nodeDomains[oldId].end(); ++it)
+ orderedDoms.push_back(it->first);
+ }
+ orderedDoms.push_back(idom); // TODO order ==> push_front or back
+ //stringstream txt;
+ //for (int i=0; i<orderedDoms.size(); i++)
+ // txt << orderedDoms[i] << " ";
+ //MESSAGE("orderedDoms " << txt.str());
+ mutipleNodes[oldId] = orderedDoms;
+ }
+ double *coords = grid->GetPoint(oldId);
+ SMDS_MeshNode *newNode = meshDS->AddNode(coords[0], coords[1], coords[2]);
+ copyPosition( meshDS->FindNodeVtk( oldId ), newNode );
+ int newId = newNode->GetVtkID();
+ nodeDomains[oldId][idom] = newId; // cloned node for other domains
+ //MESSAGE("-+-+-c oldNode " << oldId << " domain " << idomain << " newNode " << newId << " domain " << idom << " size=" <<nodeDomains[oldId].size());
+ }
+ }
+ }
+ }
+ }
+
+ //MESSAGE(".. Creation of elements");
+ for (int idomain = idom0; idomain < nbDomains; idomain++)
+ {
+ itface = faceDomains.begin();
+ for (; itface != faceDomains.end(); ++itface)
+ {
+ std::map<int, int> domvol = itface->second;
+ if (!domvol.count(idomain))
+ continue;
+ DownIdType face = itface->first;
+ //MESSAGE(" --- face " << face.cellId);
+ std::set<int> oldNodes;
+ oldNodes.clear();
+ grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+ int nbMultipleNodes = 0;
+ std::set<int>::iterator itn = oldNodes.begin();
+ for (; itn != oldNodes.end(); ++itn)
+ {
+ int oldId = *itn;
+ if (mutipleNodes.count(oldId))
+ nbMultipleNodes++;
+ }
+ if (nbMultipleNodes > 1) // check if an edge of the face is shared between 3 or more domains
+ {
+ //MESSAGE("multiple Nodes detected on a shared face");
+ int downId = itface->first.cellId;
+ unsigned char cellType = itface->first.cellType;
+ // --- shared edge or shared face ?
+ if ((cellType == VTK_LINE) || (cellType == VTK_QUADRATIC_EDGE)) // shared edge (between two faces)
+ {
+ int nodes[3];
+ int nbNodes = grid->getDownArray(cellType)->getNodes(downId, nodes);
+ for (int i=0; i< nbNodes; i=i+nbNodes-1) // i=0 , i=nbNodes-1
+ if (mutipleNodes.count(nodes[i]))
+ if (!mutipleNodesToFace.count(nodes[i]))
+ mutipleNodesToFace[nodes[i]] = mutipleNodes[nodes[i]];
+ }
+ else // shared face (between two volumes)
+ {
+ int nbEdges = grid->getDownArray(cellType)->getNumberOfDownCells(downId);
+ const int* downEdgeIds = grid->getDownArray(cellType)->getDownCells(downId);
+ const unsigned char* edgeType = grid->getDownArray(cellType)->getDownTypes(downId);
+ for (int ie =0; ie < nbEdges; ie++)
+ {
+ int nodes[3];
+ int nbNodes = grid->getDownArray(edgeType[ie])->getNodes(downEdgeIds[ie], nodes);
+ if ( mutipleNodes.count(nodes[0]) && mutipleNodes.count( nodes[ nbNodes-1 ]))
+ {
+ vector<int> vn0 = mutipleNodes[nodes[0]];
+ vector<int> vn1 = mutipleNodes[nodes[nbNodes - 1]];
+ vector<int> doms;
+ for ( size_t i0 = 0; i0 < vn0.size(); i0++ )
+ for ( size_t i1 = 0; i1 < vn1.size(); i1++ )
+ if ( vn0[i0] == vn1[i1] )
+ doms.push_back( vn0[ i0 ]);
+ if ( doms.size() > 2 )
+ {
+ //MESSAGE(" detect edgesMultiDomains " << nodes[0] << " " << nodes[nbNodes - 1]);
+ double *coords = grid->GetPoint(nodes[0]);
+ gp_Pnt p0(coords[0], coords[1], coords[2]);
+ coords = grid->GetPoint(nodes[nbNodes - 1]);
+ gp_Pnt p1(coords[0], coords[1], coords[2]);
+ gp_Pnt gref;
+ int vtkVolIds[1000]; // an edge can belong to a lot of volumes
+ map<int, SMDS_MeshVolume*> domvol; // domain --> a volume with the edge
+ map<int, double> angleDom; // oriented angles between planes defined by edge and volume centers
+ int nbvol = grid->GetParentVolumes(vtkVolIds, downEdgeIds[ie], edgeType[ie]);
+ for ( size_t id = 0; id < doms.size(); id++ )
+ {
+ int idom = doms[id];
+ const TIDSortedElemSet& domain = (idom == iRestDom) ? theRestDomElems : theElems[idom];
+ for ( int ivol = 0; ivol < nbvol; ivol++ )
+ {
+ int smdsId = meshDS->FromVtkToSmds(vtkVolIds[ivol]);
+ const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
+ if (domain.count(elem))
+ {
+ const SMDS_MeshVolume* svol = SMDS_Mesh::DownCast<SMDS_MeshVolume>(elem);
+ domvol[idom] = (SMDS_MeshVolume*) svol;
+ //MESSAGE(" domain " << idom << " volume " << elem->GetID());
+ double values[3] = { 0,0,0 };
+ vtkIdType npts = 0;
+ vtkIdType* pts = 0;
+ grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
+ for ( vtkIdType i = 0; i < npts; ++i )
+ {
+ double *coords = grid->GetPoint( pts[i] );
+ for ( int j = 0; j < 3; ++j )
+ values[j] += coords[j] / npts;
+ }
+ if ( id == 0 )
+ {
+ gref.SetCoord( values[0], values[1], values[2] );
+ angleDom[idom] = 0;
+ }
+ else
+ {
+ gp_Pnt g( values[0], values[1], values[2] );
+ angleDom[idom] = OrientedAngle(p0, p1, gref, g); // -pi<angle<+pi
+ //MESSAGE(" angle=" << angleDom[idom]);
+ }
+ break;
+ }
+ }
+ }
+ map<double, int> sortedDom; // sort domains by angle
+ for (map<int, double>::iterator ia = angleDom.begin(); ia != angleDom.end(); ++ia)
+ sortedDom[ia->second] = ia->first;
+ vector<int> vnodes;
+ vector<int> vdom;
+ for (map<double, int>::iterator ib = sortedDom.begin(); ib != sortedDom.end(); ++ib)
+ {
+ vdom.push_back(ib->second);
+ //MESSAGE(" ordered domain " << ib->second << " angle " << ib->first);
+ }
+ for (int ino = 0; ino < nbNodes; ino++)
+ vnodes.push_back(nodes[ino]);
+ edgesMultiDomains[vnodes] = vdom; // nodes vector --> ordered domains
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // --- iterate on shared faces (volumes to modify, face to extrude)
+ // get node id's of the face (id SMDS = id VTK)
+ // create flat element with old and new nodes if requested
+
+ // --- new quad nodes on flat quad elements: oldId --> ((domain1 X domain2) --> newId)
+ // (domain1 X domain2) = domain1 + MAXINT*domain2
+
+ std::map<int, std::map<long,int> > nodeQuadDomains;
+ std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
+
+ //MESSAGE(".. Creation of elements: simple junction");
+ if (createJointElems)
+ {
+ int idg;
+ string joints2DName = "joints2D";
+ mapOfJunctionGroups[joints2DName] = this->myMesh->AddGroup(SMDSAbs_Face, joints2DName.c_str(), idg);
+ SMESHDS_Group *joints2DGrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[joints2DName]->GetGroupDS());
+ string joints3DName = "joints3D";
+ mapOfJunctionGroups[joints3DName] = this->myMesh->AddGroup(SMDSAbs_Volume, joints3DName.c_str(), idg);
+ SMESHDS_Group *joints3DGrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[joints3DName]->GetGroupDS());
+
+ itface = faceDomains.begin();
+ for (; itface != faceDomains.end(); ++itface)
+ {
+ DownIdType face = itface->first;
+ std::set<int> oldNodes;
+ std::set<int>::iterator itn;
+ oldNodes.clear();
+ grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+
+ std::map<int, int> domvol = itface->second;
+ std::map<int, int>::iterator itdom = domvol.begin();
+ int dom1 = itdom->first;
+ int vtkVolId = itdom->second;
+ itdom++;
+ int dom2 = itdom->first;
+ SMDS_MeshCell *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
+ nodeQuadDomains);
+ stringstream grpname;
+ grpname << "j_";
+ if (dom1 < dom2)
+ grpname << dom1 << "_" << dom2;
+ else
+ grpname << dom2 << "_" << dom1;
+ string namegrp = grpname.str();
+ if (!mapOfJunctionGroups.count(namegrp))
+ mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(vol->GetType(), namegrp.c_str(), idg);
+ SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+ if (sgrp)
+ sgrp->Add(vol->GetID());
+ if (vol->GetType() == SMDSAbs_Volume)
+ joints3DGrp->Add(vol->GetID());
+ else if (vol->GetType() == SMDSAbs_Face)
+ joints2DGrp->Add(vol->GetID());
+ }
+ }
+
+ // --- create volumes on multiple domain intersection if requested
+ // iterate on mutipleNodesToFace
+ // iterate on edgesMultiDomains
+
+ //MESSAGE(".. Creation of elements: multiple junction");
+ if (createJointElems)
+ {
+ // --- iterate on mutipleNodesToFace
+
+ std::map<int, std::vector<int> >::iterator itn = mutipleNodesToFace.begin();
+ for (; itn != mutipleNodesToFace.end(); ++itn)
+ {
+ int node = itn->first;
+ vector<int> orderDom = itn->second;
+ vector<vtkIdType> orderedNodes;
+ for ( size_t idom = 0; idom < orderDom.size(); idom++ )
+ orderedNodes.push_back( nodeDomains[ node ][ orderDom[ idom ]]);
+ SMDS_MeshFace* face = this->GetMeshDS()->AddFaceFromVtkIds(orderedNodes);
+
+ stringstream grpname;
+ grpname << "m2j_";
+ grpname << 0 << "_" << 0;
+ int idg;
+ string namegrp = grpname.str();
+ if (!mapOfJunctionGroups.count(namegrp))
+ mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Face, namegrp.c_str(), idg);
+ SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+ if (sgrp)
+ sgrp->Add(face->GetID());
+ }
+
+ // --- iterate on edgesMultiDomains
+
+ std::map<std::vector<int>, std::vector<int> >::iterator ite = edgesMultiDomains.begin();
+ for (; ite != edgesMultiDomains.end(); ++ite)
+ {
+ vector<int> nodes = ite->first;
+ vector<int> orderDom = ite->second;
+ vector<vtkIdType> orderedNodes;
+ if (nodes.size() == 2)
+ {
+ //MESSAGE(" use edgesMultiDomains " << nodes[0] << " " << nodes[1]);
+ for ( size_t ino = 0; ino < nodes.size(); ino++ )
+ if ( orderDom.size() == 3 )
+ for ( size_t idom = 0; idom < orderDom.size(); idom++ )
+ orderedNodes.push_back( nodeDomains[ nodes[ ino ]][ orderDom[ idom ]]);
+ else
+ for (int idom = orderDom.size()-1; idom >=0; idom--)
+ orderedNodes.push_back( nodeDomains[ nodes[ ino ]][ orderDom[ idom ]]);
+ SMDS_MeshVolume* vol = this->GetMeshDS()->AddVolumeFromVtkIds(orderedNodes);
+
+ int idg;
+ string namegrp = "jointsMultiples";
+ if (!mapOfJunctionGroups.count(namegrp))
+ mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+ SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+ if (sgrp)
+ sgrp->Add(vol->GetID());
+ }
+ else
+ {
+ //INFOS("Quadratic multiple joints not implemented");
+ // TODO quadratic nodes
+ }
+ }
+ }
+
+ // --- list the explicit faces and edges of the mesh that need to be modified,
+ // i.e. faces and edges built with one or more duplicated nodes.
+ // associate these faces or edges to their corresponding domain.
+ // only the first domain found is kept when a face or edge is shared
+
+ std::map<DownIdType, std::map<int,int>, DownIdCompare> faceOrEdgeDom; // cellToModify --> (id domain --> id cell)
+ std::map<int,int> feDom; // vtk id of cell to modify --> id domain
+ faceOrEdgeDom.clear();
+ feDom.clear();
+
+ //MESSAGE(".. Modification of elements");
+ for (int idomain = idom0; idomain < nbDomains; idomain++)
+ {
+ std::map<int, std::map<int, int> >::const_iterator itnod = nodeDomains.begin();
+ for (; itnod != nodeDomains.end(); ++itnod)
+ {
+ int oldId = itnod->first;
+ //MESSAGE(" node " << oldId);
+ vtkCellLinks::Link l = grid->GetCellLinks()->GetLink(oldId);
+ for (int i = 0; i < l.ncells; i++)
+ {
+ int vtkId = l.cells[i];
+ int vtkType = grid->GetCellType(vtkId);
+ int downId = grid->CellIdToDownId(vtkId);
+ if (downId < 0)
+ continue; // new cells: not to be modified
+ DownIdType aCell(downId, vtkType);
+ int volParents[1000];
+ int nbvol = grid->GetParentVolumes(volParents, vtkId);
+ for (int j = 0; j < nbvol; j++)
+ if (celldom.count(volParents[j]) && (celldom[volParents[j]] == idomain))
+ if (!feDom.count(vtkId))
+ {
+ feDom[vtkId] = idomain;
+ faceOrEdgeDom[aCell] = emptyMap;
+ faceOrEdgeDom[aCell][idomain] = vtkId; // affect face or edge to the first domain only
+ //MESSAGE("affect cell " << this->GetMeshDS()->FromVtkToSmds(vtkId) << " domain " << idomain
+ // << " type " << vtkType << " downId " << downId);
+ }
+ }
+ }
+ }
+
+ // --- iterate on shared faces (volumes to modify, face to extrude)
+ // get node id's of the face
+ // replace old nodes by new nodes in volumes, and update inverse connectivity
+
+ std::map<DownIdType, std::map<int,int>, DownIdCompare>* maps[3] = {&faceDomains, &cellDomains, &faceOrEdgeDom};
+ for (int m=0; m<3; m++)
+ {
+ std::map<DownIdType, std::map<int,int>, DownIdCompare>* amap = maps[m];
+ itface = (*amap).begin();
+ for (; itface != (*amap).end(); ++itface)
+ {
+ DownIdType face = itface->first;
+ std::set<int> oldNodes;
+ std::set<int>::iterator itn;
+ oldNodes.clear();
+ grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+ //MESSAGE("examine cell, downId " << face.cellId << " type " << int(face.cellType));
+ std::map<int, int> localClonedNodeIds;
+
+ std::map<int, int> domvol = itface->second;
+ std::map<int, int>::iterator itdom = domvol.begin();
+ for (; itdom != domvol.end(); ++itdom)
+ {
+ int idom = itdom->first;
+ int vtkVolId = itdom->second;
+ //MESSAGE("modify nodes of cell " << this->GetMeshDS()->FromVtkToSmds(vtkVolId) << " domain " << idom);
+ localClonedNodeIds.clear();
+ for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn)
+ {
+ int oldId = *itn;
+ if (nodeDomains[oldId].count(idom))
+ {
+ localClonedNodeIds[oldId] = nodeDomains[oldId][idom];
+ //MESSAGE(" node " << oldId << " --> " << localClonedNodeIds[oldId]);
+ }
+ }
+ meshDS->ModifyCellNodes(vtkVolId, localClonedNodeIds);
+ }
+ }
+ }
+
+ // Remove empty groups (issue 0022812)
+ std::map<std::string, SMESH_Group*>::iterator name_group = mapOfJunctionGroups.begin();
+ for ( ; name_group != mapOfJunctionGroups.end(); ++name_group )
+ {
+ if ( name_group->second && name_group->second->GetGroupDS()->IsEmpty() )
+ myMesh->RemoveGroup( name_group->second->GetGroupDS()->GetID() );
+ }
+
+ meshDS->CleanDownWardConnectivity(); // Mesh has been modified, downward connectivity is no more usable, free memory
+ grid->DeleteLinks();
+
+ CHRONOSTOP(50);
+ counters::stats();
+ return true;
+}
+
+/*!
+ * \brief Double nodes on some external faces and create flat elements.
+ * Flat elements are mainly used by some types of mechanic calculations.
+ *
+ * Each group of the list must be constituted of faces.
+ * Triangles are transformed in prisms, and quadrangles in hexahedrons.
+ * @param theElems - list of groups of faces, where a group of faces is a set of
+ * SMDS_MeshElements sorted by Id.
+ * @return TRUE if operation has been completed successfully, FALSE otherwise
+ */
+bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSortedElemSet>& theElems)
+{
+ // MESSAGE("-------------------------------------------------");
+ // MESSAGE("SMESH_MeshEditor::CreateFlatElementsOnFacesGroups");
+ // MESSAGE("-------------------------------------------------");
+
+ SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
+
+ // --- For each group of faces
+ // duplicate the nodes, create a flat element based on the face
+ // replace the nodes of the faces by their clones
+
+ std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> clonedNodes;
+ std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> intermediateNodes;
+ clonedNodes.clear();
+ intermediateNodes.clear();
+ std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
+ mapOfJunctionGroups.clear();
+
+ for ( size_t idom = 0; idom < theElems.size(); idom++ )
+ {
+ const TIDSortedElemSet& domain = theElems[idom];
+ TIDSortedElemSet::const_iterator elemItr = domain.begin();
+ for ( ; elemItr != domain.end(); ++elemItr )
+ {
+ const SMDS_MeshFace* aFace = meshDS->DownCast<SMDS_MeshFace> ( *elemItr );
+ if (!aFace)
+ continue;
+ // MESSAGE("aFace=" << aFace->GetID());
+ bool isQuad = aFace->IsQuadratic();
+ vector<const SMDS_MeshNode*> ln0, ln1, ln2, ln3, ln4;
+
+ // --- clone the nodes, create intermediate nodes for non medium nodes of a quad face
+
+ SMDS_NodeIteratorPtr nodeIt = aFace->nodeIterator();
+ while (nodeIt->more())
+ {
+ const SMDS_MeshNode* node = nodeIt->next();
+ bool isMedium = ( isQuad && aFace->IsMediumNode( node ));
+ if (isMedium)
+ ln2.push_back(node);
+ else
+ ln0.push_back(node);
+
+ const SMDS_MeshNode* clone = 0;
+ if (!clonedNodes.count(node))
+ {
+ clone = meshDS->AddNode(node->X(), node->Y(), node->Z());
+ copyPosition( node, clone );
+ clonedNodes[node] = clone;
+ }
+ else
+ clone = clonedNodes[node];
+
+ if (isMedium)
+ ln3.push_back(clone);
+ else
+ ln1.push_back(clone);
+
+ const SMDS_MeshNode* inter = 0;
+ if (isQuad && (!isMedium))
+ {
+ if (!intermediateNodes.count(node))
+ {
+ inter = meshDS->AddNode(node->X(), node->Y(), node->Z());
+ copyPosition( node, inter );
+ intermediateNodes[node] = inter;
+ }
+ else
+ inter = intermediateNodes[node];
+ ln4.push_back(inter);
+ }
+ }
+
+ // --- extrude the face
+
+ vector<const SMDS_MeshNode*> ln;
+ SMDS_MeshVolume* vol = 0;
+ vtkIdType aType = aFace->GetVtkType();
+ switch (aType)
+ {
+ case VTK_TRIANGLE:
+ vol = meshDS->AddVolume(ln0[2], ln0[1], ln0[0], ln1[2], ln1[1], ln1[0]);
+ // MESSAGE("vol prism " << vol->GetID());
+ ln.push_back(ln1[0]);
+ ln.push_back(ln1[1]);
+ ln.push_back(ln1[2]);
+ break;
+ case VTK_QUAD:
+ vol = meshDS->AddVolume(ln0[3], ln0[2], ln0[1], ln0[0], ln1[3], ln1[2], ln1[1], ln1[0]);
+ // MESSAGE("vol hexa " << vol->GetID());
+ ln.push_back(ln1[0]);
+ ln.push_back(ln1[1]);
+ ln.push_back(ln1[2]);
+ ln.push_back(ln1[3]);
+ break;
+ case VTK_QUADRATIC_TRIANGLE:
+ vol = meshDS->AddVolume(ln1[0], ln1[1], ln1[2], ln0[0], ln0[1], ln0[2], ln3[0], ln3[1], ln3[2],
+ ln2[0], ln2[1], ln2[2], ln4[0], ln4[1], ln4[2]);
+ // MESSAGE("vol quad prism " << vol->GetID());
+ ln.push_back(ln1[0]);
+ ln.push_back(ln1[1]);
+ ln.push_back(ln1[2]);
+ ln.push_back(ln3[0]);
+ ln.push_back(ln3[1]);
+ ln.push_back(ln3[2]);
+ break;
+ case VTK_QUADRATIC_QUAD:
+ // vol = meshDS->AddVolume(ln0[0], ln0[1], ln0[2], ln0[3], ln1[0], ln1[1], ln1[2], ln1[3],
+ // ln2[0], ln2[1], ln2[2], ln2[3], ln3[0], ln3[1], ln3[2], ln3[3],
+ // ln4[0], ln4[1], ln4[2], ln4[3]);
+ vol = meshDS->AddVolume(ln1[0], ln1[1], ln1[2], ln1[3], ln0[0], ln0[1], ln0[2], ln0[3],
+ ln3[0], ln3[1], ln3[2], ln3[3], ln2[0], ln2[1], ln2[2], ln2[3],
+ ln4[0], ln4[1], ln4[2], ln4[3]);
+ // MESSAGE("vol quad hexa " << vol->GetID());
+ ln.push_back(ln1[0]);
+ ln.push_back(ln1[1]);
+ ln.push_back(ln1[2]);
+ ln.push_back(ln1[3]);
+ ln.push_back(ln3[0]);
+ ln.push_back(ln3[1]);
+ ln.push_back(ln3[2]);
+ ln.push_back(ln3[3]);
+ break;
+ case VTK_POLYGON:
+ break;
+ default:
+ break;
+ }
+
+ if (vol)
+ {
+ stringstream grpname;
+ grpname << "jf_";
+ grpname << idom;
+ int idg;
+ string namegrp = grpname.str();
+ if (!mapOfJunctionGroups.count(namegrp))
+ mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+ SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+ if (sgrp)
+ sgrp->Add(vol->GetID());
+ }
+
+ // --- modify the face
+
+ const_cast<SMDS_MeshFace*>( aFace )->ChangeNodes( &ln[0], ln.size() );
+ }
+ }
+ return true;
+}
+
+/*!
+ * \brief identify all the elements around a geom shape, get the faces delimiting the hole
+ * Build groups of volume to remove, groups of faces to replace on the skin of the object,
+ * groups of faces to remove inside the object, (idem edges).
+ * Build ordered list of nodes at the border of each group of faces to replace (to be used to build a geom subshape)
+ */
+void SMESH_MeshEditor::CreateHoleSkin(double radius,
+ const TopoDS_Shape& theShape,
+ SMESH_NodeSearcher* theNodeSearcher,
+ const char* groupName,
+ std::vector<double>& nodesCoords,
+ std::vector<std::vector<int> >& listOfListOfNodes)
+{
+ // MESSAGE("--------------------------------");
+ // MESSAGE("SMESH_MeshEditor::CreateHoleSkin");
+ // MESSAGE("--------------------------------");
+
+ // --- zone of volumes to remove is given :
+ // 1 either by a geom shape (one or more vertices) and a radius,
+ // 2 either by a group of nodes (representative of the shape)to use with the radius,
+ // 3 either by a group of nodes where all the elements build on one of this nodes are to remove,
+ // In the case 2, the group of nodes is an external group of nodes from another mesh,
+ // In the case 3, the group of nodes is an internal group of the mesh (obtained for instance by a filter),
+ // defined by it's name.
+
+ SMESHDS_GroupBase* groupDS = 0;
+ SMESH_Mesh::GroupIteratorPtr groupIt = this->myMesh->GetGroups();
+ while ( groupIt->more() )
+ {
+ groupDS = 0;
+ SMESH_Group * group = groupIt->next();
+ if ( !group ) continue;
+ groupDS = group->GetGroupDS();
+ if ( !groupDS || groupDS->IsEmpty() ) continue;
+ std::string grpName = group->GetName();
+ //MESSAGE("grpName=" << grpName);
+ if (grpName == groupName)
+ break;
+ else
+ groupDS = 0;
+ }
+
+ bool isNodeGroup = false;
+ bool isNodeCoords = false;
+ if (groupDS)
+ {
+ if (groupDS->GetType() != SMDSAbs_Node)
+ return;
+ isNodeGroup = true; // a group of nodes exists and it is in this mesh
+ }
+
+ if (nodesCoords.size() > 0)
+ isNodeCoords = true; // a list o nodes given by their coordinates
+ //MESSAGE("---" << isNodeGroup << " " << isNodeCoords);
+
+ // --- define groups to build
+
+ int idg; // --- group of SMDS volumes
+ string grpvName = groupName;
+ grpvName += "_vol";
+ SMESH_Group *grp = this->myMesh->AddGroup(SMDSAbs_Volume, grpvName.c_str(), idg);
+ if (!grp)
+ {
+ MESSAGE("group not created " << grpvName);
+ return;
+ }
+ SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(grp->GetGroupDS());
+
+ int idgs; // --- group of SMDS faces on the skin
+ string grpsName = groupName;
+ grpsName += "_skin";
+ SMESH_Group *grps = this->myMesh->AddGroup(SMDSAbs_Face, grpsName.c_str(), idgs);
+ if (!grps)
+ {
+ MESSAGE("group not created " << grpsName);
+ return;
+ }
+ SMESHDS_Group *sgrps = dynamic_cast<SMESHDS_Group*>(grps->GetGroupDS());
+
+ int idgi; // --- group of SMDS faces internal (several shapes)
+ string grpiName = groupName;
+ grpiName += "_internalFaces";
+ SMESH_Group *grpi = this->myMesh->AddGroup(SMDSAbs_Face, grpiName.c_str(), idgi);
+ if (!grpi)
+ {
+ MESSAGE("group not created " << grpiName);
+ return;
+ }
+ SMESHDS_Group *sgrpi = dynamic_cast<SMESHDS_Group*>(grpi->GetGroupDS());
+
+ int idgei; // --- group of SMDS faces internal (several shapes)
+ string grpeiName = groupName;
+ grpeiName += "_internalEdges";
+ SMESH_Group *grpei = this->myMesh->AddGroup(SMDSAbs_Edge, grpeiName.c_str(), idgei);
+ if (!grpei)
+ {
+ MESSAGE("group not created " << grpeiName);
+ return;
+ }
+ SMESHDS_Group *sgrpei = dynamic_cast<SMESHDS_Group*>(grpei->GetGroupDS());
+
+ // --- build downward connectivity
+
+ SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
+ meshDS->BuildDownWardConnectivity(true);
+ SMDS_UnstructuredGrid* grid = meshDS->GetGrid();
+
+ // --- set of volumes detected inside
+
+ std::set<int> setOfInsideVol;
+ std::set<int> setOfVolToCheck;
+
+ std::vector<gp_Pnt> gpnts;
+ gpnts.clear();
+
+ if (isNodeGroup) // --- a group of nodes is provided : find all the volumes using one or more of this nodes
+ {
+ //MESSAGE("group of nodes provided");
+ SMDS_ElemIteratorPtr elemIt = groupDS->GetElements();
+ while ( elemIt->more() )
+ {
+ const SMDS_MeshElement* elem = elemIt->next();
+ if (!elem)
+ continue;
+ const SMDS_MeshNode* node = dynamic_cast<const SMDS_MeshNode*>(elem);
+ if (!node)
+ continue;
+ SMDS_MeshElement* vol = 0;
+ SMDS_ElemIteratorPtr volItr = node->GetInverseElementIterator(SMDSAbs_Volume);
+ while (volItr->more())
+ {
+ vol = (SMDS_MeshElement*)volItr->next();
+ setOfInsideVol.insert(vol->GetVtkID());
+ sgrp->Add(vol->GetID());
+ }
+ }
+ }
+ else if (isNodeCoords)
+ {
+ //MESSAGE("list of nodes coordinates provided");
+ size_t i = 0;
+ int k = 0;
+ while ( i < nodesCoords.size()-2 )
+ {
+ double x = nodesCoords[i++];
+ double y = nodesCoords[i++];
+ double z = nodesCoords[i++];
+ gp_Pnt p = gp_Pnt(x, y ,z);
+ gpnts.push_back(p);
+ //MESSAGE("TopoDS_Vertex " << k << " " << p.X() << " " << p.Y() << " " << p.Z());
+ k++;
+ }
+ }
+ else // --- no group, no coordinates : use the vertices of the geom shape provided, and radius
+ {
+ //MESSAGE("no group of nodes provided, using vertices from geom shape, and radius");
+ TopTools_IndexedMapOfShape vertexMap;
+ TopExp::MapShapes( theShape, TopAbs_VERTEX, vertexMap );
+ gp_Pnt p = gp_Pnt(0,0,0);
+ if (vertexMap.Extent() < 1)
+ return;
+
+ for ( int i = 1; i <= vertexMap.Extent(); ++i )
+ {
+ const TopoDS_Vertex& vertex = TopoDS::Vertex( vertexMap( i ));
+ p = BRep_Tool::Pnt(vertex);
+ gpnts.push_back(p);
+ //MESSAGE("TopoDS_Vertex " << i << " " << p.X() << " " << p.Y() << " " << p.Z());
+ }
+ }
+
+ if (gpnts.size() > 0)
+ {
+ const SMDS_MeshNode* startNode = theNodeSearcher->FindClosestTo(gpnts[0]);
+ //MESSAGE("startNode->nodeId " << nodeId);
+
+ double radius2 = radius*radius;
+ //MESSAGE("radius2 " << radius2);
+
+ // --- volumes on start node
+
+ setOfVolToCheck.clear();
+ SMDS_MeshElement* startVol = 0;
+ SMDS_ElemIteratorPtr volItr = startNode->GetInverseElementIterator(SMDSAbs_Volume);
+ while (volItr->more())
+ {
+ startVol = (SMDS_MeshElement*)volItr->next();
+ setOfVolToCheck.insert(startVol->GetVtkID());
+ }
+ if (setOfVolToCheck.empty())
+ {
+ MESSAGE("No volumes found");
+ return;
+ }
+
+ // --- starting with central volumes then their neighbors, check if they are inside
+ // or outside the domain, until no more new neighbor volume is inside.
+ // Fill the group of inside volumes
+
+ std::map<int, double> mapOfNodeDistance2;
+ mapOfNodeDistance2.clear();
+ std::set<int> setOfOutsideVol;
+ while (!setOfVolToCheck.empty())
+ {
+ std::set<int>::iterator it = setOfVolToCheck.begin();
+ int vtkId = *it;
+ //MESSAGE("volume to check, vtkId " << vtkId << " smdsId " << meshDS->FromVtkToSmds(vtkId));
+ bool volInside = false;
+ vtkIdType npts = 0;
+ vtkIdType* pts = 0;
+ grid->GetCellPoints(vtkId, npts, pts);
+ for (int i=0; i<npts; i++)
+ {
+ double distance2 = 0;
+ if (mapOfNodeDistance2.count(pts[i]))
+ {
+ distance2 = mapOfNodeDistance2[pts[i]];
+ //MESSAGE("point " << pts[i] << " distance2 " << distance2);
+ }
+ else
+ {
+ double *coords = grid->GetPoint(pts[i]);
+ gp_Pnt aPoint = gp_Pnt(coords[0], coords[1], coords[2]);
+ distance2 = 1.E40;
+ for ( size_t j = 0; j < gpnts.size(); j++ )
+ {
+ double d2 = aPoint.SquareDistance( gpnts[ j ]);
+ if (d2 < distance2)
+ {
+ distance2 = d2;
+ if (distance2 < radius2)
+ break;
+ }
+ }
+ mapOfNodeDistance2[pts[i]] = distance2;
+ //MESSAGE(" point " << pts[i] << " distance2 " << distance2 << " coords " << coords[0] << " " << coords[1] << " " << coords[2]);
+ }
+ if (distance2 < radius2)
+ {
+ volInside = true; // one or more nodes inside the domain
+ sgrp->Add(meshDS->FromVtkToSmds(vtkId));
+ break;
+ }
+ }
+ if (volInside)
+ {
+ setOfInsideVol.insert(vtkId);
+ //MESSAGE(" volume inside, vtkId " << vtkId << " smdsId " << meshDS->FromVtkToSmds(vtkId));
+ int neighborsVtkIds[NBMAXNEIGHBORS];
+ int downIds[NBMAXNEIGHBORS];
+ unsigned char downTypes[NBMAXNEIGHBORS];
+ int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId);
+ for (int n = 0; n < nbNeighbors; n++)
+ if (!setOfInsideVol.count(neighborsVtkIds[n]) ||setOfOutsideVol.count(neighborsVtkIds[n]))
+ setOfVolToCheck.insert(neighborsVtkIds[n]);