+/*!
+ * \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 describe a partition of the mesh volumes.
+ * The nodes of the internal faces at the boundaries of the groups are doubled.
+ * In option, the internal faces are replaced by flat elements.
+ * Triangles are transformed in prisms, and quadrangles in hexahedrons.
+ * The flat elements are stored in groups of volumes.
+ * @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();
+
+ 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("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
+ {
+ 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)
+
+ 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());
+ }
+ }
+ }
+ }
+ }
+
+ 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;
+
+ 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
+
+ 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();
+
+ 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;
+}
+
+