From: prascle <prascle>
Date: Thu, 2 Feb 2012 17:06:56 +0000 (+0000)
Subject: PR: flat elements on 2D meshes
X-Git-Tag: V6_5_0a1~154
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=fda2d7b99e7b57cfec883aba8c9be0b8d2bd0eff;p=modules%2Fsmesh.git

PR: flat elements on 2D meshes
---

diff --git a/src/SMDS/SMDS_Downward.hxx b/src/SMDS/SMDS_Downward.hxx
index 4f64cfbb5..89207ae08 100644
--- a/src/SMDS/SMDS_Downward.hxx
+++ b/src/SMDS/SMDS_Downward.hxx
@@ -78,6 +78,7 @@ public:
   virtual const unsigned char* getUpTypes(int cellId) = 0;
   virtual void getNodeIds(int cellId, std::set<int>& nodeSet) = 0;
   virtual int getNodes(int cellId, int* nodevec) {return 0; }
+  virtual void getOrderedNodesOfFace(int cellId, std::vector<vtkIdType>& orderedNodes) {};
   int getVtkCellId(int cellId)
   {
     return _vtkCellIds[cellId];
@@ -178,7 +179,6 @@ public:
   virtual const int* getUpCells(int cellId);
   virtual const unsigned char* getUpTypes(int cellId);
   virtual void getNodeIds(int cellId, std::set<int>& nodeSet);
-  virtual void getOrderedNodesOfFace(int cellId, std::vector<vtkIdType>& orderedNodes) = 0;
 protected:
   SMDS_Down3D(SMDS_UnstructuredGrid *grid, int nbDownCells);
   ~SMDS_Down3D();
diff --git a/src/SMDS/SMDS_Mesh.cxx b/src/SMDS/SMDS_Mesh.cxx
index 02ac1026c..4afb4da46 100644
--- a/src/SMDS/SMDS_Mesh.cxx
+++ b/src/SMDS/SMDS_Mesh.cxx
@@ -1507,6 +1507,53 @@ SMDS_MeshVolume* SMDS_Mesh::AddVolumeFromVtkIdsWithID(const std::vector<vtkIdTyp
   return volvtk;
 }
 
+SMDS_MeshFace* SMDS_Mesh::AddFaceFromVtkIds(const std::vector<vtkIdType>& vtkNodeIds)
+{
+  int ID = myElementIDFactory->GetFreeID();
+  SMDS_MeshFace * f = SMDS_Mesh::AddFaceFromVtkIdsWithID(vtkNodeIds, ID);
+  if (f == NULL) myElementIDFactory->ReleaseID(ID);
+  return f;
+}
+
+SMDS_MeshFace* SMDS_Mesh::AddFaceFromVtkIdsWithID(const std::vector<vtkIdType>& vtkNodeIds, const int ID)
+{
+  SMDS_VtkFace *facevtk = myFacePool->getNew();
+  facevtk->init(vtkNodeIds, this);
+  if (!this->registerElement(ID,facevtk))
+    {
+      this->myGrid->GetCellTypesArray()->SetValue(facevtk->getVtkId(), VTK_EMPTY_CELL);
+      myFacePool->destroy(facevtk);
+      return 0;
+    }
+  adjustmyCellsCapacity(ID);
+  myCells[ID] = facevtk;
+  vtkIdType aVtkType = facevtk->GetVtkType();
+  switch (aVtkType)
+  {
+    case VTK_TRIANGLE:
+      myInfo.myNbTriangles++;
+      break;
+    case VTK_QUAD:
+      myInfo.myNbQuadrangles++;
+      break;
+    case VTK_QUADRATIC_TRIANGLE:
+      myInfo.myNbQuadTriangles++;
+      break;
+    case VTK_QUADRATIC_QUAD:
+      myInfo.myNbQuadQuadrangles++;
+      break;
+    case VTK_BIQUADRATIC_QUAD:
+      myInfo.myNbBiQuadQuadrangles++;
+      break;
+    case VTK_POLYGON:
+      myInfo.myNbPolygons++;
+      break;
+     default:
+      myInfo.myNbPolygons++;
+  }
+  return facevtk;
+}
+
 ///////////////////////////////////////////////////////////////////////////////
 /// Registers element with the given ID, maintains inverse connections
 ///////////////////////////////////////////////////////////////////////////////
diff --git a/src/SMDS/SMDS_Mesh.hxx b/src/SMDS/SMDS_Mesh.hxx
index 91f75d9b7..c9b6468fc 100644
--- a/src/SMDS/SMDS_Mesh.hxx
+++ b/src/SMDS/SMDS_Mesh.hxx
@@ -568,6 +568,11 @@ public:
   virtual SMDS_MeshVolume* AddVolumeFromVtkIdsWithID(const std::vector<vtkIdType>& vtkNodeIds,
                                                      const int ID);
 
+  virtual SMDS_MeshFace* AddFaceFromVtkIds(const std::vector<vtkIdType>& vtkNodeIds);
+
+  virtual SMDS_MeshFace* AddFaceFromVtkIdsWithID(const std::vector<vtkIdType>& vtkNodeIds,
+                                                     const int ID);
+
   virtual void RemoveElement(const SMDS_MeshElement *        elem,
                              std::list<const SMDS_MeshElement *>& removedElems,
                              std::list<const SMDS_MeshElement *>& removedNodes,
diff --git a/src/SMDS/SMDS_UnstructuredGrid.cxx b/src/SMDS/SMDS_UnstructuredGrid.cxx
index 2902856bc..2ac8afe63 100644
--- a/src/SMDS/SMDS_UnstructuredGrid.cxx
+++ b/src/SMDS/SMDS_UnstructuredGrid.cxx
@@ -929,15 +929,17 @@ void SMDS_UnstructuredGrid::ModifyCellNodes(int vtkVolId, std::map<int, int> loc
  * @param orderedNodes list of nodes to reorder (in out)
  * @return size of the list
  */
-int SMDS_UnstructuredGrid::getOrderedNodesOfFace(int vtkVolId, std::vector<vtkIdType>& orderedNodes)
+int SMDS_UnstructuredGrid::getOrderedNodesOfFace(int vtkVolId, int& dim, std::vector<vtkIdType>& orderedNodes)
 {
   int vtkType = this->GetCellType(vtkVolId);
-  int cellDim = SMDS_Downward::getCellDimension(vtkType);
-  if (cellDim != 3)
-    return 0;
-  SMDS_Down3D *downvol = static_cast<SMDS_Down3D*> (_downArray[vtkType]);
-  int downVolId = this->_cellIdToDownId[vtkVolId];
-  downvol->getOrderedNodesOfFace(downVolId, orderedNodes);
+  dim = SMDS_Downward::getCellDimension(vtkType);
+  if (dim == 3)
+    {
+      SMDS_Down3D *downvol = static_cast<SMDS_Down3D*> (_downArray[vtkType]);
+      int downVolId = this->_cellIdToDownId[vtkVolId];
+      downvol->getOrderedNodesOfFace(downVolId, orderedNodes);
+    }
+  // else nothing to do;
   return orderedNodes.size();
 }
 
@@ -970,7 +972,7 @@ void SMDS_UnstructuredGrid::BuildLinks()
  * @param nodeDomains: map(original id --> map(domain --> duplicated node id))
  * @return ok if success.
  */
-SMDS_MeshVolume* SMDS_UnstructuredGrid::extrudeVolumeFromFace(int vtkVolId,
+SMDS_MeshCell* SMDS_UnstructuredGrid::extrudeVolumeFromFace(int vtkVolId,
                                                   int domain1,
                                                   int domain2,
                                                   std::set<int>& originalNodes,
@@ -984,65 +986,89 @@ SMDS_MeshVolume* SMDS_UnstructuredGrid::extrudeVolumeFromFace(int vtkVolId,
   for (; it != originalNodes.end(); ++it)
     orderedOriginals.push_back(*it);
 
-  int nbNodes = this->getOrderedNodesOfFace(vtkVolId, orderedOriginals);
+  int dim = 0;
+  int nbNodes = this->getOrderedNodesOfFace(vtkVolId, dim, orderedOriginals);
   vector<vtkIdType> orderedNodes;
 
+  bool isQuadratic = false;
   switch (orderedOriginals.size())
   {
     case 3:
-    case 4:
-      for (int i = 0; i < nbNodes; i++)
-        orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
-      for (int i = 0; i < nbNodes; i++)
-        orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
+      if (dim == 2)
+        isQuadratic = true;
       break;
     case 6:
     case 8:
-      {
-        long dom1 = domain1;
-        long dom2 = domain2;
-        long dom1_2; // for nodeQuadDomains
-        if (domain1 < domain2)
-          dom1_2 = dom1 + INT_MAX * dom2;
-        else
-          dom1_2 = dom2 + INT_MAX * dom1;
-        //cerr << "dom1=" << dom1 << " dom2=" << dom2 << " dom1_2=" << dom1_2 << endl;
-        int ima = orderedOriginals.size();
-        int mid = orderedOriginals.size() / 2;
-        //cerr << "ima=" << ima << " mid=" << mid << endl;
-        for (int i = 0; i < mid; i++)
-          orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
-        for (int i = 0; i < mid; i++)
-          orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
-        for (int i = mid; i < ima; i++)
-          orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
-        for (int i = mid; i < ima; i++)
-          orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
-        for (int i = 0; i < mid; i++)
-          {
-            int oldId = orderedOriginals[i];
-            int newId;
-            if (nodeQuadDomains.count(oldId) && nodeQuadDomains[oldId].count(dom1_2))
-              newId = nodeQuadDomains[oldId][dom1_2];
-            else
-              {
-                double *coords = this->GetPoint(oldId);
-                SMDS_MeshNode *newNode = _mesh->AddNode(coords[0], coords[1], coords[2]);
-                newId = newNode->getVtkId();
-                std::map<long, int> emptyMap;
-                nodeQuadDomains[oldId] = emptyMap;
-                nodeQuadDomains[oldId][dom1_2] = newId;
-              }
-            orderedNodes.push_back(newId);
-          }
-      }
+      isQuadratic = true;
       break;
     default:
-      ASSERT(0);
+      isQuadratic = false;
+      break;
   }
 
-  SMDS_MeshVolume *vol = _mesh->AddVolumeFromVtkIds(orderedNodes);
+  if (isQuadratic)
+    {
+      long dom1 = domain1;
+      long dom2 = domain2;
+      long dom1_2; // for nodeQuadDomains
+      if (domain1 < domain2)
+        dom1_2 = dom1 + INT_MAX * dom2;
+      else
+        dom1_2 = dom2 + INT_MAX * dom1;
+      //cerr << "dom1=" << dom1 << " dom2=" << dom2 << " dom1_2=" << dom1_2 << endl;
+      int ima = orderedOriginals.size();
+      int mid = orderedOriginals.size() / 2;
+      //cerr << "ima=" << ima << " mid=" << mid << endl;
+      for (int i = 0; i < mid; i++)
+        orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
+      for (int i = 0; i < mid; i++)
+        orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
+      for (int i = mid; i < ima; i++)
+        orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
+      for (int i = mid; i < ima; i++)
+        orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
+      for (int i = 0; i < mid; i++)
+        {
+          int oldId = orderedOriginals[i];
+          int newId;
+          if (nodeQuadDomains.count(oldId) && nodeQuadDomains[oldId].count(dom1_2))
+            newId = nodeQuadDomains[oldId][dom1_2];
+          else
+            {
+              double *coords = this->GetPoint(oldId);
+              SMDS_MeshNode *newNode = _mesh->AddNode(coords[0], coords[1], coords[2]);
+              newId = newNode->getVtkId();
+              std::map<long, int> emptyMap;
+              nodeQuadDomains[oldId] = emptyMap;
+              nodeQuadDomains[oldId][dom1_2] = newId;
+            }
+          orderedNodes.push_back(newId);
+        }
+    }
+  else
+    {
+      for (int i = 0; i < nbNodes; i++)
+        orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
+      if (dim == 3)
+        for (int i = 0; i < nbNodes; i++)
+          orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
+      else
+        for (int i = nbNodes-1; i >= 0; i--)
+          orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
+
+    }
+
+  if (dim == 3)
+    {
+      SMDS_MeshVolume *vol = _mesh->AddVolumeFromVtkIds(orderedNodes);
+      return vol;
+    }
+  else if (dim == 2)
+    {
+      SMDS_MeshFace *face = _mesh->AddFaceFromVtkIds(orderedNodes);
+      return face;
+    }
 
   // TODO update sub-shape list of elements and nodes
-  return vol;
+  return 0;
 }
diff --git a/src/SMDS/SMDS_UnstructuredGrid.hxx b/src/SMDS/SMDS_UnstructuredGrid.hxx
index d2489a020..370aec2e9 100644
--- a/src/SMDS/SMDS_UnstructuredGrid.hxx
+++ b/src/SMDS/SMDS_UnstructuredGrid.hxx
@@ -48,6 +48,7 @@
 
 class SMDS_Downward;
 class SMDS_Mesh;
+class SMDS_MeshCell;
 class SMDS_MeshVolume;
 
 class SMDS_EXPORT SMDS_CellLinks: public vtkCellLinks
@@ -86,11 +87,11 @@ public:
   int GetParentVolumes(int* volVtkIds, int downId, unsigned char downType);
   void GetNodeIds(std::set<int>& nodeSet, int downId, unsigned char downType);
   void ModifyCellNodes(int vtkVolId, std::map<int, int> localClonedNodeIds);
-  int getOrderedNodesOfFace(int vtkVolId, std::vector<vtkIdType>& orderedNodes);
+  int getOrderedNodesOfFace(int vtkVolId, int& dim, std::vector<vtkIdType>& orderedNodes);
   void BuildLinks();
-  SMDS_MeshVolume* extrudeVolumeFromFace(int vtkVolId, int domain1, int domain2, std::set<int>& originalNodes,
-                                         std::map<int, std::map<int, int> >& nodeDomains,
-                                         std::map<int, std::map<long,int> >& nodeQuadDomains);
+  SMDS_MeshCell* extrudeVolumeFromFace(int vtkVolId, int domain1, int domain2, std::set<int>& originalNodes,
+                                       std::map<int, std::map<int, int> >& nodeDomains,
+                                       std::map<int, std::map<long,int> >& nodeQuadDomains);
   vtkCellLinks* GetLinks()
   {
     return Links;
diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx
index a466c3c82..9e9768530 100644
--- a/src/SMESH/SMESH_MeshEditor.cxx
+++ b/src/SMESH/SMESH_MeshEditor.cxx
@@ -10727,6 +10727,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       //     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 is domain.
 
+      //MESSAGE("Domain " << idom);
       const TIDSortedElemSet& domain = theElems[idom];
       TIDSortedElemSet::const_iterator elemItr = domain.begin();
       for (; elemItr != domain.end(); ++elemItr)
@@ -10735,6 +10736,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
           if (!anElem)
             continue;
           int vtkId = anElem->getVtkId();
+          //MESSAGE("  vtkId " << vtkId << " smdsId " << anElem->GetID());
           int neighborsVtkIds[NBMAXNEIGHBORS];
           int downIds[NBMAXNEIGHBORS];
           unsigned char downTypes[NBMAXNEIGHBORS];
@@ -10752,6 +10754,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                     {
                       faceDomains[face][idom] = vtkId; // volume associated to face in this domain
                       celldom[vtkId] = idom;
+                      //MESSAGE("       cell with a border " << vtkId << " domain " << idom);
                     }
                 }
             }
@@ -10767,6 +10770,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
 
   for (int idomain = 0; idomain < theElems.size(); idomain++)
     {
+      //MESSAGE("Domain " << idomain);
       const TIDSortedElemSet& domain = theElems[idomain];
       itface = faceDomains.begin();
       for (; itface != faceDomains.end(); ++itface)
@@ -10806,6 +10810,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                     continue;
                   cellDomains[aCell][idomain] = vtkId;
                   celldom[vtkId] = idomain;
+                  //MESSAGE("       cell " << vtkId << " domain " << idomain);
                 }
             }
         }
@@ -10820,7 +10825,8 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   //     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 muti domains with domain order
+  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++)
     {
@@ -10890,76 +10896,89 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
               //MESSAGE("multiple Nodes detected on a shared face");
               int downId = itface->first.cellId;
               unsigned char cellType = itface->first.cellType;
-              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++)
+              // --- 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(edgeType[ie])->getNodes(downEdgeIds[ie], nodes);
-                  if (mutipleNodes.count(nodes[0]) && mutipleNodes.count(nodes[nbNodes-1]))
+                  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++)
                     {
-                      vector<int> vn0 = mutipleNodes[nodes[0]];
-                      vector<int> vn1 = mutipleNodes[nodes[nbNodes - 1]];
-                      sort( vn0.begin(), vn0.end() );
-                      sort( vn1.begin(), vn1.end() );
-                      if (vn0 == vn1)
+                      int nodes[3];
+                      int nbNodes = grid->getDownArray(edgeType[ie])->getNodes(downEdgeIds[ie], nodes);
+                      if (mutipleNodes.count(nodes[0]) && mutipleNodes.count(nodes[nbNodes-1]))
                         {
-                          //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 < vn0.size(); id++)
+                          vector<int> vn0 = mutipleNodes[nodes[0]];
+                          vector<int> vn1 = mutipleNodes[nodes[nbNodes - 1]];
+                          sort( vn0.begin(), vn0.end() );
+                          sort( vn1.begin(), vn1.end() );
+                          if (vn0 == vn1)
                             {
-                              int idom = vn0[id];
-                              for (int ivol=0; ivol<nbvol; ivol++)
+                              //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 < vn0.size(); id++)
                                 {
-                                  int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
-                                  SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
-                                  if (theElems[idom].count(elem))
+                                  int idom = vn0[id];
+                                  for (int ivol=0; ivol<nbvol; ivol++)
                                     {
-                                      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
+                                      int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
+                                      SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
+                                      if (theElems[idom].count(elem))
                                         {
-                                          gp_Pnt g(values[0], values[1], values[2]);
-                                          angleDom[idom] = OrientedAngle(p0, p1, gref, g); // -pi<angle<+pi
-                                          //MESSAGE("  angle=" << angleDom[idom]);
+                                          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;
                                         }
-                                      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
                             }
-                          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
                         }
                     }
                 }
@@ -10994,7 +11013,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
           int vtkVolId = itdom->second;
           itdom++;
           int dom2 = itdom->first;
-          SMDS_MeshVolume *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
+          SMDS_MeshCell *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
                                                              nodeQuadDomains);
           stringstream grpname;
           grpname << "j_";
@@ -11005,7 +11024,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
           int idg;
           string namegrp = grpname.str();
           if (!mapOfJunctionGroups.count(namegrp))
-            mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+            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());
@@ -11013,10 +11032,37 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
     }
 
   // --- 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)
         {
@@ -11048,7 +11094,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
             }
           else
             {
-              //MESSAGE("Quadratic multiple joints not implemented");
+              MESSAGE("Quadratic multiple joints not implemented");
               // TODO quadratic nodes
             }
         }
diff --git a/src/SMESH_I/SMESH_MeshEditor_i.cxx b/src/SMESH_I/SMESH_MeshEditor_i.cxx
index ac3973702..3fcfabc58 100644
--- a/src/SMESH_I/SMESH_MeshEditor_i.cxx
+++ b/src/SMESH_I/SMESH_MeshEditor_i.cxx
@@ -5739,8 +5739,8 @@ CORBA::Boolean SMESH_MeshEditor_i::DoubleNodesOnGroupBoundaries( const SMESH::Li
     SMESH::SMESH_GroupBase_var aGrp = theDomains[ i ];
     if ( !CORBA::is_nil( aGrp ) /*&& ( aGrp->GetType() != SMESH::NODE )*/ )
     {
-      if ( aGrp->GetType() != SMESH::VOLUME )
-        THROW_SALOME_CORBA_EXCEPTION("Not a volume group", SALOME::BAD_PARAM);
+//      if ( aGrp->GetType() != SMESH::VOLUME )
+//        THROW_SALOME_CORBA_EXCEPTION("Not a volume group", SALOME::BAD_PARAM);
       TIDSortedElemSet domain;
       domain.clear();
       domains.push_back(domain);