Salome HOME
0021347: [CEA 497] Visualisation into SMESH and VISU of hexagonal prism cells (MED_OC...
[modules/smesh.git] / src / SMDS / SMDS_UnstructuredGrid.cxx
index bbe371c8dd0450e156b4e25aa71568f8a38a2535..8015105e7fe6a7aa85d9f721f8f2da7d9857913d 100644 (file)
@@ -1,8 +1,28 @@
+// Copyright (C) 2010-2011  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;
 
@@ -92,7 +113,7 @@ int SMDS_UnstructuredGrid::InsertNextLinkedCell(int type, int npts, vtkIdType *p
     return vtkUnstructuredGrid::InsertNextLinkedCell(type, npts, pts);
 
   // --- type = VTK_POLYHEDRON
-  MESSAGE("InsertNextLinkedCell VTK_POLYHEDRON");
+  //MESSAGE("InsertNextLinkedCell VTK_POLYHEDRON");
   int cellid = this->InsertNextCell(type, npts, pts);
 
   set<vtkIdType> setOfNodes;
@@ -105,7 +126,7 @@ int SMDS_UnstructuredGrid::InsertNextLinkedCell(int type, int npts, vtkIdType *p
       i++;
       for (int k = 0; k < nbnodes; k++)
         {
-          MESSAGE(" cell " << cellid << " face " << nf << " node " << pts[i]);
+          //MESSAGE(" cell " << cellid << " face " << nf << " node " << pts[i]);
           setOfNodes.insert(pts[i]);
           i++;
         }
@@ -114,7 +135,7 @@ int SMDS_UnstructuredGrid::InsertNextLinkedCell(int type, int npts, vtkIdType *p
   set<vtkIdType>::iterator it = setOfNodes.begin();
   for (; it != setOfNodes.end(); ++it)
     {
-      MESSAGE("reverse link for node " << *it << " cell " << cellid);
+      //MESSAGE("reverse link for node " << *it << " cell " << cellid);
       this->Links->ResizeCellList(*it, 1);
       this->Links->AddCellReference(cellid, *it);
     }
@@ -131,90 +152,34 @@ 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 startHole = 0;
-  int endHole = 0;
-  int startBloc = 0;
-  int endBloc = 0;
   int alreadyCopied = 0;
-  int holes = 0;
-
-  typedef enum
-  {
-    lookHoleStart, lookHoleEnd, lookBlocEnd
-  } enumState;
-  enumState compactState = lookHoleStart;
-
-  //   if (this->Links)
-  //   {
-  //           this->Links->UnRegister(this);
-  //           this->Links = 0;
-  //   }
 
   // --- if newNodeSize, create a new compacted vtkPoints
 
-  vtkPoints *newPoints = 0;
+  vtkPoints *newPoints = vtkPoints::New();
+  newPoints->SetDataType(VTK_DOUBLE);
+  newPoints->SetNumberOfPoints(newNodeSize);
   if (newNodeSize)
     {
       MESSAGE("-------------- compactGrid, newNodeSize " << newNodeSize);
-      newPoints = vtkPoints::New();
-      newPoints->Initialize();
-      newPoints->Allocate(newNodeSize);
-      newPoints->SetNumberOfPoints(newNodeSize);
+      // 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();
 
-      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);
-        }
+      int i = 0;
+      while ( i < oldNodeSize )
+      {
+        // skip a hole if any
+        while ( i < oldNodeSize && idNodesOldToNew[i] < 0 )
+          ++i;
+        int startBloc = i;
+        // look for a block end
+        while ( i < oldNodeSize && idNodesOldToNew[i] >= 0 )
+          ++i;
+        int endBloc = i;
+        copyNodes(newPoints, idNodesOldToNew, alreadyCopied, startBloc, endBloc);
+      }
       newPoints->Squeeze();
     }
 
@@ -236,99 +201,83 @@ void SMDS_UnstructuredGrid::compactGrid(std::vector<int>& idNodesOldToNew, int n
   newLocations->Initialize();
   newLocations->SetNumberOfValues(newCellSize);
 
-  startHole = 0;
-  endHole = 0;
-  startBloc = 0;
-  endBloc = 0;
-  alreadyCopied = 0;
-  holes = 0;
-  compactState = lookHoleStart;
-
   // TODO some polyhedron may be huge (only in some tests)
   vtkIdType tmpid[NBMAXNODESINCELL];
   vtkIdType *pointsCell = &tmpid[0]; // --- points id to fill a new cell
 
-  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);
-    }
+  alreadyCopied = 0;
+  int i = 0;
+  while ( i < oldCellSize )
+  {
+    // skip a hole if any
+    while ( i < oldCellSize && this->Types->GetValue(i) == VTK_EMPTY_CELL )
+      ++i;
+    int startBloc = i;
+    // look for a block end
+    while ( i < oldCellSize && this->Types->GetValue(i) != VTK_EMPTY_CELL )
+      ++i;
+    int endBloc = i;
+    if ( endBloc > startBloc )
+      copyBloc(newTypes, idCellsOldToNew, idNodesOldToNew, newConnectivity, newLocations, pointsCell,
+               alreadyCopied, startBloc, endBloc);
+  }
 
   newConnectivity->Squeeze();
-  //newTypes->Squeeze();
-  //newLocations->Squeeze();
 
-  if (newNodeSize)
+  if (1/*newNodeSize*/)
     {
       MESSAGE("------- newNodeSize, setPoints");
       this->SetPoints(newPoints);
       MESSAGE("NumberOfPoints: " << this->GetNumberOfPoints());
     }
-//#ifdef VTK_HAVE_POLYHEDRON
-  // TODO compact faces for Polyhedrons
-  // refaire completement faces et faceLocation
-  // pour chaque cell, recup oldCellId, oldFacesId, recopie dans newFaces de la faceStream
-  // en changeant les numeros de noeuds
-//  vtkIdTypeArray *newFaceLocations = vtkIdTypeArray::New();
-//  newFaceLocations->Initialize();
-//  vtkIdTypeArray *newFaces = vtkIdTypeArray::New();
-//  newFaces->Initialize();
-//  newFaceLocations->DeepCopy(this->FaceLocations);
-//  newFaces->DeepCopy(this->Faces);
-//  this->SetCells(newTypes, newLocations, newConnectivity, newFaceLocations, newFaces);
-//  newFaceLocations->Delete();
-//  newFaces->Delete();
-  if (this->FaceLocations) this->FaceLocations->Register(this);
-  if (this->Faces) this->Faces->Register(this);
-  this->SetCells(newTypes, newLocations, newConnectivity, FaceLocations, Faces);
-//#else
-//  this->SetCells(newTypes, newLocations, newConnectivity);
-//#endif
+
+  if (this->FaceLocations)
+    {
+      vtkIdTypeArray *newFaceLocations = vtkIdTypeArray::New();
+      newFaceLocations->Initialize();
+      newFaceLocations->Allocate(newTypes->GetSize());
+      vtkIdTypeArray *newFaces = vtkIdTypeArray::New();
+      newFaces->Initialize();
+      newFaces->Allocate(this->Faces->GetSize());
+      for (int i = 0; i < oldCellSize; i++)
+        {
+          if (this->Types->GetValue(i) == VTK_EMPTY_CELL)
+            continue;
+          int newCellId = idCellsOldToNew[i];
+          if (newTypes->GetValue(newCellId) == VTK_POLYHEDRON)
+            {
+               newFaceLocations->InsertNextValue(newFaces->GetMaxId()+1);
+               int oldFaceLoc = this->FaceLocations->GetValue(i);
+               int nCellFaces = this->Faces->GetValue(oldFaceLoc++);
+               newFaces->InsertNextValue(nCellFaces);
+               for (int n=0; n<nCellFaces; n++)
+                 {
+                   int nptsInFace = this->Faces->GetValue(oldFaceLoc++);
+                   newFaces->InsertNextValue(nptsInFace);
+                   for (int k=0; k<nptsInFace; k++)
+                     {
+                       int oldpt = this->Faces->GetValue(oldFaceLoc++);
+                       newFaces->InsertNextValue(idNodesOldToNew[oldpt]);
+                     }
+                 }
+            }
+          else
+            {
+               newFaceLocations->InsertNextValue(-1);
+            }
+        }
+      newFaceLocations->Squeeze();
+      newFaces->Squeeze();
+      newFaceLocations->Register(this);
+      newFaces->Register(this);
+      this->SetCells(newTypes, newLocations, newConnectivity, newFaceLocations, newFaces);
+      newFaceLocations->Delete();
+      newFaces->Delete();
+    }
+  else
+    this->SetCells(newTypes, newLocations, newConnectivity, FaceLocations, Faces);
+
+  newPoints->Delete();
   newTypes->Delete();
   newLocations->Delete();
   newConnectivity->Delete();
@@ -344,10 +293,9 @@ void SMDS_UnstructuredGrid::copyNodes(vtkPoints *newPoints, std::vector<int>& id
   int nbPoints = end - start;
   if (nbPoints > 0)
     {
-      memcpy(target, source, 3 * sizeof(float) * nbPoints);
+      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
     }
 }
 
@@ -373,7 +321,7 @@ void SMDS_UnstructuredGrid::copyBloc(vtkUnsignedCharArray *newTypes, std::vector
           pointsCell[l] = idNodesOldToNew[oldval];
           //MESSAGE("   " << oldval << " " << pointsCell[l]);
         }
-      int newcnt = newConnectivity->InsertNextCell(nbpts, pointsCell);
+      /*int newcnt = */newConnectivity->InsertNextCell(nbpts, pointsCell);
       int newLoc = newConnectivity->GetInsertLocation(nbpts);
       //MESSAGE(newcnt << " " << newLoc);
       newLocations->SetValue(alreadyCopied, newLoc);
@@ -383,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];
 }
 
@@ -393,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.
  *
@@ -404,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:
@@ -544,7 +512,7 @@ void SMDS_UnstructuredGrid::BuildDownwardConnectivity(bool withEdges)
           int vtkVolId = i;
           // MESSAGE("vtk volume " << vtkVolId);
           //ASSERT(_downArray[vtkType]);
-          int connVolId = _downArray[vtkType]->addCell(vtkVolId);
+          /*int connVolId = */_downArray[vtkType]->addCell(vtkVolId);
 
           // --- find all the faces of the volume, describe the faces by their nodes
 
@@ -775,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++)
@@ -794,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++)
         {
@@ -812,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.
@@ -846,13 +924,12 @@ 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<int>& orderedNodes)
+int SMDS_UnstructuredGrid::getOrderedNodesOfFace(int vtkVolId, std::vector<vtkIdType>& orderedNodes)
 {
   int vtkType = this->GetCellType(vtkVolId);
   int cellDim = SMDS_Downward::getCellDimension(vtkType);
@@ -878,3 +955,94 @@ 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_MeshVolume* 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 nbNodes = this->getOrderedNodesOfFace(vtkVolId, orderedOriginals);
+  vector<vtkIdType> orderedNodes;
+
+  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]);
+      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);
+          }
+      }
+      break;
+    default:
+      ASSERT(0);
+  }
+
+  SMDS_MeshVolume *vol = _mesh->AddVolumeFromVtkIds(orderedNodes);
+
+  // TODO update subshape list of elements and nodes
+  return vol;
+}