Salome HOME
Merge from V6_5_BR 05/06/2012
[modules/smesh.git] / src / SMDS / SMDS_UnstructuredGrid.cxx
index a7837d4934857ced24809dfa38aaab838c339938..099b8cbe85235b332d2cdfb89df257bde51bfde0 100644 (file)
@@ -1,8 +1,28 @@
+// Copyright (C) 2010-2012  CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
 #define CHRONODEF
 #include "SMDS_UnstructuredGrid.hxx"
 #include "SMDS_Mesh.hxx"
 #include "SMDS_MeshInfo.hxx"
 #include "SMDS_Downward.hxx"
+#include "SMDS_MeshVolume.hxx"
 
 #include "utilities.h"
 
@@ -12,6 +32,7 @@
 #include <vtkUnsignedCharArray.h>
 
 #include <list>
+#include <climits>
 
 using namespace std;
 
@@ -310,7 +331,12 @@ void SMDS_UnstructuredGrid::copyBloc(vtkUnsignedCharArray *newTypes, std::vector
 
 int SMDS_UnstructuredGrid::CellIdToDownId(int vtkCellId)
 {
-  // ASSERT((vtkCellId >= 0) && (vtkCellId < _cellIdToDownId.size()));
+  if((vtkCellId < 0) || (vtkCellId >= _cellIdToDownId.size()))
+    {
+      //MESSAGE("SMDS_UnstructuredGrid::CellIdToDownId structure not up to date: vtkCellId="
+      //    << vtkCellId << " max="<< _cellIdToDownId.size());
+      return -1;
+    }
   return _cellIdToDownId[vtkCellId];
 }
 
@@ -320,6 +346,17 @@ void SMDS_UnstructuredGrid::setCellIdToDownId(int vtkCellId, int downId)
   _cellIdToDownId[vtkCellId] = downId;
 }
 
+void SMDS_UnstructuredGrid::CleanDownwardConnectivity()
+{
+  for (int i = 0; i < _downArray.size(); i++)
+    {
+      if (_downArray[i])
+        delete _downArray[i];
+      _downArray[i] = 0;
+    }
+  _cellIdToDownId.clear();
+}
+
 /*! Build downward connectivity: to do only when needed because heavy memory load.
  *  Downward connectivity is no more valid if vtkUnstructuredGrid is modified.
  *
@@ -331,83 +368,87 @@ void SMDS_UnstructuredGrid::BuildDownwardConnectivity(bool withEdges)
 
   // --- erase previous data if any
 
-  for (int i = 0; i < _downArray.size(); i++)
-    {
-      if (_downArray[i])
-        delete _downArray[i];
-      _downArray[i] = 0;
-    }
-  _cellIdToDownId.clear();
+  this->CleanDownwardConnectivity();
 
   // --- create SMDS_Downward structures (in _downArray vector[vtkCellType])
 
-  _downArray.resize(VTK_MAXTYPE + 1, 0); // --- max. type value = VTK_QUADRATIC_PYRAMID
-
-  _downArray[VTK_LINE] = new SMDS_DownEdge(this);
-  _downArray[VTK_QUADRATIC_EDGE] = new SMDS_DownQuadEdge(this);
-  _downArray[VTK_TRIANGLE] = new SMDS_DownTriangle(this);
-  _downArray[VTK_QUADRATIC_TRIANGLE] = new SMDS_DownQuadTriangle(this);
-  _downArray[VTK_QUAD] = new SMDS_DownQuadrangle(this);
-  _downArray[VTK_QUADRATIC_QUAD] = new SMDS_DownQuadQuadrangle(this);
-  _downArray[VTK_TETRA] = new SMDS_DownTetra(this);
-  _downArray[VTK_QUADRATIC_TETRA] = new SMDS_DownQuadTetra(this);
-  _downArray[VTK_PYRAMID] = new SMDS_DownPyramid(this);
-  _downArray[VTK_QUADRATIC_PYRAMID] = new SMDS_DownQuadPyramid(this);
-  _downArray[VTK_WEDGE] = new SMDS_DownPenta(this);
-  _downArray[VTK_QUADRATIC_WEDGE] = new SMDS_DownQuadPenta(this);
-  _downArray[VTK_HEXAHEDRON] = new SMDS_DownHexa(this);
-  _downArray[VTK_QUADRATIC_HEXAHEDRON] = new SMDS_DownQuadHexa(this);
+  _downArray.resize(VTK_MAXTYPE + 1, 0);
+
+  _downArray[VTK_LINE]                    = new SMDS_DownEdge(this);
+  _downArray[VTK_QUADRATIC_EDGE]          = new SMDS_DownQuadEdge(this);
+  _downArray[VTK_TRIANGLE]                = new SMDS_DownTriangle(this);
+  _downArray[VTK_QUADRATIC_TRIANGLE]      = new SMDS_DownQuadTriangle(this);
+  _downArray[VTK_QUAD]                    = new SMDS_DownQuadrangle(this);
+  _downArray[VTK_QUADRATIC_QUAD]          = new SMDS_DownQuadQuadrangle(this);
+  _downArray[VTK_BIQUADRATIC_QUAD]        = new SMDS_DownQuadQuadrangle(this);
+  _downArray[VTK_TETRA]                   = new SMDS_DownTetra(this);
+  _downArray[VTK_QUADRATIC_TETRA]         = new SMDS_DownQuadTetra(this);
+  _downArray[VTK_PYRAMID]                 = new SMDS_DownPyramid(this);
+  _downArray[VTK_QUADRATIC_PYRAMID]       = new SMDS_DownQuadPyramid(this);
+  _downArray[VTK_WEDGE]                   = new SMDS_DownPenta(this);
+  _downArray[VTK_QUADRATIC_WEDGE]         = new SMDS_DownQuadPenta(this);
+  _downArray[VTK_HEXAHEDRON]              = new SMDS_DownHexa(this);
+  _downArray[VTK_QUADRATIC_HEXAHEDRON]    = new SMDS_DownQuadHexa(this);
+  _downArray[VTK_TRIQUADRATIC_HEXAHEDRON] = new SMDS_DownQuadHexa(this);
+  _downArray[VTK_HEXAGONAL_PRISM]         = new SMDS_DownPenta(this);
 
   // --- get detailed info of number of cells of each type, allocate SMDS_downward structures
 
   const SMDS_MeshInfo &meshInfo = _mesh->GetMeshInfo();
 
-  int nbLinTetra = meshInfo.NbTetras(ORDER_LINEAR);
-  int nbQuadTetra = meshInfo.NbTetras(ORDER_QUADRATIC);
-  int nbLinPyra = meshInfo.NbPyramids(ORDER_LINEAR);
-  int nbQuadPyra = meshInfo.NbPyramids(ORDER_QUADRATIC);
-  int nbLinPrism = meshInfo.NbPrisms(ORDER_LINEAR);
-  int nbQuadPrism = meshInfo.NbPrisms(ORDER_QUADRATIC);
-  int nbLinHexa = meshInfo.NbHexas(ORDER_LINEAR);
-  int nbQuadHexa = meshInfo.NbHexas(ORDER_QUADRATIC);
-
-  int nbLineGuess = int((4.0 / 3.0) * nbLinTetra + 2 * nbLinPrism + 2.5 * nbLinPyra + 3 * nbLinHexa);
+  int nbLinTetra  = meshInfo.NbTetras  (ORDER_LINEAR);
+  int nbQuadTetra = meshInfo.NbTetras  (ORDER_QUADRATIC);
+  int nbLinPyra   = meshInfo.NbPyramids(ORDER_LINEAR);
+  int nbQuadPyra  = meshInfo.NbPyramids(ORDER_QUADRATIC);
+  int nbLinPrism  = meshInfo.NbPrisms  (ORDER_LINEAR);
+  int nbQuadPrism = meshInfo.NbPrisms  (ORDER_QUADRATIC);
+  int nbLinHexa   = meshInfo.NbHexas   (ORDER_LINEAR);
+  int nbQuadHexa  = meshInfo.NbHexas   (ORDER_QUADRATIC);
+  int nbHexPrism  = meshInfo.NbHexPrisms();
+
+  int nbLineGuess     = int((4.0 / 3.0) * nbLinTetra + 2 * nbLinPrism + 2.5 * nbLinPyra + 3 * nbLinHexa);
   int nbQuadEdgeGuess = int((4.0 / 3.0) * nbQuadTetra + 2 * nbQuadPrism + 2.5 * nbQuadPyra + 3 * nbQuadHexa);
-  int nbLinTriaGuess = 2 * nbLinTetra + nbLinPrism + 2 * nbLinPyra;
+  int nbLinTriaGuess  = 2 * nbLinTetra + nbLinPrism + 2 * nbLinPyra;
   int nbQuadTriaGuess = 2 * nbQuadTetra + nbQuadPrism + 2 * nbQuadPyra;
-  int nbLinQuadGuess = int((2.0 / 3.0) * nbLinPrism + (1.0 / 2.0) * nbLinPyra + 3 * nbLinHexa);
+  int nbLinQuadGuess  = int((2.0 / 3.0) * nbLinPrism + (1.0 / 2.0) * nbLinPyra + 3 * nbLinHexa);
   int nbQuadQuadGuess = int((2.0 / 3.0) * nbQuadPrism + (1.0 / 2.0) * nbQuadPyra + 3 * nbQuadHexa);
 
-  int GuessSize[VTK_QUADRATIC_TETRA];
-  GuessSize[VTK_LINE] = nbLineGuess;
-  GuessSize[VTK_QUADRATIC_EDGE] = nbQuadEdgeGuess;
-  GuessSize[VTK_TRIANGLE] = nbLinTriaGuess;
-  GuessSize[VTK_QUADRATIC_TRIANGLE] = nbQuadTriaGuess;
-  GuessSize[VTK_QUAD] = nbLinQuadGuess;
-  GuessSize[VTK_QUADRATIC_QUAD] = nbQuadQuadGuess;
-  GuessSize[VTK_TETRA] = nbLinTetra;
-  GuessSize[VTK_QUADRATIC_TETRA] = nbQuadTetra;
-  GuessSize[VTK_PYRAMID] = nbLinPyra;
-  GuessSize[VTK_QUADRATIC_PYRAMID] = nbQuadPyra;
-  GuessSize[VTK_WEDGE] = nbLinPrism;
-  GuessSize[VTK_QUADRATIC_WEDGE] = nbQuadPrism;
-  GuessSize[VTK_HEXAHEDRON] = nbLinHexa;
-  GuessSize[VTK_QUADRATIC_HEXAHEDRON] = nbQuadHexa;
-
-  _downArray[VTK_LINE]->allocate(nbLineGuess);
-  _downArray[VTK_QUADRATIC_EDGE]->allocate(nbQuadEdgeGuess);
-  _downArray[VTK_TRIANGLE]->allocate(nbLinTriaGuess);
-  _downArray[VTK_QUADRATIC_TRIANGLE]->allocate(nbQuadTriaGuess);
-  _downArray[VTK_QUAD]->allocate(nbLinQuadGuess);
-  _downArray[VTK_QUADRATIC_QUAD]->allocate(nbQuadQuadGuess);
-  _downArray[VTK_TETRA]->allocate(nbLinTetra);
-  _downArray[VTK_QUADRATIC_TETRA]->allocate(nbQuadTetra);
-  _downArray[VTK_PYRAMID]->allocate(nbLinPyra);
-  _downArray[VTK_QUADRATIC_PYRAMID]->allocate(nbQuadPyra);
-  _downArray[VTK_WEDGE]->allocate(nbLinPrism);
-  _downArray[VTK_QUADRATIC_WEDGE]->allocate(nbQuadPrism);
-  _downArray[VTK_HEXAHEDRON]->allocate(nbLinHexa);
-  _downArray[VTK_QUADRATIC_HEXAHEDRON]->allocate(nbQuadHexa);
+  int GuessSize[VTK_MAXTYPE];
+  GuessSize[VTK_LINE]                    = nbLineGuess;
+  GuessSize[VTK_QUADRATIC_EDGE]          = nbQuadEdgeGuess;
+  GuessSize[VTK_TRIANGLE]                = nbLinTriaGuess;
+  GuessSize[VTK_QUADRATIC_TRIANGLE]      = nbQuadTriaGuess;
+  GuessSize[VTK_QUAD]                    = nbLinQuadGuess;
+  GuessSize[VTK_QUADRATIC_QUAD]          = nbQuadQuadGuess;
+  GuessSize[VTK_BIQUADRATIC_QUAD]        = nbQuadQuadGuess;
+  GuessSize[VTK_TETRA]                   = nbLinTetra;
+  GuessSize[VTK_QUADRATIC_TETRA]         = nbQuadTetra;
+  GuessSize[VTK_PYRAMID]                 = nbLinPyra;
+  GuessSize[VTK_QUADRATIC_PYRAMID]       = nbQuadPyra;
+  GuessSize[VTK_WEDGE]                   = nbLinPrism;
+  GuessSize[VTK_QUADRATIC_WEDGE]         = nbQuadPrism;
+  GuessSize[VTK_HEXAHEDRON]              = nbLinHexa;
+  GuessSize[VTK_QUADRATIC_HEXAHEDRON]    = nbQuadHexa;
+  GuessSize[VTK_TRIQUADRATIC_HEXAHEDRON] = nbQuadHexa;
+  GuessSize[VTK_HEXAGONAL_PRISM]         = nbHexPrism;
+
+  _downArray[VTK_LINE]                   ->allocate(nbLineGuess);
+  _downArray[VTK_QUADRATIC_EDGE]         ->allocate(nbQuadEdgeGuess);
+  _downArray[VTK_TRIANGLE]               ->allocate(nbLinTriaGuess);
+  _downArray[VTK_QUADRATIC_TRIANGLE]     ->allocate(nbQuadTriaGuess);
+  _downArray[VTK_QUAD]                   ->allocate(nbLinQuadGuess);
+  _downArray[VTK_QUADRATIC_QUAD]         ->allocate(nbQuadQuadGuess);
+  _downArray[VTK_BIQUADRATIC_QUAD]       ->allocate(nbQuadQuadGuess);
+  _downArray[VTK_TETRA]                  ->allocate(nbLinTetra);
+  _downArray[VTK_QUADRATIC_TETRA]        ->allocate(nbQuadTetra);
+  _downArray[VTK_PYRAMID]                ->allocate(nbLinPyra);
+  _downArray[VTK_QUADRATIC_PYRAMID]      ->allocate(nbQuadPyra);
+  _downArray[VTK_WEDGE]                  ->allocate(nbLinPrism);
+  _downArray[VTK_QUADRATIC_WEDGE]        ->allocate(nbQuadPrism);
+  _downArray[VTK_HEXAHEDRON]             ->allocate(nbLinHexa);
+  _downArray[VTK_QUADRATIC_HEXAHEDRON]   ->allocate(nbQuadHexa);
+  _downArray[VTK_TRIQUADRATIC_HEXAHEDRON]->allocate(nbQuadHexa);
+  _downArray[VTK_HEXAGONAL_PRISM]        ->allocate(nbHexPrism);
 
   // --- iteration on vtkUnstructuredGrid cells, only faces
   //     for each vtk face:
@@ -740,6 +781,115 @@ int SMDS_UnstructuredGrid::GetNeighbors(int* neighborsVtkIds, int* downIds, unsi
   return nb;
 }
 
+/*! get the volumes containing a face or an edge of the grid
+ * The edge or face belongs to the vtkUnstructuredGrid
+ * @param volVtkIds vector of parent volume ids to fill (reserve enough space!)
+ * @param vtkId vtk id of the face or edge
+ */
+int SMDS_UnstructuredGrid::GetParentVolumes(int* volVtkIds, int vtkId)
+{
+  int vtkType = this->GetCellType(vtkId);
+  int dim = SMDS_Downward::getCellDimension(vtkType);
+  int nbFaces = 0;
+  unsigned char cellTypes[1000];
+  int downCellId[1000];
+  if (dim == 1)
+    {
+      int downId = this->CellIdToDownId(vtkId);
+      if (downId < 0)
+        {
+          MESSAGE("Downward structure not up to date: new edge not taken into account");
+          return 0;
+        }
+      nbFaces = _downArray[vtkType]->getNumberOfUpCells(downId);
+      const int *upCells = _downArray[vtkType]->getUpCells(downId);
+      const unsigned char* upTypes = _downArray[vtkType]->getUpTypes(downId);
+      for (int i=0; i< nbFaces; i++)
+        {
+          cellTypes[i] = upTypes[i];
+          downCellId[i] = upCells[i];
+        }
+    }
+  else if (dim == 2)
+    {
+      nbFaces = 1;
+      cellTypes[0] = this->GetCellType(vtkId);
+      int downId = this->CellIdToDownId(vtkId);
+      if (downId < 0)
+        {
+          MESSAGE("Downward structure not up to date: new face not taken into account");
+          return 0;
+        }
+      downCellId[0] = downId;
+    }
+
+  int nbvol =0;
+  for (int i=0; i<nbFaces; i++)
+    {
+      int vtkTypeFace = cellTypes[i];
+      int downId = downCellId[i];
+      int nv = _downArray[vtkTypeFace]->getNumberOfUpCells(downId);
+      const int *upCells = _downArray[vtkTypeFace]->getUpCells(downId);
+      const unsigned char* upTypes = _downArray[vtkTypeFace]->getUpTypes(downId);
+       for (int j=0; j<nv; j++)
+        {
+          int vtkVolId = _downArray[upTypes[j]]->getVtkCellId(upCells[j]);
+          if (vtkVolId >= 0)
+            volVtkIds[nbvol++] = vtkVolId;
+        }
+    }
+  return nbvol;
+}
+
+/*! get the volumes containing a face or an edge of the downward structure
+ * The edge or face does not necessary belong to the vtkUnstructuredGrid
+ * @param volVtkIds vector of parent volume ids to fill (reserve enough space!)
+ * @param downId id in the downward structure
+ * @param downType type of cell
+ */
+int SMDS_UnstructuredGrid::GetParentVolumes(int* volVtkIds, int downId, unsigned char downType)
+{
+  int vtkType = downType;
+  int dim = SMDS_Downward::getCellDimension(vtkType);
+  int nbFaces = 0;
+  unsigned char cellTypes[1000];
+  int downCellId[1000];
+  if (dim == 1)
+    {
+      nbFaces = _downArray[vtkType]->getNumberOfUpCells(downId);
+      const int *upCells = _downArray[vtkType]->getUpCells(downId);
+      const unsigned char* upTypes = _downArray[vtkType]->getUpTypes(downId);
+      for (int i=0; i< nbFaces; i++)
+        {
+          cellTypes[i] = upTypes[i];
+          downCellId[i] = upCells[i];
+        }
+    }
+  else if (dim == 2)
+    {
+      nbFaces = 1;
+      cellTypes[0] = vtkType;
+      downCellId[0] = downId;
+    }
+
+  int nbvol =0;
+  for (int i=0; i<nbFaces; i++)
+    {
+      int vtkTypeFace = cellTypes[i];
+      int downId = downCellId[i];
+      int nv = _downArray[vtkTypeFace]->getNumberOfUpCells(downId);
+      const int *upCells = _downArray[vtkTypeFace]->getUpCells(downId);
+      const unsigned char* upTypes = _downArray[vtkTypeFace]->getUpTypes(downId);
+       for (int j=0; j<nv; j++)
+        {
+          int vtkVolId = _downArray[upTypes[j]]->getVtkCellId(upCells[j]);
+          if (vtkVolId >= 0)
+            volVtkIds[nbvol++] = vtkVolId;
+        }
+    }
+  return nbvol;
+}
+
 /*! get the node id's of a cell.
  * The cell is defined by it's downward connectivity id and type.
  * @param nodeSet set of of vtk node id's to fill.
@@ -779,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();
 }
 
@@ -805,3 +957,118 @@ void SMDS_UnstructuredGrid::BuildLinks()
   this->Links->BuildLinks(this, this->Connectivity);
   this->Links->Delete();
 }
+
+/*! Create a volume (prism or hexahedron) by duplication of a face.
+ * Designed for use in creation of flat elements separating volume domains.
+ * A face separating two domains is shared by two volume cells.
+ * All the nodes are already created (for the two faces).
+ * Each original Node is associated to corresponding nodes in the domains.
+ * Some nodes may be duplicated for more than two domains, when domain separations intersect.
+ * In that case, even some of the nodes to use for the original face may be changed.
+ * @param vtkVolId: vtk id of a volume containing the face, to get an orientation for the face.
+ * @param domain1: domain of the original face
+ * @param domain2: domain of the duplicated face
+ * @param originalNodes: the vtk node ids of the original face
+ * @param nodeDomains: map(original id --> map(domain --> duplicated node id))
+ * @return ok if success.
+ */
+SMDS_MeshCell* SMDS_UnstructuredGrid::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)
+{
+  //MESSAGE("extrudeVolumeFromFace " << vtkVolId);
+  vector<vtkIdType> orderedOriginals;
+  orderedOriginals.clear();
+  set<int>::const_iterator it = originalNodes.begin();
+  for (; it != originalNodes.end(); ++it)
+    orderedOriginals.push_back(*it);
+
+  int dim = 0;
+  int nbNodes = this->getOrderedNodesOfFace(vtkVolId, dim, orderedOriginals);
+  vector<vtkIdType> orderedNodes;
+
+  bool isQuadratic = false;
+  switch (orderedOriginals.size())
+  {
+    case 3:
+      if (dim == 2)
+        isQuadratic = true;
+      break;
+    case 6:
+    case 8:
+      isQuadratic = true;
+      break;
+    default:
+      isQuadratic = false;
+      break;
+  }
+
+  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 0;
+}