Salome HOME
Merge from V6_5_BR 05/06/2012
[modules/smesh.git] / src / SMDS / SMDS_UnstructuredGrid.cxx
index 8c394c8d7589e14174b8e10b296801d5a202a04d..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;
 
@@ -131,17 +152,9 @@ void SMDS_UnstructuredGrid::setSMDS_mesh(SMDS_Mesh *mesh)
 void SMDS_UnstructuredGrid::compactGrid(std::vector<int>& idNodesOldToNew, int newNodeSize,
                                         std::vector<int>& idCellsOldToNew, int newCellSize)
 {
-  // TODO utiliser mieux vtk pour faire plus simple (plus couteux ?)
-
   MESSAGE("------------------------- SMDS_UnstructuredGrid::compactGrid " << newNodeSize << " " << newCellSize);CHRONO(1);
   int alreadyCopied = 0;
 
-  //    if (this->Links)
-  //    {
-  //            this->Links->UnRegister(this);
-  //            this->Links = 0;
-  //    }
-
   // --- if newNodeSize, create a new compacted vtkPoints
 
   vtkPoints *newPoints = vtkPoints::New();
@@ -150,7 +163,7 @@ void SMDS_UnstructuredGrid::compactGrid(std::vector<int>& idNodesOldToNew, int n
   if (newNodeSize)
     {
       MESSAGE("-------------- compactGrid, newNodeSize " << newNodeSize);
-      // rnv: to fix bug "21125: EDF 1233 SMESH: Degrardation of precision in a test case for quadratic conversion"
+      // rnv: to fix bug "21125: EDF 1233 SMESH: Degradation of precision in a test case for quadratic conversion"
       // using double type for storing coordinates of nodes instead float.
       int oldNodeSize = idNodesOldToNew.size();
 
@@ -167,56 +180,6 @@ void SMDS_UnstructuredGrid::compactGrid(std::vector<int>& idNodesOldToNew, int n
         int endBloc = i;
         copyNodes(newPoints, idNodesOldToNew, alreadyCopied, startBloc, endBloc);
       }
-//       for (int i = 0; i < oldNodeSize; i++)
-//         {
-//           //MESSAGE("                                    " << i << " " << idNodesOldToNew[i]);
-//           switch (compactState)
-//           {
-//             case lookHoleStart:
-//               if (idNodesOldToNew[i] < 0)
-//                 {
-//                   MESSAGE("-------------- newNodeSize, startHole " << i << " " << oldNodeSize);
-//                   startHole = i;
-//                   if (!alreadyCopied) // copy the first bloc
-//                     {
-//                       MESSAGE("--------- copy first nodes before hole " << i << " " << oldNodeSize);
-//                       copyNodes(newPoints, idNodesOldToNew, alreadyCopied, 0, startHole);
-//                     }
-//                   compactState = lookHoleEnd;
-//                 }
-//               break;
-//             case lookHoleEnd:
-//               if (idNodesOldToNew[i] >= 0)
-//                 {
-//                   MESSAGE("-------------- newNodeSize, endHole " << i << " " << oldNodeSize);
-//                   endHole = i;
-//                   startBloc = i;
-//                   compactState = lookBlocEnd;
-//                 }
-//               break;
-//             case lookBlocEnd:
-//               if (idNodesOldToNew[i] < 0)
-//                 endBloc = i; // see nbPoints below
-//               else if (i == (oldNodeSize - 1))
-//                 endBloc = i + 1;
-//               if (endBloc)
-//                 {
-//                   MESSAGE("-------------- newNodeSize, endbloc " << endBloc << " " << oldNodeSize);
-//                   copyNodes(newPoints, idNodesOldToNew, alreadyCopied, startBloc, endBloc);
-//                   compactState = lookHoleEnd;
-//                   startHole = i;
-//                   endHole = 0;
-//                   startBloc = 0;
-//                   endBloc = 0;
-//                 }
-//               break;
-//           }
-//         }
-//       if (!alreadyCopied) // no hole, but shorter, no need to modify idNodesOldToNew
-//         {
-//           MESSAGE("------------- newNodeSize, shorter " << oldNodeSize);
-//           copyNodes(newPoints, idNodesOldToNew, alreadyCopied, 0, newNodeSize);
-//         }
       newPoints->Squeeze();
     }
 
@@ -258,60 +221,8 @@ void SMDS_UnstructuredGrid::compactGrid(std::vector<int>& idNodesOldToNew, int n
       copyBloc(newTypes, idCellsOldToNew, idNodesOldToNew, newConnectivity, newLocations, pointsCell,
                alreadyCopied, startBloc, endBloc);
   }
-//   for (int i = 0; i < oldCellSize; i++)
-//     {
-//       switch (compactState)
-//       {
-//         case lookHoleStart:
-//           if (this->Types->GetValue(i) == VTK_EMPTY_CELL)
-//             {
-//               MESSAGE(" -------- newCellSize, startHole " << i << " " << oldCellSize);
-//               startHole = i;
-//               compactState = lookHoleEnd;
-//               if (!alreadyCopied) // copy the first bloc
-//                 {
-//                   MESSAGE("--------- copy first bloc before hole " << i << " " << oldCellSize);
-//                   copyBloc(newTypes, idCellsOldToNew, idNodesOldToNew, newConnectivity, newLocations, pointsCell,
-//                            alreadyCopied, 0, startHole);
-//                 }
-//             }
-//           break;
-//         case lookHoleEnd:
-//           if (this->Types->GetValue(i) != VTK_EMPTY_CELL)
-//             {
-//               MESSAGE(" -------- newCellSize, EndHole " << i << " " << oldCellSize);
-//               endHole = i;
-//               startBloc = i;
-//               compactState = lookBlocEnd;
-//               holes += endHole - startHole;
-//             }
-//           break;
-//         case lookBlocEnd:
-//           endBloc = 0;
-//           if (this->Types->GetValue(i) == VTK_EMPTY_CELL)
-//             endBloc = i;
-//           else if (i == (oldCellSize - 1))
-//             endBloc = i + 1;
-//           if (endBloc)
-//             {
-//               MESSAGE(" -------- newCellSize, endBloc " << endBloc << " " << oldCellSize);
-//               copyBloc(newTypes, idCellsOldToNew, idNodesOldToNew, newConnectivity, newLocations, pointsCell,
-//                        alreadyCopied, startBloc, endBloc);
-//               compactState = lookHoleEnd;
-//             }
-//           break;
-//       }
-//     }
-//   if (!alreadyCopied) // no hole, but shorter
-//     {
-//       MESSAGE(" -------- newCellSize, shorter " << oldCellSize);
-//       copyBloc(newTypes, idCellsOldToNew, idNodesOldToNew, newConnectivity, newLocations, pointsCell, alreadyCopied, 0,
-//                oldCellSize);
-//     }
 
   newConnectivity->Squeeze();
-  //newTypes->Squeeze();
-  //newLocations->Squeeze();
 
   if (1/*newNodeSize*/)
     {
@@ -385,7 +296,6 @@ void SMDS_UnstructuredGrid::copyNodes(vtkPoints *newPoints, std::vector<int>& id
       memcpy(target, source, 3 * sizeof(double) * nbPoints);
       for (int j = start; j < end; j++)
         idNodesOldToNew[j] = alreadyCopied++; // old vtkId --> new vtkId
-        //idNodesOldToNew[alreadyCopied++] = idNodesOldToNew[j]; // new vtkId --> old SMDS id
     }
 }
 
@@ -421,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];
 }
 
@@ -431,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.
  *
@@ -442,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:
@@ -813,15 +743,15 @@ int SMDS_UnstructuredGrid::GetNeighbors(int* neighborsVtkIds, int* downIds, unsi
 {
   int vtkType = this->GetCellType(vtkId);
   int cellDim = SMDS_Downward::getCellDimension(vtkType);
-  if (cellDim != 3)
-    return 0; // TODO voisins des faces ou edges
+  if (cellDim <2)
+    return 0; // TODO voisins des edges = edges connectees
   int cellId = this->_cellIdToDownId[vtkId];
 
   int nbCells = _downArray[vtkType]->getNumberOfDownCells(cellId);
   const int *downCells = _downArray[vtkType]->getDownCells(cellId);
   const unsigned char* downTyp = _downArray[vtkType]->getDownTypes(cellId);
 
-  // --- iteration on faces of the 3D cell.
+  // --- iteration on faces of the 3D cell (or edges on the 2D cell).
 
   int nb = 0;
   for (int i = 0; i < nbCells; i++)
@@ -832,7 +762,8 @@ int SMDS_UnstructuredGrid::GetNeighbors(int* neighborsVtkIds, int* downIds, unsi
       const int *upCells = _downArray[cellType]->getUpCells(downId);
       const unsigned char* upTypes = _downArray[cellType]->getUpTypes(downId);
 
-      // --- max 2 upCells, one is this cell, the other is a neighbor
+      // ---for a volume, max 2 upCells, one is this cell, the other is a neighbor
+      //    for a face, number of neighbors (connected faces) not known
 
       for (int j = 0; j < nbUp; j++)
         {
@@ -850,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.
@@ -884,21 +924,22 @@ void SMDS_UnstructuredGrid::ModifyCellNodes(int vtkVolId, std::map<int, int> loc
     }
 }
 
-/*! Create a volume (prism or hexahedron) by duplication of a face.
- * the nodes of the new face are already created.
+/*! reorder the nodes of a face
  * @param vtkVolId vtk id of a volume containing the face, to get an orientation for the face.
- * @param localClonedNodeIds map old node id to new node id.
- * @return vtk id of the new volume.
+ * @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();
 }
 
@@ -916,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;
+}