+
+//================================================================================
+/*!
+ * \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;
+}
+
+//================================================================================
+/*!
+ \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 )
+{
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ if ( theElems.size() == 0 )
+ return false;
+
+ SMESHDS_Mesh* aMeshDS = GetMeshDS();
+ if ( !aMeshDS )
+ return false;
+
+ bool res = false;
+ std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > 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,
+ std::map< const SMDS_MeshNode*,
+ const SMDS_MeshNode* >& theNodeNodeMap,
+ const bool theIsDoubleElem )
+{
+ MESSAGE("doubleNodes");
+ // iterate on through element and duplicate them (by nodes duplication)
+ bool res = false;
+ TIDSortedElemSet::const_iterator elemItr = theElems.begin();
+ for ( ; elemItr != theElems.end(); ++elemItr )
+ {
+ const SMDS_MeshElement* anElem = *elemItr;
+ if (!anElem)
+ continue;
+
+ bool isDuplicate = false;
+ // duplicate nodes to duplicate element
+ std::vector<const SMDS_MeshNode*> newNodes( anElem->NbNodes() );
+ SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
+ int ind = 0;
+ while ( anIter->more() )
+ {
+
+ SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
+ SMDS_MeshNode* aNewNode = aCurrNode;
+ if ( theNodeNodeMap.find( aCurrNode ) != theNodeNodeMap.end() )
+ aNewNode = (SMDS_MeshNode*)theNodeNodeMap[ aCurrNode ];
+ else if ( theIsDoubleElem && theNodesNot.find( aCurrNode ) == theNodesNot.end() )
+ {
+ // duplicate node
+ aNewNode = theMeshDS->AddNode( aCurrNode->X(), aCurrNode->Y(), aCurrNode->Z() );
+ theNodeNodeMap[ aCurrNode ] = aNewNode;
+ myLastCreatedNodes.Append( aNewNode );
+ }
+ isDuplicate |= (aCurrNode != aNewNode);
+ newNodes[ ind++ ] = aNewNode;
+ }
+ if ( !isDuplicate )
+ continue;
+
+ if ( theIsDoubleElem )
+ AddElement(newNodes, anElem->GetType(), anElem->IsPoly());
+ else
+ {
+ MESSAGE("ChangeElementNodes");
+ theMeshDS->ChangeElementNodes( anElem, &newNodes[ 0 ], anElem->NbNodes() );
+ }
+ 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 )
+{
+ MESSAGE("DoubleNodes");
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ 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 )
+ {
+ int aCurr = *aNodeIter;
+ SMDS_MeshNode* aNode = (SMDS_MeshNode*)aMeshDS->FindNode( aCurr );
+ if ( !aNode )
+ continue;
+
+ // duplicate node
+
+ const SMDS_MeshNode* aNewNode = aMeshDS->AddNode( aNode->X(), aNode->Y(), aNode->Z() );
+ if ( aNewNode )
+ {
+ anOldNodeToNewNode[ aNode ] = aNewNode;
+ myLastCreatedNodes.Append( aNewNode );
+ }
+ }
+
+ // Create map of new nodes for modified elements
+
+ std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> > anElemToNodes;
+
+ std::list< int >::const_iterator anElemIter;
+ for ( anElemIter = theListOfModifiedElems.begin();
+ anElemIter != theListOfModifiedElems.end(); ++anElemIter )
+ {
+ int aCurr = *anElemIter;
+ SMDS_MeshElement* anElem = (SMDS_MeshElement*)aMeshDS->FindElement( aCurr );
+ if ( !anElem )
+ continue;
+
+ vector<const SMDS_MeshNode*> aNodeArr( anElem->NbNodes() );
+
+ SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
+ int ind = 0;
+ while ( anIter->more() )
+ {
+ SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
+ if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() )
+ {
+ const SMDS_MeshNode* aNewNode = anOldNodeToNewNode[ aCurrNode ];
+ aNodeArr[ ind++ ] = aNewNode;
+ }
+ else
+ aNodeArr[ ind++ ] = aCurrNode;
+ }
+ anElemToNodes[ anElem ] = aNodeArr;
+ }
+
+ // Change nodes of elements
+
+ std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> >::iterator
+ anElemToNodesIter = anElemToNodes.begin();
+ for ( ; anElemToNodesIter != anElemToNodes.end(); ++anElemToNodesIter )
+ {
+ const SMDS_MeshElement* anElem = anElemToNodesIter->first;
+ vector<const SMDS_MeshNode*> aNodeArr = anElemToNodesIter->second;
+ if ( anElem )
+ {
+ MESSAGE("ChangeElementNodes");
+ aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], anElem->NbNodes() );
+ }
+ }
+
+ 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);
+ SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator();
+ while (aNodeItr->more())
+ centerXYZ += SMESH_TNodeXYZ(cast2Node( 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)
+ {
+ _state = TopAbs_OUT;
+ _extremum.Perform(aPnt);
+ if ( _extremum.IsDone() )
+ for ( int iSol = 1; iSol <= _extremum.NbExt() && _state == TopAbs_OUT; ++iSol)
+#if OCC_VERSION_LARGE > 0x06040000 // Porting to OCCT6.5.1
+ _state = ( _extremum.SquareDistance(iSol) <= theTol ? TopAbs_IN : TopAbs_OUT );
+#else
+ _state = ( _extremum.Value(iSol) <= theTol ? TopAbs_IN : TopAbs_OUT );
+#endif
+ }
+ 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).
+ The replicated nodes should be associated to affected elements.
+ \return groups of affected elements
+ \sa DoubleNodeElemGroupsInRegion()
+ */
+//================================================================================
+
+bool SMESH_MeshEditor::AffectedElemGroupsInRegion( const TIDSortedElemSet& theElems,
+ const TIDSortedElemSet& theNodesNot,
+ const TopoDS_Shape& theShape,
+ TIDSortedElemSet& theAffectedElems)
+{
+ if ( theShape.IsNull() )
+ return false;
+
+ const double aTol = Precision::Confusion();
+ auto_ptr< BRepClass3d_SolidClassifier> bsc3d;
+ auto_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;
+ 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.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();
+ auto_ptr< BRepClass3d_SolidClassifier> bsc3d;
+ auto_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 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.get() ?
+ 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)
+{
+// MESSAGE(" p0: " << p0.X() << " " << p0.Y() << " " << p0.Z());
+// MESSAGE(" p1: " << p1.X() << " " << p1.Y() << " " << p1.Z());
+// MESSAGE(" g1: " << g1.X() << " " << g1.Y() << " " << g1.Z());
+// MESSAGE(" g2: " << g2.X() << " " << g2.Y() << " " << g2.Z());
+ 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);
+ return n2.AngleWithRef(n1, vref);
+}
+
+/*!
+ * \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
+ * @return TRUE if operation has been completed successfully, FALSE otherwise
+ */
+bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSortedElemSet>& theElems,
+ bool createJointElems)
+{
+ 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());
+
+ // Check if the domains do not share an element
+ for (int idom = 0; idom < theElems.size()-1; idom++)
+ {
+// MESSAGE("... Check of domain #" << idom);
+ const TIDSortedElemSet& domain = theElems[idom];
+ TIDSortedElemSet::const_iterator elemItr = domain.begin();
+ for (; elemItr != domain.end(); ++elemItr)
+ {
+ SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr;
+ int idombisdeb = idom + 1 ;
+ for (int idombis = idombisdeb; idombis < theElems.size(); idombis++) // check if the element belongs to a domain further in the list
+ {
+ 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 < theElems.size(); 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)
+ {
+ SMDS_MeshElement* anElem = (SMDS_MeshElement*) *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 (! domain.count(elem)) // neighbor is in another domain : face is shared
+ {
+ bool ok = false ;
+ for (int idombis = 0; idombis < theElems.size(); 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 ) // the characteristics of the face is stored
+ {
+ DownIdType face(downIds[n], downTypes[n]);
+ if (!faceDomains.count(face))
+ faceDomains[face] = emptyMap; // create an empty entry for face
+ 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);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ //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 = 0; idomain < theElems.size(); idomain++)
+ {
+ //MESSAGE("Domain " << idomain);
+ const TIDSortedElemSet& domain = theElems[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);
+ 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);
+ if (!cellDomains.count(aCell))
+ cellDomains[aCell] = emptyMap; // create an empty entry for cell
+ 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 = 0; idomain < theElems.size(); 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);
+ std::set<int>::iterator itn = oldNodes.begin();
+ for (; itn != oldNodes.end(); ++itn)
+ {
+ int oldId = *itn;
+ //MESSAGE("-+-+-a node " << oldId);
+ if (!nodeDomains.count(oldId))
+ nodeDomains[oldId] = emptyMap; // create an empty entry for node
+ 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>::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]);
+ 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 = 0; idomain < theElems.size(); 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 (int i0 = 0; i0 < vn0.size(); i0++)
+ for (int 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_VtkVolume*> 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 (int id=0; id < doms.size(); id++)
+ {
+ int idom = doms[id];
+ for (int ivol=0; ivol<nbvol; ivol++)
+ {
+ int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
+ SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
+ if (theElems[idom].count(elem))
+ {
+ SMDS_VtkVolume* svol = dynamic_cast<SMDS_VtkVolume*>(elem);
+ domvol[idom] = svol;
+ //MESSAGE(" domain " << idom << " volume " << elem->GetID());
+ double values[3];
+ vtkIdType npts = 0;
+ vtkIdType* pts = 0;
+ grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
+ SMDS_VtkVolume::gravityCenter(grid, pts, npts, values);
+ if (id ==0)
+ {
+ gref.SetXYZ(gp_XYZ(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 (int 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 (int ino=0; ino < nodes.size(); ino++)
+ if (orderDom.size() == 3)
+ for (int 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 = 0; idomain < theElems.size(); 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);
+ }
+ }
+ }
+
+ meshDS->CleanDownWardConnectivity(); // Mesh has been modified, downward connectivity is no more usable, free memory
+ grid->BuildLinks();
+
+ 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 (int idom = 0; idom < theElems.size(); idom++)
+ {
+ const TIDSortedElemSet& domain = theElems[idom];
+ TIDSortedElemSet::const_iterator elemItr = domain.begin();
+ for (; elemItr != domain.end(); ++elemItr)
+ {
+ SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr;
+ SMDS_MeshFace* aFace = dynamic_cast<SMDS_MeshFace*> (anElem);
+ 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_ElemIteratorPtr nodeIt = aFace->nodesIterator();
+ while (nodeIt->more())
+ {
+ const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*> (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());
+ 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());
+ 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
+
+ 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");
+ int 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());
+ }
+ }
+ 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)
+ {
+ int nodeId = 0;
+ const SMDS_MeshNode* startNode = theNodeSearcher->FindClosestTo(gpnts[0]);
+ if (startNode)
+ nodeId = startNode->GetID();
+ MESSAGE("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 (int 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]);
+ }
+ else
+ {
+ setOfOutsideVol.insert(vtkId);
+ MESSAGE(" volume outside, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+ }
+ setOfVolToCheck.erase(vtkId);
+ }
+ }
+
+ // --- for outside hexahedrons, check if they have more than one neighbor volume inside
+ // If yes, add the volume to the inside set
+
+ bool addedInside = true;
+ std::set<int> setOfVolToReCheck;
+ while (addedInside)
+ {
+ MESSAGE(" --------------------------- re check");
+ addedInside = false;
+ std::set<int>::iterator itv = setOfInsideVol.begin();
+ for (; itv != setOfInsideVol.end(); ++itv)
+ {
+ int vtkId = *itv;
+ 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]))
+ setOfVolToReCheck.insert(neighborsVtkIds[n]);
+ }
+ setOfVolToCheck = setOfVolToReCheck;
+ setOfVolToReCheck.clear();
+ while (!setOfVolToCheck.empty())
+ {
+ std::set<int>::iterator it = setOfVolToCheck.begin();
+ int vtkId = *it;
+ if (grid->GetCellType(vtkId) == VTK_HEXAHEDRON)
+ {
+ MESSAGE("volume to recheck, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+ int countInside = 0;
+ 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]))
+ countInside++;
+ MESSAGE("countInside " << countInside);
+ if (countInside > 1)
+ {
+ MESSAGE(" volume inside, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+ setOfInsideVol.insert(vtkId);
+ sgrp->Add(meshDS->fromVtkToSmds(vtkId));
+ addedInside = true;
+ }
+ else
+ setOfVolToReCheck.insert(vtkId);
+ }
+ setOfVolToCheck.erase(vtkId);
+ }
+ }
+
+ // --- map of Downward faces at the boundary, inside the global volume
+ // map of Downward faces on the skin of the global volume (equivalent to SMDS faces on the skin)
+ // fill group of SMDS faces inside the volume (when several volume shapes)
+ // fill group of SMDS faces on the skin of the global volume (if skin)
+
+ std::map<DownIdType, int, DownIdCompare> boundaryFaces; // boundary faces inside the volume --> corresponding cell
+ std::map<DownIdType, int, DownIdCompare> skinFaces; // faces on the skin of the global volume --> corresponding cell
+ std::set<int>::iterator it = setOfInsideVol.begin();
+ for (; it != setOfInsideVol.end(); ++it)
+ {
+ int vtkId = *it;
+ //MESSAGE(" vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+ int neighborsVtkIds[NBMAXNEIGHBORS];
+ int downIds[NBMAXNEIGHBORS];
+ unsigned char downTypes[NBMAXNEIGHBORS];
+ int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId, true);
+ for (int n = 0; n < nbNeighbors; n++)
+ {
+ int neighborDim = SMDS_Downward::getCellDimension(grid->GetCellType(neighborsVtkIds[n]));
+ if (neighborDim == 3)
+ {
+ if (! setOfInsideVol.count(neighborsVtkIds[n])) // neighbor volume is not inside : face is boundary
+ {
+ DownIdType face(downIds[n], downTypes[n]);
+ boundaryFaces[face] = vtkId;
+ }
+ // if the face between to volumes is in the mesh, get it (internal face between shapes)
+ int vtkFaceId = grid->getDownArray(downTypes[n])->getVtkCellId(downIds[n]);
+ if (vtkFaceId >= 0)
+ {
+ sgrpi->Add(meshDS->fromVtkToSmds(vtkFaceId));
+ // find also the smds edges on this face
+ int nbEdges = grid->getDownArray(downTypes[n])->getNumberOfDownCells(downIds[n]);
+ const int* dEdges = grid->getDownArray(downTypes[n])->getDownCells(downIds[n]);
+ const unsigned char* dTypes = grid->getDownArray(downTypes[n])->getDownTypes(downIds[n]);
+ for (int i = 0; i < nbEdges; i++)
+ {
+ int vtkEdgeId = grid->getDownArray(dTypes[i])->getVtkCellId(dEdges[i]);
+ if (vtkEdgeId >= 0)
+ sgrpei->Add(meshDS->fromVtkToSmds(vtkEdgeId));
+ }
+ }
+ }
+ else if (neighborDim == 2) // skin of the volume
+ {
+ DownIdType face(downIds[n], downTypes[n]);
+ skinFaces[face] = vtkId;
+ int vtkFaceId = grid->getDownArray(downTypes[n])->getVtkCellId(downIds[n]);
+ if (vtkFaceId >= 0)
+ sgrps->Add(meshDS->fromVtkToSmds(vtkFaceId));
+ }
+ }
+ }
+
+ // --- identify the edges constituting the wire of each subshape on the skin
+ // define polylines with the nodes of edges, equivalent to wires
+ // project polylines on subshapes, and partition, to get geom faces
+
+ std::map<int, std::set<int> > shapeIdToVtkIdSet; // shapeId --> set of vtkId on skin
+ std::set<int> emptySet;
+ emptySet.clear();
+ std::set<int> shapeIds;
+
+ SMDS_ElemIteratorPtr itelem = sgrps->GetElements();
+ while (itelem->more())
+ {
+ const SMDS_MeshElement *elem = itelem->next();
+ int shapeId = elem->getshapeId();
+ int vtkId = elem->getVtkId();
+ if (!shapeIdToVtkIdSet.count(shapeId))
+ {
+ shapeIdToVtkIdSet[shapeId] = emptySet;
+ shapeIds.insert(shapeId);
+ }
+ shapeIdToVtkIdSet[shapeId].insert(vtkId);
+ }
+
+ std::map<int, std::set<DownIdType, DownIdCompare> > shapeIdToEdges; // shapeId --> set of downward edges
+ std::set<DownIdType, DownIdCompare> emptyEdges;
+ emptyEdges.clear();
+
+ std::map<int, std::set<int> >::iterator itShape = shapeIdToVtkIdSet.begin();
+ for (; itShape != shapeIdToVtkIdSet.end(); ++itShape)
+ {
+ int shapeId = itShape->first;
+ MESSAGE(" --- Shape ID --- "<< shapeId);
+ shapeIdToEdges[shapeId] = emptyEdges;
+
+ std::vector<int> nodesEdges;
+
+ std::set<int>::iterator its = itShape->second.begin();
+ for (; its != itShape->second.end(); ++its)
+ {
+ int vtkId = *its;
+ MESSAGE(" " << 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 (neighborsVtkIds[n]<0) // only smds faces are considered as neighbors here
+ continue;
+ int smdsId = meshDS->fromVtkToSmds(neighborsVtkIds[n]);
+ const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
+ if ( shapeIds.count(elem->getshapeId()) && !sgrps->Contains(elem)) // edge : neighbor in the set of shape, not in the group
+ {
+ DownIdType edge(downIds[n], downTypes[n]);
+ if (!shapeIdToEdges[shapeId].count(edge))
+ {
+ shapeIdToEdges[shapeId].insert(edge);
+ int vtkNodeId[3];
+ int nbNodes = grid->getDownArray(downTypes[n])->getNodes(downIds[n],vtkNodeId);
+ nodesEdges.push_back(vtkNodeId[0]);
+ nodesEdges.push_back(vtkNodeId[nbNodes-1]);
+ MESSAGE(" --- nodes " << vtkNodeId[0]+1 << " " << vtkNodeId[nbNodes-1]+1);
+ }
+ }
+ }
+ }
+
+ std::list<int> order;
+ order.clear();
+ if (nodesEdges.size() > 0)
+ {
+ order.push_back(nodesEdges[0]); MESSAGE(" --- back " << order.back()+1); // SMDS id = VTK id + 1;
+ nodesEdges[0] = -1;
+ order.push_back(nodesEdges[1]); MESSAGE(" --- back " << order.back()+1);
+ nodesEdges[1] = -1; // do not reuse this edge
+ bool found = true;
+ while (found)
+ {
+ int nodeTofind = order.back(); // try first to push back
+ int i = 0;
+ for (i = 0; i<nodesEdges.size(); i++)
+ if (nodesEdges[i] == nodeTofind)
+ break;
+ if (i == nodesEdges.size())
+ found = false; // no follower found on back
+ else
+ {
+ if (i%2) // odd ==> use the previous one
+ if (nodesEdges[i-1] < 0)
+ found = false;
+ else
+ {
+ order.push_back(nodesEdges[i-1]); MESSAGE(" --- back " << order.back()+1);
+ nodesEdges[i-1] = -1;
+ }
+ else // even ==> use the next one
+ if (nodesEdges[i+1] < 0)
+ found = false;
+ else
+ {
+ order.push_back(nodesEdges[i+1]); MESSAGE(" --- back " << order.back()+1);
+ nodesEdges[i+1] = -1;
+ }
+ }
+ if (found)
+ continue;
+ // try to push front
+ found = true;
+ nodeTofind = order.front(); // try to push front
+ for (i = 0; i<nodesEdges.size(); i++)
+ if (nodesEdges[i] == nodeTofind)
+ break;
+ if (i == nodesEdges.size())
+ {
+ found = false; // no predecessor found on front
+ continue;
+ }
+ if (i%2) // odd ==> use the previous one
+ if (nodesEdges[i-1] < 0)
+ found = false;
+ else
+ {
+ order.push_front(nodesEdges[i-1]); MESSAGE(" --- front " << order.front()+1);
+ nodesEdges[i-1] = -1;
+ }
+ else // even ==> use the next one
+ if (nodesEdges[i+1] < 0)
+ found = false;
+ else
+ {
+ order.push_front(nodesEdges[i+1]); MESSAGE(" --- front " << order.front()+1);
+ nodesEdges[i+1] = -1;
+ }
+ }
+ }
+
+
+ std::vector<int> nodes;
+ nodes.push_back(shapeId);
+ std::list<int>::iterator itl = order.begin();
+ for (; itl != order.end(); itl++)
+ {
+ nodes.push_back((*itl) + 1); // SMDS id = VTK id + 1;
+ MESSAGE(" ordered node " << nodes[nodes.size()-1]);
+ }
+ listOfListOfNodes.push_back(nodes);
+ }
+
+ // partition geom faces with blocFissure
+ // mesh blocFissure and geom faces of the skin (external wires given, triangle algo to choose)
+ // mesh volume around blocFissure (skin triangles and quadrangle given, tetra algo to choose)
+
+ return;
+}
+
+
+//================================================================================
+/*!
+ * \brief Generates skin mesh (containing 2D cells) from 3D mesh
+ * The created 2D mesh elements based on nodes of free faces of boundary volumes
+ * \return TRUE if operation has been completed successfully, FALSE otherwise
+ */
+//================================================================================
+
+bool SMESH_MeshEditor::Make2DMeshFrom3D()
+{
+ // iterates on volume elements and detect all free faces on them
+ SMESHDS_Mesh* aMesh = GetMeshDS();
+ if (!aMesh)
+ return false;
+ //bool res = false;
+ int nbFree = 0, nbExisted = 0, nbCreated = 0;
+ SMDS_VolumeIteratorPtr vIt = aMesh->volumesIterator();
+ while(vIt->more())
+ {
+ const SMDS_MeshVolume* volume = vIt->next();
+ SMDS_VolumeTool vTool( volume, /*ignoreCentralNodes=*/false );
+ vTool.SetExternalNormal();
+ //const bool isPoly = volume->IsPoly();
+ const int iQuad = volume->IsQuadratic();
+ for ( int iface = 0, n = vTool.NbFaces(); iface < n; iface++ )
+ {
+ if (!vTool.IsFreeFace(iface))
+ continue;
+ nbFree++;
+ vector<const SMDS_MeshNode *> nodes;
+ int nbFaceNodes = vTool.NbFaceNodes(iface);
+ const SMDS_MeshNode** faceNodes = vTool.GetFaceNodes(iface);
+ int inode = 0;
+ for ( ; inode < nbFaceNodes; inode += iQuad+1)
+ nodes.push_back(faceNodes[inode]);
+ if (iQuad) { // add medium nodes
+ for ( inode = 1; inode < nbFaceNodes; inode += 2)
+ nodes.push_back(faceNodes[inode]);
+ if ( nbFaceNodes == 9 ) // bi-quadratic quad
+ nodes.push_back(faceNodes[8]);
+ }
+ // add new face based on volume nodes
+ if (aMesh->FindElement( nodes, SMDSAbs_Face, /*noMedium=*/false) ) {
+ nbExisted++;
+ continue; // face already exsist
+ }
+ AddElement(nodes, SMDSAbs_Face, ( !iQuad && nbFaceNodes/(iQuad+1) > 4 ));
+ nbCreated++;
+ }
+ }
+ return ( nbFree==(nbExisted+nbCreated) );
+}
+
+namespace
+{
+ inline const SMDS_MeshNode* getNodeWithSameID(SMESHDS_Mesh* mesh, const SMDS_MeshNode* node)
+ {
+ if ( const SMDS_MeshNode* n = mesh->FindNode( node->GetID() ))
+ return n;
+ return mesh->AddNodeWithID( node->X(),node->Y(),node->Z(), node->GetID() );
+ }
+}
+//================================================================================
+/*!
+ * \brief Creates missing boundary elements
+ * \param elements - elements whose boundary is to be checked
+ * \param dimension - defines type of boundary elements to create
+ * \param group - a group to store created boundary elements in
+ * \param targetMesh - a mesh to store created boundary elements in
+ * \param toCopyElements - if true, the checked elements will be copied into the targetMesh
+ * \param toCopyExistingBoundary - if true, not only new but also pre-existing
+ * boundary elements will be copied into the targetMesh
+ * \param toAddExistingBondary - if true, not only new but also pre-existing
+ * boundary elements will be added into the new group
+ * \param aroundElements - if true, elements will be created on boundary of given
+ * elements else, on boundary of the whole mesh.
+ * \return nb of added boundary elements
+ */
+//================================================================================
+
+int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
+ Bnd_Dimension dimension,
+ SMESH_Group* group/*=0*/,
+ SMESH_Mesh* targetMesh/*=0*/,
+ bool toCopyElements/*=false*/,
+ bool toCopyExistingBoundary/*=false*/,
+ bool toAddExistingBondary/*= false*/,
+ bool aroundElements/*= false*/)
+{
+ SMDSAbs_ElementType missType = (dimension == BND_2DFROM3D) ? SMDSAbs_Face : SMDSAbs_Edge;
+ SMDSAbs_ElementType elemType = (dimension == BND_1DFROM2D) ? SMDSAbs_Face : SMDSAbs_Volume;
+ // hope that all elements are of the same type, do not check them all
+ if ( !elements.empty() && (*elements.begin())->GetType() != elemType )
+ throw SALOME_Exception(LOCALIZED("wrong element type"));
+
+ if ( !targetMesh )
+ toCopyElements = toCopyExistingBoundary = false;
+
+ SMESH_MeshEditor tgtEditor( targetMesh ? targetMesh : myMesh );
+ SMESHDS_Mesh* aMesh = GetMeshDS(), *tgtMeshDS = tgtEditor.GetMeshDS();
+ int nbAddedBnd = 0;
+
+ // editor adding present bnd elements and optionally holding elements to add to the group
+ SMESH_MeshEditor* presentEditor;
+ SMESH_MeshEditor tgtEditor2( tgtEditor.GetMesh() );
+ presentEditor = toAddExistingBondary ? &tgtEditor : &tgtEditor2;
+
+ SMESH_MesherHelper helper( *myMesh );
+ const TopAbs_ShapeEnum missShapeType = ( missType==SMDSAbs_Face ? TopAbs_FACE : TopAbs_EDGE );
+ SMDS_VolumeTool vTool;
+ TIDSortedElemSet avoidSet;
+ const TIDSortedElemSet emptySet, *elemSet = aroundElements ? &elements : &emptySet;
+ int inode;
+
+ typedef vector<const SMDS_MeshNode*> TConnectivity;
+
+ SMDS_ElemIteratorPtr eIt;
+ if (elements.empty())
+ eIt = aMesh->elementsIterator(elemType);
+ else
+ eIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
+
+ while (eIt->more())
+ {
+ const SMDS_MeshElement* elem = eIt->next();
+ const int iQuad = elem->IsQuadratic();
+
+ // ------------------------------------------------------------------------------------
+ // 1. For an elem, get present bnd elements and connectivities of missing bnd elements
+ // ------------------------------------------------------------------------------------
+ vector<const SMDS_MeshElement*> presentBndElems;
+ vector<TConnectivity> missingBndElems;
+ TConnectivity nodes, elemNodes;
+ if ( vTool.Set(elem, /*ignoreCentralNodes=*/true) ) // elem is a volume --------------
+ {
+ vTool.SetExternalNormal();
+ const SMDS_MeshElement* otherVol = 0;
+ for ( int iface = 0, n = vTool.NbFaces(); iface < n; iface++ )
+ {
+ if ( !vTool.IsFreeFace(iface, &otherVol) &&
+ ( !aroundElements || elements.count( otherVol )))
+ continue;
+ const SMDS_MeshNode** nn = vTool.GetFaceNodes(iface);
+ const int nbFaceNodes = vTool.NbFaceNodes (iface);
+ if ( missType == SMDSAbs_Edge ) // boundary edges
+ {
+ nodes.resize( 2+iQuad );
+ for ( int i = 0; i < nbFaceNodes; i += 1+iQuad)
+ {
+ for ( int j = 0; j < nodes.size(); ++j )
+ nodes[j] =nn[i+j];
+ if ( const SMDS_MeshElement* edge =
+ aMesh->FindElement(nodes,SMDSAbs_Edge,/*noMedium=*/false))
+ presentBndElems.push_back( edge );
+ else
+ missingBndElems.push_back( nodes );
+ }
+ }
+ else // boundary face
+ {
+ nodes.clear();
+ for ( inode = 0; inode < nbFaceNodes; inode += 1+iQuad)
+ nodes.push_back( nn[inode] ); // add corner nodes
+ if (iQuad)
+ for ( inode = 1; inode < nbFaceNodes; inode += 2)
+ nodes.push_back( nn[inode] ); // add medium nodes
+ int iCenter = vTool.GetCenterNodeIndex(iface); // for HEX27
+ if ( iCenter > 0 )
+ nodes.push_back( vTool.GetNodes()[ iCenter ] );
+
+ if (const SMDS_MeshElement * f = aMesh->FindElement( nodes,
+ SMDSAbs_Face, /*noMedium=*/false ))
+ presentBndElems.push_back( f );
+ else
+ missingBndElems.push_back( nodes );
+
+ if ( targetMesh != myMesh )
+ {
+ // add 1D elements on face boundary to be added to a new mesh
+ const SMDS_MeshElement* edge;
+ for ( inode = 0; inode < nbFaceNodes; inode += 1+iQuad)
+ {
+ if ( iQuad )
+ edge = aMesh->FindEdge( nn[inode], nn[inode+1], nn[inode+2]);
+ else
+ edge = aMesh->FindEdge( nn[inode], nn[inode+1]);
+ if ( edge && avoidSet.insert( edge ).second )
+ presentBndElems.push_back( edge );
+ }
+ }
+ }
+ }
+ }
+ else if ( elem->GetType() == SMDSAbs_Face ) // elem is a face ------------------------
+ {
+ avoidSet.clear(), avoidSet.insert( elem );
+ elemNodes.assign( SMDS_MeshElement::iterator( elem->interlacedNodesElemIterator() ),
+ SMDS_MeshElement::iterator() );
+ elemNodes.push_back( elemNodes[0] );
+ nodes.resize( 2 + iQuad );
+ const int nbLinks = elem->NbCornerNodes();
+ for ( int i = 0, iN = 0; i < nbLinks; i++, iN += 1+iQuad )
+ {
+ nodes[0] = elemNodes[iN];
+ nodes[1] = elemNodes[iN+1+iQuad];
+ if ( SMESH_MeshAlgos::FindFaceInSet( nodes[0], nodes[1], *elemSet, avoidSet))
+ continue; // not free link
+
+ if ( iQuad ) nodes[2] = elemNodes[iN+1];
+ if ( const SMDS_MeshElement* edge =
+ aMesh->FindElement(nodes,SMDSAbs_Edge,/*noMedium=*/false))
+ presentBndElems.push_back( edge );
+ else
+ missingBndElems.push_back( nodes );
+ }
+ }
+
+ // ---------------------------------
+ // 2. Add missing boundary elements
+ // ---------------------------------
+ if ( targetMesh != myMesh )
+ // instead of making a map of nodes in this mesh and targetMesh,
+ // we create nodes with same IDs.
+ for ( int i = 0; i < missingBndElems.size(); ++i )
+ {
+ TConnectivity& srcNodes = missingBndElems[i];
+ TConnectivity nodes( srcNodes.size() );
+ for ( inode = 0; inode < nodes.size(); ++inode )
+ nodes[inode] = getNodeWithSameID( tgtMeshDS, srcNodes[inode] );
+ if ( aroundElements && tgtEditor.GetMeshDS()->FindElement( nodes,
+ missType,
+ /*noMedium=*/false))
+ continue;
+ tgtEditor.AddElement(nodes, missType, !iQuad && nodes.size()/(iQuad+1)>4);
+ ++nbAddedBnd;
+ }
+ else
+ for ( int i = 0; i < missingBndElems.size(); ++i )
+ {
+ TConnectivity& nodes = missingBndElems[i];
+ if ( aroundElements && tgtEditor.GetMeshDS()->FindElement( nodes,
+ missType,
+ /*noMedium=*/false))
+ continue;
+ SMDS_MeshElement* elem =
+ tgtEditor.AddElement(nodes, missType, !iQuad && nodes.size()/(iQuad+1)>4);
+ ++nbAddedBnd;
+
+ // try to set a new element to a shape
+ if ( myMesh->HasShapeToMesh() )
+ {
+ bool ok = true;
+ set< pair<TopAbs_ShapeEnum, int > > mediumShapes;
+ const int nbN = nodes.size() / (iQuad+1 );
+ for ( inode = 0; inode < nbN && ok; ++inode )
+ {
+ pair<int, TopAbs_ShapeEnum> i_stype =
+ helper.GetMediumPos( nodes[inode], nodes[(inode+1)%nbN]);
+ if (( ok = ( i_stype.first > 0 && i_stype.second >= TopAbs_FACE )))
+ mediumShapes.insert( make_pair ( i_stype.second, i_stype.first ));
+ }
+ if ( ok && mediumShapes.size() > 1 )
+ {
+ set< pair<TopAbs_ShapeEnum, int > >::iterator stype_i = mediumShapes.begin();
+ pair<TopAbs_ShapeEnum, int> stype_i_0 = *stype_i;
+ for ( ++stype_i; stype_i != mediumShapes.end() && ok; ++stype_i )
+ {
+ if (( ok = ( stype_i->first != stype_i_0.first )))
+ ok = helper.IsSubShape( aMesh->IndexToShape( stype_i->second ),
+ aMesh->IndexToShape( stype_i_0.second ));
+ }
+ }
+ if ( ok && mediumShapes.begin()->first == missShapeType )
+ aMesh->SetMeshElementOnShape( elem, mediumShapes.begin()->second );
+ }
+ }
+
+ // ----------------------------------
+ // 3. Copy present boundary elements
+ // ----------------------------------
+ if ( toCopyExistingBoundary )
+ for ( int i = 0 ; i < presentBndElems.size(); ++i )
+ {
+ const SMDS_MeshElement* e = presentBndElems[i];
+ TConnectivity nodes( e->NbNodes() );
+ for ( inode = 0; inode < nodes.size(); ++inode )
+ nodes[inode] = getNodeWithSameID( tgtMeshDS, e->GetNode(inode) );
+ presentEditor->AddElement(nodes, e->GetType(), e->IsPoly());
+ }
+ else // store present elements to add them to a group
+ for ( int i = 0 ; i < presentBndElems.size(); ++i )
+ {
+ presentEditor->myLastCreatedElems.Append(presentBndElems[i]);
+ }
+
+ } // loop on given elements
+
+ // ---------------------------------------------
+ // 4. Fill group with boundary elements
+ // ---------------------------------------------
+ if ( group )
+ {
+ if ( SMESHDS_Group* g = dynamic_cast<SMESHDS_Group*>( group->GetGroupDS() ))
+ for ( int i = 0; i < tgtEditor.myLastCreatedElems.Size(); ++i )
+ g->SMDSGroup().Add( tgtEditor.myLastCreatedElems( i+1 ));
+ }
+ tgtEditor.myLastCreatedElems.Clear();
+ tgtEditor2.myLastCreatedElems.Clear();
+
+ // -----------------------
+ // 5. Copy given elements
+ // -----------------------
+ if ( toCopyElements && targetMesh != myMesh )
+ {
+ if (elements.empty())
+ eIt = aMesh->elementsIterator(elemType);
+ else
+ eIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
+ while (eIt->more())
+ {
+ const SMDS_MeshElement* elem = eIt->next();
+ TConnectivity nodes( elem->NbNodes() );
+ for ( inode = 0; inode < nodes.size(); ++inode )
+ nodes[inode] = getNodeWithSameID( tgtMeshDS, elem->GetNode(inode) );
+ tgtEditor.AddElement(nodes, elemType, elem->IsPoly());
+
+ tgtEditor.myLastCreatedElems.Clear();
+ }
+ }
+ return nbAddedBnd;
+}