Salome HOME
remove a reference to the $MED_ROOT_DIR in the Makefile.in wich is useless
[modules/med.git] / src / MEDMEM / MEDMEM_Connectivity.cxx
index 9d57d54dbe1b85d967b8c8ca298ca3040a35f0ae..2c8917e540936506e9642a688b2663620412ed5d 100644 (file)
@@ -1,4 +1,3 @@
-using namespace std;
 #include "MEDMEM_Connectivity.hxx"
 #include "MEDMEM_Family.hxx"
 #include "MEDMEM_Group.hxx"
@@ -10,7 +9,22 @@ using namespace std;
 #include "MEDMEM_STRING.hxx"
 #include <iomanip>
 
+using namespace std;
 using namespace MEDMEM;
+using namespace MED_EN;
+
+// Enlarge the vector if necessary to insert the element
+static inline void insert_vector(vector<int> &Vect, int Indice, int Element)
+{
+  if (Indice >= Vect.capacity())
+    Vect.reserve(Indice + 1000);
+
+  if (Indice >= Vect.size())
+    Vect.resize(Indice+1);
+
+  Vect[Indice] = Element;
+}
+
 /*!
    Default Constructor. /n
    Default for Entity is MED_CELL and type of Connectivity is MED_NODAL */
@@ -23,9 +37,14 @@ CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity /* =MED_CELL */) :
                                _geometricTypes((medGeometryElement*)NULL),
                                _type((CELLMODEL*)NULL),
                                _entityDimension(0),
+                               _numberOfNodes(0),
                                _count((int*)NULL),
                                _nodal((MEDSKYLINEARRAY*)NULL),
+                               _polygonsNodal((MEDSKYLINEARRAY*)NULL),
+                               _polyhedronNodal((POLYHEDRONARRAY*)NULL),
                                _descending((MEDSKYLINEARRAY*)NULL),
+                               _polygonsDescending((MEDSKYLINEARRAY*)NULL),
+                               _polyhedronDescending((MEDSKYLINEARRAY*)NULL),
                                _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
                                _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
                                _neighbourhood((MEDSKYLINEARRAY*)NULL),
@@ -44,8 +63,13 @@ CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity /* =MED_CELL *
                                _typeConnectivity(MED_NODAL),
                                _numberOfTypes(numberOfTypes),
                                _entityDimension(0),
+                               _numberOfNodes(0),
                                _nodal((MEDSKYLINEARRAY*)NULL),
+                               _polygonsNodal((MEDSKYLINEARRAY*)NULL),
+                               _polyhedronNodal((POLYHEDRONARRAY*)NULL),
                                _descending((MEDSKYLINEARRAY*)NULL),
+                               _polygonsDescending((MEDSKYLINEARRAY*)NULL),
+                               _polyhedronDescending((MEDSKYLINEARRAY*)NULL),
                                _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
                                _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
                                _neighbourhood((MEDSKYLINEARRAY*)NULL),
@@ -77,7 +101,7 @@ CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
  else
     _geometricTypes = (medGeometryElement *) NULL;
 
- if (m._type != NULL) 
+ if (m._type != NULL)
  {
     _type = new CELLMODEL[_numberOfTypes];
     for (int i=0;i<_numberOfTypes;i++)
@@ -88,22 +112,42 @@ CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
 
  if (m._count != NULL)
  {
-      _count = new med_int[_numberOfTypes+1];
-      memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(med_int));
+      _count = new int[_numberOfTypes+1];
+      memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(int));
  }
  else
-    _count = (med_int *) NULL;
+    _count = (int *) NULL;
 
  if (m._nodal != NULL)
     _nodal = new MEDSKYLINEARRAY(* m._nodal);
  else
     _nodal = (MEDSKYLINEARRAY *) NULL;
 
+ if (m._polygonsNodal != NULL)
+   _polygonsNodal = new MEDSKYLINEARRAY(* m._polygonsNodal);
+ else
+   _polygonsNodal = (MEDSKYLINEARRAY *) NULL;
+
+ if (m._polyhedronNodal != NULL)
+   _polyhedronNodal = new POLYHEDRONARRAY(* m._polyhedronNodal);
+ else
+   _polyhedronNodal = (POLYHEDRONARRAY *) NULL;
+
  if (m._descending != NULL)
     _descending = new MEDSKYLINEARRAY(* m._descending);
  else
     _descending = (MEDSKYLINEARRAY *) NULL;
 
+ if (m._polygonsDescending != NULL)
+   _polygonsDescending = new MEDSKYLINEARRAY(* m._polygonsDescending);
+ else
+   _polygonsDescending = (MEDSKYLINEARRAY *) NULL;
+
+ if (m._polyhedronDescending != NULL)
+   _polyhedronDescending = new MEDSKYLINEARRAY(* m._polyhedronDescending);
+ else
+   _polyhedronDescending = (MEDSKYLINEARRAY *) NULL;
+
  if (m._reverseNodalConnectivity != NULL)
     _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
  else
@@ -142,8 +186,16 @@ CONNECTIVITY::~CONNECTIVITY()
      delete [] _count;
   if (_nodal != NULL)
      delete _nodal;
+  if (_polygonsNodal != NULL)
+     delete _polygonsNodal;
+  if (_polyhedronNodal != NULL)
+     delete _polyhedronNodal;
   if (_descending != NULL)
      delete _descending;
+  if (_polygonsDescending != NULL)
+     delete _polygonsDescending;
+  if (_polyhedronDescending != NULL)
+    delete _polyhedronDescending;
   if (_reverseNodalConnectivity != NULL)
      delete _reverseNodalConnectivity;
   if (_reverseDescendingConnectivity != NULL)
@@ -167,7 +219,7 @@ void CONNECTIVITY::setConstituent(CONNECTIVITY * Constituent)
   if (Entity == MED_CELL)
     throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : could not set constituent on MED_CELL !"));
 
-  if ((Entity == MED_EDGE)&(_entityDimension == 3)) 
+  if ((Entity == MED_EDGE)&(_entityDimension == 3))
   {
     if (_constituent == NULL)
       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : Entity not found !"));
@@ -185,13 +237,13 @@ void CONNECTIVITY::setGeometricTypes(const medGeometryElement * Types,
 //--------------------------------------------------------------------//
 {
   if (Entity == _entity)
-    for (int i=0; i<_numberOfTypes; i++) 
+    for (int i=0; i<_numberOfTypes; i++)
     {
       _geometricTypes[i] = Types[i];
       _type[i] = CELLMODEL(_geometricTypes[i]);
     }
-  else 
-    { 
+  else
+    {
     if (_constituent == NULL)
         throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setGeometricTypes : Entity not found !"));
     _constituent->setGeometricTypes(Types,Entity);
@@ -204,7 +256,7 @@ void CONNECTIVITY::setCount(const int * Count, const medEntityMesh Entity)
                            throw (MEDEXCEPTION)
 //--------------------------------------------------------------------//
 {
-  if (Entity == _entity) 
+  if (Entity == _entity)
   {
     int * index = new int[Count[_numberOfTypes]];
     index[0]=1;
@@ -236,7 +288,7 @@ void CONNECTIVITY::setNodal(const int * Connectivity,
                            throw (MEDEXCEPTION)
 //--------------------------------------------------------------------//
 {
-  if (_entity == Entity) 
+  if (_entity == Entity)
   {
     // find geometric type
     bool find = false;
@@ -250,7 +302,7 @@ void CONNECTIVITY::setNodal(const int * Connectivity,
       }
     if (!find)
       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : geometric type not found"));
-  } else 
+  } else
   {
     if (_constituent == NULL)
       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : Entity not found !"));
@@ -258,6 +310,74 @@ void CONNECTIVITY::setNodal(const int * Connectivity,
   }
 }
 
+
+//--------------------------------------------------------------------//
+void CONNECTIVITY::setPolygonsConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, const int* PolygonsConnectivity, const int* PolygonsConnectivityIndex, int ConnectivitySize, int NumberOfPolygons)
+//--------------------------------------------------------------------//
+{
+  const char* LOC = "CONNECTIVITY::setPolygonsConnectivity";
+  BEGIN_OF(LOC);
+
+  if (_entity == Entity)
+    {
+      MEDSKYLINEARRAY* Connectivity = new MEDSKYLINEARRAY(NumberOfPolygons,ConnectivitySize,PolygonsConnectivityIndex,PolygonsConnectivity);
+
+      if (ConnectivityType == MED_NODAL)
+       {
+         if (_polygonsNodal != (MEDSKYLINEARRAY*) NULL)
+           delete _polygonsNodal;
+         _polygonsNodal = Connectivity;
+       }
+      else
+       {
+         if (_typeConnectivity != MED_DESCENDING)
+           _typeConnectivity = MED_DESCENDING; //by default it is set to MED_NODAL
+         if (_polygonsDescending != (MEDSKYLINEARRAY*) NULL)
+           delete _polygonsDescending;
+         _polygonsDescending = Connectivity;
+       }
+    }
+  else
+    {
+      if (_constituent == (CONNECTIVITY*) NULL)
+       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not found !"));
+      _constituent->setPolygonsConnectivity(ConnectivityType, Entity, PolygonsConnectivity, PolygonsConnectivityIndex, ConnectivitySize, NumberOfPolygons);
+    }
+}
+
+
+//--------------------------------------------------------------------//
+void CONNECTIVITY::setPolyhedronConnectivity(medConnectivity ConnectivityType, const int* PolyhedronConnectivity, const int* PolyhedronIndex, int ConnectivitySize, int NumberOfPolyhedron, const int* PolyhedronFacesIndex /* =(int*)NULL */, int NumberOfFaces /* =0 */)
+//--------------------------------------------------------------------//
+{
+  const char* LOC = "CONNECTIVITY::setPolyhedronConnectivity";
+  BEGIN_OF(LOC);
+
+  if (_entity == MED_CELL)
+    {
+      if (ConnectivityType == MED_NODAL)
+       {
+         if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
+           delete _polyhedronNodal;
+         _polyhedronNodal = new POLYHEDRONARRAY(NumberOfPolyhedron,NumberOfFaces,ConnectivitySize);
+         _polyhedronNodal->setPolyhedronIndex(PolyhedronIndex);
+         _polyhedronNodal->setFacesIndex(PolyhedronFacesIndex);
+         _polyhedronNodal->setNodes(PolyhedronConnectivity);
+       }
+      else
+       {
+         if (_typeConnectivity != MED_DESCENDING)
+           _typeConnectivity = MED_DESCENDING; //by default it is set to MED_NODAL
+         if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
+           delete _polyhedronDescending;
+         _polyhedronDescending = new MEDSKYLINEARRAY(NumberOfPolyhedron,ConnectivitySize,PolyhedronIndex,PolyhedronConnectivity);
+       }
+    }
+  else
+    throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : _entity must be MED_CELL to set polyhedron !"));
+}
+
+
 /*! A DOCUMENTER */
 //------------------------------------------------------------------------------------------//
 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
@@ -268,11 +388,11 @@ void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEn
   // a temporary limitation due to calculteDescendingConnectivity function !!!
   if ((_entityDimension==3) & (Entity==MED_EDGE))
     throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
-  
+
   if (Entity==_entity)
     if (ConnectivityType==MED_NODAL)
       calculateNodalConnectivity();
-    else 
+    else
       if (Entity==MED_CELL)
        calculateDescendingConnectivity();
       else
@@ -287,526 +407,161 @@ void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEn
 /*!  Give, in full or no interlace mode (for nodal connectivity),
      descending or nodal connectivity.
 
-      You must give a <medEntityMesh> (ie:MED_EDGE) 
-      and a <medGeometryElement> (ie:MED_SEG3).  */
+      You must give a %medEntityMesh (ie:MED_EDGE)
+      and a %medGeometryElement (ie:MED_SEG3).  */
 
 //------------------------------------------------------------//
-void CONNECTIVITY::updateFamily(vector<FAMILY*> myFamilies) 
+void CONNECTIVITY::updateFamily(const vector<FAMILY*>& myFamilies)
 //------------------------------------------------------------//
 {
   const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
-  BEGIN_OF(LOC);
-
   int numberOfFamilies = myFamilies.size();
-  if (numberOfFamilies == 0 ) {
-    MESSAGE(LOC<<"No family");
+  if (numberOfFamilies == 0 || _constituent == NULL)
     return;
-  }
   // does we do an update ?
-  if ((_constituent != NULL)&(_descending != NULL)) {
-    MESSAGE(LOC<<"Constituent is already defined");
+  if ((_constituent != NULL) && (_descending != NULL))
     return;
-  }
-
-  if ((_constituent != NULL)&(_descending == NULL)) {
-    if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
-      MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
-      return;
+  if (myFamilies[0]->getEntity() != _constituent->getEntity())
+    return;
+  CONNECTIVITY * oldConstituent = _constituent;
+  _constituent = (CONNECTIVITY *)NULL;
+  if (oldConstituent->_nodal==NULL)
+    throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
+
+  //Loc vars defined to treat polygons exactly the same as classic types. Not nice but necessary.
+  int oldNumberOfFaceTab[2];
+  const int * oldConstituentValueTab[2];
+  const int * oldConstituentIndexTab[2];
+  int * renumberingFromOldToNewTab[2];//Final mapping array between old numbers and new numbers.;
+
+  int oldNumberOfFace = oldConstituent->_nodal->getNumberOf(); oldNumberOfFaceTab[0]=oldNumberOfFace;
+  const int * oldConstituentValue = oldConstituent->_nodal->getValue(); oldConstituentValueTab[0]=oldConstituentValue;
+  const int * oldConstituentIndex = oldConstituent->_nodal->getIndex(); oldConstituentIndexTab[0]=oldConstituentIndex;
+  int * renumberingFromOldToNew= new int [oldNumberOfFace]; renumberingFromOldToNewTab[0]=renumberingFromOldToNew;
+
+  int oldNumberOfFacePoly = oldConstituent->getNumberOfPolygons();
+  const int * oldConstituentValuePoly=0;
+  const int * oldConstituentIndexPoly=0;
+  int * renumberingFromOldToNewPoly=0;
+
+  int nbOfTurnInGlobalLoop=1;//Defined to know if a second search on polygons is needed.
+  if(oldNumberOfFacePoly>0)
+    {
+      oldNumberOfFaceTab[1]=oldNumberOfFacePoly;
+      oldConstituentValuePoly = oldConstituent->_polygonsNodal->getValue(); oldConstituentValueTab[1]=oldConstituentValuePoly;
+      oldConstituentIndexPoly = oldConstituent->_polygonsNodal->getIndex(); oldConstituentIndexTab[1]=oldConstituentIndexPoly;
+      renumberingFromOldToNewPoly=new int[oldNumberOfFacePoly]; renumberingFromOldToNewTab[1]=renumberingFromOldToNewPoly;
+      nbOfTurnInGlobalLoop++;
     }
 
-    // well we could go !
-    CONNECTIVITY * oldConstituent = _constituent;
-
-//      for(int i=0; i<numberOfFamilies; i++) {
-//        FAMILY * myFamily = myFamilies[i];
-//        MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
-//      }
-
-    _constituent = (CONNECTIVITY *)NULL;
-    // for instance we must have nodal connectivity in constituent :
-    if (oldConstituent->_nodal == NULL)
-      {
-       MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
-       throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
-      }
-    int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
-    const int * oldConstituentValue = oldConstituent->_nodal->getValue();
-    const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
-
-    calculateDescendingConnectivity();
+  calculateDescendingConnectivity();//perform calculation of descending connectivity to update all connectivities after having taken into account of direction of faces...
+  _constituent->calculateReverseNodalConnectivity();//getting d-1 nodal connectivity to find get new face numbers from nodes numbers...
 
-    //    int newNumberOfFace = _constituent->_nodal->getNumberOf();
-    const int * newConstituentValue = _constituent->_nodal->getValue();
-    const int * newConstituentIndex = _constituent->_nodal->getIndex();
-    
-    const int * newReverseDescendingIndex =
-      _reverseDescendingConnectivity->getIndex();
-    
-    const int * newDescendingIndex = _descending->getIndex();
-    //    const int * newDescendingValue = _descending->getValue();
-
-    // loop on all family,
-    //   for all constituent in family, we get it's old connectivity 
-    //   (with oldCconstituentValue and oldConstituentIndex)
-    //   and search the constituent in newConstituentValue with class 
-    //   ModulusArry
-    //
-    //   when a new face is found, replace old constituent 
-    //   number in family with new one
-    //   If face is not rigth oriented, we must change _descending attribute 
-    //   and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
-
-    // Voila a toi de jouer Nadir :-)
-
-    // First we get the renumbering from the oldCconstituentValue and
-    // oldConstituentIndex in the the new one, newConstituentValue and
-    // newConstituentIndex with a negative sign if the face is not
-    // right orented
-
-    int * renumberingFromOldToNew = new int [oldNumberOfFace];
-    int index1 = 0;
-    int indexm1 = 0;
-
-    _constituent->calculateReverseNodalConnectivity();
-    
-    for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
-      {
-       int index = 0;
-
-       renumberingFromOldToNew[iOldFace] = iOldFace+1;
-       //        renumberingFromOldToNew[iOldFace] = 999999;
-       
-       int face_it_beginOld = oldConstituentIndex[iOldFace];
-       int face_it_endOld = oldConstituentIndex[iOldFace+1];
-       int face_size_itOld = face_it_endOld - face_it_beginOld;
-
-       const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
-       int face_size_itNew;
-       
-       const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
-       const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
-
-       // set an array wich contains faces numbers arround first node 
-       int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
-       int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
-       int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
-
-       int * FacesList = new int[NumberOfFacesInList];
-
-       for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
-         FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
-       }
-       // foreach node in sub cell, we search elements which are in common
-       // at the end, we must have only one !
+  const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity(); //Common to polygons and classic geometric types
+  const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex(); //Common to polygons and classic geometric types
 
-       for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
-         {
-           int NewNumberOfFacesInList = 0;
-           int * NewFacesList = new int[NumberOfFacesInList];
-
-           for (int l1=0; l1<NumberOfFacesInList; l1++) {
-             for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
-               if (FacesList[l1]<reverseFaceNodal[l2-1])
-                 // increasing order : FacesList[l1] are not in elements list of node l
-                 break;
-               if (FacesList[l1]==reverseFaceNodal[l2-1]) {
-                 // we have found one
-                 NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
-                 NewNumberOfFacesInList++;
-                 break;
+  for(int loop=0;loop<nbOfTurnInGlobalLoop;loop++)
+    {
+      int oldNumberOfFaceLoop=oldNumberOfFaceTab[loop];
+      const int * oldConstituentValueLoop=oldConstituentValueTab[loop];
+      const int * oldConstituentIndexLoop= oldConstituentIndexTab[loop];
+      int * renumberingFromOldToNewLoop=renumberingFromOldToNewTab[loop];
+      for(int iOldFace=0;iOldFace<oldNumberOfFaceLoop;iOldFace++)
+       {
+         const int *nodesOfCurrentFaceOld=oldConstituentValueLoop+oldConstituentIndexLoop[iOldFace]-1;
+         int nbOfNodesOfCurrentFaceOld=oldConstituentIndexLoop[iOldFace+1]-oldConstituentIndexLoop[iOldFace];
+
+         //retrieving new number of polygon face...
+         int sizeOfNewFaceNb=reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1];
+         int *newFaceNb=new int[sizeOfNewFaceNb];
+         memcpy(newFaceNb,reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1]-1,sizeOfNewFaceNb*sizeof(int));//all faces candidates are in newFaceNb !!!
+         for(int curNode=1;curNode<nbOfNodesOfCurrentFaceOld && sizeOfNewFaceNb>1;curNode++)
+           {
+             const int *newFacesNbForCurNode=reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1]-1;
+             int sizeOfNewFacesNbForCurNode=reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1];
+             for(int i=0;i<sizeOfNewFaceNb;)
+               {
+                 bool found=false;
+                 for(int j=0;j<sizeOfNewFacesNbForCurNode && !found;j++)
+                   {
+                     if(newFacesNbForCurNode[j]==newFaceNb[i])
+                       found=true;
+                   }
+                 if(found)
+                   i++;
+                 else
+                   newFaceNb[i]=newFaceNb[--sizeOfNewFaceNb];
                }
-             }
            }
-           NumberOfFacesInList = NewNumberOfFacesInList;
-           delete [] FacesList;
-           FacesList = NewFacesList;
-         }
-
-       if (!NumberOfFacesInList==0)
-         {
-           if (NumberOfFacesInList>1)
-             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
-       
-           MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
-
-           int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
-           int face_it_endNew = newConstituentIndex[FacesList[0]];
-           face_size_itNew = face_it_endNew - face_it_beginNew;
-
-           const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
-           MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
-       
-           int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
-
-           //SCRUTE(retCompareNewOld);
-
-           // Real new face found
-       
-           if(retCompareNewOld == 1)
-             {
-               renumberingFromOldToNew[iOldFace] = FacesList[0];
-               index = 1;
-               index1++;
-             }
-       
-           // Reverse new face found
-       
-           if(retCompareNewOld == -1)
-             {
-               renumberingFromOldToNew[iOldFace] = FacesList[0];
-               index = 1;
-               indexm1++;
-           
-               int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
-               int face_it_end = newReverseDescendingIndex[FacesList[0]];
-               int face_size_it = face_it_end - face_it_begin;
-
-               if (face_size_it == 1)
-                 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
-           
-               if (face_size_it > 2)
-                 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
-           
-               // we have always 2 neighbourings
-               int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
-               int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
-               // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
-               //                  if (cell2 == 0)
-               //                    throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
-       
-               if (cell2 != 0) { // we are not on border !!!!
-             
-                 _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
-                 // Updating _constituent->_nodal because of reversity
-                 const int * oldArray = oldConstituentValue+face_it_beginOld-1;
-                 for(int iarray=1;iarray<=face_size_itNew;iarray++){
-                   _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
-                 }
-             
-                 // Updating _reverseDescendingConnectivity
-             
-             
-                 _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
-             
-                 // Updating _descending for cell1 and cell2
-                 for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
-                   if (_descending->getIndexValue(iface)==FacesList[0])
-                     _descending->setIndexValue(iface,-FacesList[0]);
-                   else if (_descending->getIndexValue(iface)==-FacesList[0])
-                     _descending->setIndexValue(iface,FacesList[0]);
-                 // else nothing to do
-             
-                 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
-                   if (_descending->getIndexValue(iface)==FacesList[0])
-                     _descending->setIndexValue(iface,-FacesList[0]);
-                   else if (_descending->getIndexValue(iface)==-FacesList[0])
-                     _descending->setIndexValue(iface,FacesList[0]);
-                 // else nothing to do
-
-               } else {// else we are on border and we do nothing !!!!!!!!
-                 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
-                 MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
-                 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
-               }
-             }
-       
-           if(index == 0)
-             {
-               MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
-               throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
-             }
-         }
-       delete [] FacesList;
-      }
-    
-    MESSAGE(LOC<<"The Renumbering is finished and the status is");
-
-    // Updating the Family
-    
-    for(int i=0; i<numberOfFamilies; i++) {
-      FAMILY * myFamily = myFamilies[i];
-      
-      MEDSKYLINEARRAY * number = myFamily->getnumber();
-      int numberOfLines_skyline = number->getNumberOf();
-      const int * index_skyline = number->getIndex();
-      
-      for (int i=0;i<numberOfLines_skyline;i++) {
-       for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
-         number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
+         if(sizeOfNewFaceNb!=1)
+           throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"unexisting face specified!!!"));
+         renumberingFromOldToNewLoop[iOldFace]=newFaceNb[0];
+         delete [] newFaceNb;
+         //end of retrieving new number of polygon face...
+         int nbOfNodesOfCurrentFaceNew;
+         const int *nodesOfCurrentFaceNew=_constituent->getConnectivityOfAnElementWithPoly(MED_NODAL,_constituent->getEntity(),
+                                                                                           renumberingFromOldToNewLoop[iOldFace],nbOfNodesOfCurrentFaceNew);
+         MEDMODULUSARRAY modulusArrayOld(nbOfNodesOfCurrentFaceOld,nodesOfCurrentFaceOld);
+         MEDMODULUSARRAY modulusArrayNew(nbOfNodesOfCurrentFaceNew,nodesOfCurrentFaceNew);
+         int retCompareNewOld=modulusArrayNew.compare(modulusArrayOld);
+         if(retCompareNewOld==0)
+           throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"Uncompatible given user face with calculated existing faces"));
+         if(retCompareNewOld==-1)
+           invertConnectivityForAFace(renumberingFromOldToNewLoop[iOldFace],nodesOfCurrentFaceOld,loop==1);
        }
-      }
-      MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
     }
-
-    delete oldConstituent ;
-    delete [] renumberingFromOldToNew;
-  }
-  
-
-  END_OF(LOC);
+  // Updating the Family
+  for(vector<FAMILY*>::const_iterator iter=myFamilies.begin();iter!=myFamilies.end();iter++)
+    (*iter)->changeElementsNbs(_constituent->getEntity(),renumberingFromOldToNew,oldNumberOfFace,renumberingFromOldToNewPoly);
+  delete oldConstituent ;
+  delete [] renumberingFromOldToNew;
+  if(oldNumberOfFacePoly>0)
+    delete [] renumberingFromOldToNewPoly;
   return;
 }
 
-// meme methode que updateFamily, mais avec des groupes. Il n'est pas possible d'utiliser
-// l'heritage car les pointeurs sont dans un conteneur.
-void CONNECTIVITY::updateGroup(vector<GROUP*> myFamilies)
-//------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------//
+const int * MEDMEM::CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
+//------------------------------------------------------------------------------------------------------------------//
 {
-  const char * LOC = "CONNECTIVITY::updateFamily(vector<GROUP*>) ";
+  const char * LOC = "CONNECTIVITY::getConnectivity";
   BEGIN_OF(LOC);
 
-  int numberOfFamilies = myFamilies.size();
-  if (numberOfFamilies == 0 ) {
-    MESSAGE(LOC<<"No family");
-    return;
-  }
-  // does we do an update ?
-  if ((_constituent != NULL)&(_descending != NULL)) {
-    MESSAGE(LOC<<"Constituent is already defined");
-    return;
-  }
-
-  if ((_constituent != NULL)&(_descending == NULL)) {
-    if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
-      MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
-      return;
-    }
-
-    // well we could go !
-    CONNECTIVITY * oldConstituent = _constituent;
-
-//      for(int i=0; i<numberOfFamilies; i++) {
-//        FAMILY * myFamily = myFamilies[i];
-//        MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
-//      }
+  MEDSKYLINEARRAY * Connectivity;
+  if (Entity==_entity) {
 
-    _constituent = (CONNECTIVITY *)NULL;
-    // for instance we must have nodal connectivity in constituent :
-    if (oldConstituent->_nodal == NULL)
+    if (ConnectivityType==MED_NODAL)
       {
-       MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
-       throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
+       calculateNodalConnectivity();
+       Connectivity=_nodal;
       }
-    int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
-    const int * oldConstituentValue = oldConstituent->_nodal->getValue();
-    const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
-
-    calculateDescendingConnectivity();
-
-    //    int newNumberOfFace = _constituent->_nodal->getNumberOf();
-    const int * newConstituentValue = _constituent->_nodal->getValue();
-    const int * newConstituentIndex = _constituent->_nodal->getIndex();
-    
-    const int * newReverseDescendingIndex =
-      _reverseDescendingConnectivity->getIndex();
-    
-    const int * newDescendingIndex = _descending->getIndex();
-    //    const int * newDescendingValue = _descending->getValue();
-
-    // loop on all family,
-    //   for all constituent in family, we get it's old connectivity 
-    //   (with oldCconstituentValue and oldConstituentIndex)
-    //   and search the constituent in newConstituentValue with class 
-    //   ModulusArry
-    //
-    //   when a new face is found, replace old constituent 
-    //   number in family with new one
-    //   If face is not rigth oriented, we must change _descending attribute 
-    //   and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
-
-    // Voila a toi de jouer Nadir :-)
-
-    // First we get the renumbering from the oldCconstituentValue and
-    // oldConstituentIndex in the the new one, newConstituentValue and
-    // newConstituentIndex with a negative sign if the face is not
-    // right orented
-
-    int * renumberingFromOldToNew = new int [oldNumberOfFace];
-    int index1 = 0;
-    int indexm1 = 0;
-
-    _constituent->calculateReverseNodalConnectivity();
-    
-    for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
+    else
       {
-       int index = 0;
-
-       renumberingFromOldToNew[iOldFace] = iOldFace+1;
-       //        renumberingFromOldToNew[iOldFace] = 999999;
-       
-       int face_it_beginOld = oldConstituentIndex[iOldFace];
-       int face_it_endOld = oldConstituentIndex[iOldFace+1];
-       int face_size_itOld = face_it_endOld - face_it_beginOld;
-
-       const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
-       int face_size_itNew;
-       
-       const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
-       const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
-
-       // set an array wich contains faces numbers arround first node 
-       int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
-       int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
-       int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
-
-       int * FacesList = new int[NumberOfFacesInList];
-
-       for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
-         FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
-       }
-       // foreach node in sub cell, we search elements which are in common
-       // at the end, we must have only one !
-
-       for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
-         {
-           int NewNumberOfFacesInList = 0;
-           int * NewFacesList = new int[NumberOfFacesInList];
-
-           for (int l1=0; l1<NumberOfFacesInList; l1++) {
-             for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
-               if (FacesList[l1]<reverseFaceNodal[l2-1])
-                 // increasing order : FacesList[l1] are not in elements list of node l
-                 break;
-               if (FacesList[l1]==reverseFaceNodal[l2-1]) {
-                 // we have found one
-                 NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
-                 NewNumberOfFacesInList++;
-                 break;
-               }
-             }
-           }
-           NumberOfFacesInList = NewNumberOfFacesInList;
-           delete [] FacesList;
-           FacesList = NewFacesList;
-         }
-
-       if (!NumberOfFacesInList==0)
-         {
-           if (NumberOfFacesInList>1)
-             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
-       
-           MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
-
-           int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
-           int face_it_endNew = newConstituentIndex[FacesList[0]];
-           face_size_itNew = face_it_endNew - face_it_beginNew;
-
-           const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
-           MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
-       
-           int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
-
-           //SCRUTE(retCompareNewOld);
-
-           // Real new face found
-       
-           if(retCompareNewOld == 1)
-             {
-               renumberingFromOldToNew[iOldFace] = FacesList[0];
-               index = 1;
-               index1++;
-             }
-       
-           // Reverse new face found
-       
-           if(retCompareNewOld == -1)
-             {
-               renumberingFromOldToNew[iOldFace] = FacesList[0];
-               index = 1;
-               indexm1++;
-           
-               int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
-               int face_it_end = newReverseDescendingIndex[FacesList[0]];
-               int face_size_it = face_it_end - face_it_begin;
-
-               if (face_size_it == 1)
-                 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
-           
-               if (face_size_it > 2)
-                 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
-           
-               // we have always 2 neighbourings
-               int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
-               int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
-               // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
-               //                  if (cell2 == 0)
-               //                    throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
-       
-               if (cell2 != 0) { // we are not on border !!!!
-             
-                 _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
-                 // Updating _constituent->_nodal because of reversity
-                 const int * oldArray = oldConstituentValue+face_it_beginOld-1;
-                 for(int iarray=1;iarray<=face_size_itNew;iarray++){
-                   _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
-                 }
-             
-                 // Updating _reverseDescendingConnectivity
-             
-             
-                 _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
-             
-                 // Updating _descending for cell1 and cell2
-                 for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
-                   if (_descending->getIndexValue(iface)==FacesList[0])
-                     _descending->setIndexValue(iface,-FacesList[0]);
-                   else if (_descending->getIndexValue(iface)==-FacesList[0])
-                     _descending->setIndexValue(iface,FacesList[0]);
-                 // else nothing to do
-             
-                 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
-                   if (_descending->getIndexValue(iface)==FacesList[0])
-                     _descending->setIndexValue(iface,-FacesList[0]);
-                   else if (_descending->getIndexValue(iface)==-FacesList[0])
-                     _descending->setIndexValue(iface,FacesList[0]);
-                 // else nothing to do
-
-               } else {// else we are on border and we do nothing !!!!!!!!
-                 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
-                 MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
-                 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
-               }
-             }
-       
-           if(index == 0)
-             {
-               MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
-               throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
-             }
-         }
-       delete [] FacesList;
-      }
-    
-    MESSAGE(LOC<<"The Renumbering is finished and the status is");
-
-    // Updating the Family
-    
-    for(int i=0; i<numberOfFamilies; i++) {
-      GROUP * myFamily = myFamilies[i];
-      
-      MEDSKYLINEARRAY * number = myFamily->getnumber();
-      int numberOfLines_skyline = number->getNumberOf();
-      const int * index_skyline = number->getIndex();
-      
-      for (int i=0;i<numberOfLines_skyline;i++) {
-       for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
-         number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
-       }
+       calculateDescendingConnectivity();
+       Connectivity=_descending;
       }
-      MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
-    }
 
-    delete oldConstituent ;
-    delete [] renumberingFromOldToNew;
-  }
-  
+    if (Connectivity!=NULL)
+      if (Type==MED_ALL_ELEMENTS)
+       return Connectivity->getValue();
+      else {
+       for (int i=0; i<_numberOfTypes; i++)
+         if (_geometricTypes[i]==Type)
+           //return Connectivity->getI(i+1);
+           return Connectivity->getI(_count[i]);
+       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
+      }
+    else
+      throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
+  } else
+    if (_constituent != NULL)
+      return _constituent->getConnectivity(ConnectivityType,Entity,Type);
 
-  END_OF(LOC);
-  return;
+  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 }
 
 //------------------------------------------------------------------------------------------------------------------//
-const med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type) 
+int CONNECTIVITY::getConnectivityLength(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
 //------------------------------------------------------------------------------------------------------------------//
 {
   const char * LOC = "CONNECTIVITY::getConnectivity";
@@ -814,7 +569,7 @@ const med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType,
 
   MEDSKYLINEARRAY * Connectivity;
   if (Entity==_entity) {
-    
+
     if (ConnectivityType==MED_NODAL)
       {
        calculateNodalConnectivity();
@@ -825,25 +580,25 @@ const med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType,
        calculateDescendingConnectivity();
        Connectivity=_descending;
       }
-    
+
     if (Connectivity!=NULL)
       if (Type==MED_ALL_ELEMENTS)
-       return Connectivity->getValue();
+       return Connectivity->getLength();
       else {
-       for (med_int i=0; i<_numberOfTypes; i++)
+       for (int i=0; i<_numberOfTypes; i++)
          if (_geometricTypes[i]==Type)
-           //return Connectivity->getI(i+1);
-           return Connectivity->getI(_count[i]);
+           return Connectivity->getNumberOfI(_count[i]);
        throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
       }
     else
       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
-  } else 
+  }
+  else
     if (_constituent != NULL)
-      return _constituent->getConnectivity(ConnectivityType,Entity,Type);
-  
+      return _constituent->getConnectivityLength(ConnectivityType,Entity,Type);
+
   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
-}  
+}
 
 /*!  Give morse index array to use with
      getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
@@ -854,105 +609,291 @@ const med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType,
       - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
       - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
 //-----------------------------------------------------------------------------------------------//
-const med_int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity) 
-//----0000000--------------------------------------------------------------------------------------------//
+const int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
+//-----------------------------------------------------------------------------------------------//
 {
   const char * LOC = "CONNECTIVITY::getConnectivityIndex";
   BEGIN_OF(LOC);
 
   MEDSKYLINEARRAY * Connectivity;
   if (Entity==_entity) {
-    
+
     if (ConnectivityType==MED_NODAL)
       Connectivity=_nodal;
     else
       Connectivity=_descending;
-    
+
     if (Connectivity!=NULL)
       return Connectivity->getIndex();
-    else 
+    else
       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 
-  } else 
+  } else
     if (_constituent != NULL)
       return _constituent->getConnectivityIndex(ConnectivityType,Entity);
 
   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 }
 
-/*! A DOCUMENTER */
-//--------------------------------------------------------------//
-const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
-//--------------------------------------------------------------//
+
+//-------------------------------------------------------------//
+const int* CONNECTIVITY::getPolygonsConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
+//-------------------------------------------------------------//
 {
-  const char * LOC = "CONNECTIVITY::getType";
+  const char* LOC = "CONNECTIVITY::getPolygonsConnectivity";
   BEGIN_OF(LOC);
 
-  if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
-    throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !"));
-  for (med_int i=0; i<_numberOfTypes; i++)
-    if (_geometricTypes[i]==Type)
-      return _type[i];
-  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
+  MEDSKYLINEARRAY* Connectivity;
+  if (Entity == _entity)
+    {
+      if (ConnectivityType == MED_NODAL)
+       {
+         calculateNodalConnectivity();
+         Connectivity = _polygonsNodal;
+       }
+      else
+       {
+         calculateDescendingConnectivity();
+         Connectivity = _polygonsDescending;
+       }
+      if (Connectivity != NULL)
+       return Connectivity->getValue();
+      else
+       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polygons Connectivity not defined !"));
+    }
+  else
+    if (_constituent != NULL)
+      return _constituent->getPolygonsConnectivity(ConnectivityType, Entity);
+  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 }
 
-/*!  Returns the number of elements of type <med geometrie element>.
-     Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
-//------------------------------------------------------------------------//
-med_int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
-//------------------------------------------------------------------------//
+
+//-------------------------------------------------------------//
+const int* CONNECTIVITY::getPolygonsConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
+//-------------------------------------------------------------//
 {
-  const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
+  const char* LOC = "CONNECTIVITY::getPolygonsConnectivityIndex";
   BEGIN_OF(LOC);
 
-  if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
-    throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
-  for (med_int i=0; i<_numberOfTypes; i++)
-    if (_geometricTypes[i]==Type)
-      return _type[i].getNumberOfNodes();
-  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
+  MEDSKYLINEARRAY* Connectivity;
+  if (Entity == _entity)
+    {
+      if (ConnectivityType == MED_NODAL)
+       {
+         //      calculateNodalConnectivity();
+         Connectivity = _polygonsNodal;
+       }
+      else
+       {
+         //      calculateDescendingConnectivity();
+         Connectivity = _polygonsDescending;
+       }
+      if (Connectivity != NULL)
+       return Connectivity->getIndex();
+      else
+       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polygons Connectivity not defined !"));
+    }
+  else
+    if (_constituent != NULL)
+      return _constituent->getPolygonsConnectivityIndex(ConnectivityType, Entity);
+  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 }
 
-/*!  Returns the number of geometric sub cells of <med geometrie element> type.
-not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
-//------------------------------------------------------------------------//
-med_int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
-//------------------------------------------------------------------------//
+
+/*! We suppose in this method that nodal and descending connectivities
+  are coherent.*/
+//-------------------------------------------------------------//
+int CONNECTIVITY::getNumberOfPolygons() const
+//-------------------------------------------------------------//
 {
-  if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
-    throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
-  for (med_int i=0; i<_numberOfTypes; i++)
-    if (_geometricTypes[i]==Type)
-      return _type[i].getNumberOfConstituents(1);
-  throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
+  if (_polygonsNodal != (MEDSKYLINEARRAY*) NULL)
+    return _polygonsNodal->getNumberOf();
+  else if (_polygonsDescending != (MEDSKYLINEARRAY*) NULL)
+    return _polygonsDescending->getNumberOf();
+  else
+    return 0;
 }
 
-/*!  Returns the number of elements of type <med geometrie element>.
 
-     Note :
-      - Implemented for MED_ALL_ELEMENTS
-      - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
-      - Not implemented for MED_NONE */
-//-----------------------------------------------------------------------------------//
-med_int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
-//-----------------------------------------------------------------------------------//
+//--------------------------------------------------------------//
+const int* CONNECTIVITY::getPolyhedronConnectivity(medConnectivity ConnectivityType) const
+//--------------------------------------------------------------//
 {
-  const char * LOC = "CONNECTIVITY::getNumberOf";
+  const char* LOC = "CONNECTIVITY::getPolyhedronConnectivity";
   BEGIN_OF(LOC);
 
-  MESSAGE(LOC<<" Entity = "<< Entity << ", _entity = "<<_entity);
+  if (_entity == MED_CELL) //polyhedron can only be MED_CELL
+    {
+      if (ConnectivityType == MED_NODAL)
+       {
+         ((CONNECTIVITY *)(this))->calculateNodalConnectivity();
+         if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
+           return _polyhedronNodal->getNodes();
+         else
+           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Nodal Connectivity not defined !"));
+       }
+      else
+       {
+         ((CONNECTIVITY *)(this))->calculateDescendingConnectivity();
+         if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
+           return _polyhedronDescending->getValue();
+         else
+           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Descending Connectivity not defined !"));
+       }
+    }
+  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
+}
 
-  if (Entity==_entity) {
-    if (Type==MED_NONE)
+
+//---------------------------------------------------------------//
+const int* CONNECTIVITY::getPolyhedronFacesIndex() const
+//---------------------------------------------------------------//
+{
+  const char* LOC = "CONNECTIVITY::getPolyhedronFacesIndex";
+  BEGIN_OF(LOC);
+
+  if (_entity == MED_CELL) //polyhedron can only be MED_CELL
+    {
+      //      calculateNodalConnectivity();
+      if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
+       return _polyhedronNodal->getFacesIndex();
+      else
+       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No Polyhedron in that Connectivity !"));
+    }
+  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
+}
+
+
+//--------------------------------------------------------------//
+const int* CONNECTIVITY::getPolyhedronIndex(medConnectivity ConnectivityType) const
+//--------------------------------------------------------------//
+{
+  const char* LOC = "CONNECTIVITY::getPolyhedronIndex";
+  BEGIN_OF(LOC);
+
+  if (_entity == MED_CELL) //polyhedron can only be MED_CELL
+    {
+      if (ConnectivityType == MED_NODAL)
+       {
+         //      calculateNodalConnectivity();
+         if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
+           return _polyhedronNodal->getPolyhedronIndex();
+         else
+           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Nodal Connectivity not defined !"));
+       }
+      else
+       {
+         //      calculateDescendingConnectivity();
+         if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
+           return _polyhedronDescending->getIndex();
+         else
+           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Descending Connectivity not defined !"));
+       }
+    }
+  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
+}
+
+
+/*! We suppose in this method that nodal and descending connectivities
+  are coherent.*/
+//-------------------------------------------------------------//
+int CONNECTIVITY::getNumberOfPolyhedronFaces() const
+//-------------------------------------------------------------//
+{
+  //  if (_polyhedronNodal == (POLYHEDRONARRAY*) NULL)
+  //    calculateNodalConnectivity();
+  if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
+    return _polyhedronNodal->getNumberOfFaces();
+  else
+    return 0;
+}
+
+
+/*! We suppose in this method that nodal and descending connectivities
+  are coherent.*/
+//--------------------------------------------------------------//
+int CONNECTIVITY::getNumberOfPolyhedron() const
+//--------------------------------------------------------------//
+{
+  if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
+    return _polyhedronNodal->getNumberOfPolyhedron();
+  else if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
+    return _polyhedronDescending->getNumberOf();
+  else
+    return 0;
+}
+
+
+/*! A DOCUMENTER */
+//--------------------------------------------------------------//
+const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
+//--------------------------------------------------------------//
+{
+  const char * LOC = "CONNECTIVITY::getType";
+  BEGIN_OF(LOC);
+
+  if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
+    throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !"));
+  for (int i=0; i<_numberOfTypes; i++)
+    if (_geometricTypes[i]==Type)
+      return _type[i];
+  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
+}
+
+/*!  Returns the number of elements of type %medGeometryElement.
+     Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
+//------------------------------------------------------------------------//
+int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
+//------------------------------------------------------------------------//
+{
+  const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
+  BEGIN_OF(LOC);
+
+  if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
+    throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
+  for (int i=0; i<_numberOfTypes; i++)
+    if (_geometricTypes[i]==Type)
+      return _type[i].getNumberOfNodes();
+  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
+}
+
+/*!  Returns the number of geometric sub cells of %medGeometryElement type.
+not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
+//------------------------------------------------------------------------//
+int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
+//------------------------------------------------------------------------//
+{
+  if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
+    throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
+  for (int i=0; i<_numberOfTypes; i++)
+    if (_geometricTypes[i]==Type)
+      return _type[i].getNumberOfConstituents(1);
+  throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
+}
+
+/*!  Returns the number of elements of type %medGeometryElement.
+
+     Note :
+      - Implemented for MED_ALL_ELEMENTS
+      - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
+      - Not implemented for MED_NONE */
+//-----------------------------------------------------------------------------------//
+int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
+//-----------------------------------------------------------------------------------//
+{
+  //const char * LOC = "CONNECTIVITY::getNumberOf";
+  if (Entity==_entity) {
+    if (Type==MED_NONE)
       return 0; // not defined !
     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
+    if (!existConnectivity(MED_NODAL,Entity) && !existConnectivity(MED_DESCENDING,Entity)) return 0; //case with only polygons for example
     if (Type==MED_ALL_ELEMENTS)
       return _count[_numberOfTypes]-1;
-    for (med_int i=0; i<_numberOfTypes; i++)
+    for (int i=0; i<_numberOfTypes; i++)
       if (_geometricTypes[i]==Type)
        return (_count[i+1] - _count[i]);
     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
-  } else 
+  } else
     if (_constituent != NULL)
       return _constituent->getNumberOf(Entity,Type);
 
@@ -962,25 +903,25 @@ med_int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type)
 
 /*! A DOCUMENTER */
 //--------------------------------------------------------------//
-const med_int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity, 
+const int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
                                medGeometryElement Type)
 //--------------------------------------------------------------//
 {
-  if (TypeConnectivity == MED_NODAL) 
+  if (TypeConnectivity == MED_NODAL)
     {
          calculateNodalConnectivity();
          if (Type==MED_ALL_ELEMENTS)
            return _nodal->getValue();
-         for (med_int i=0; i<_numberOfTypes; i++)
+         for (int i=0; i<_numberOfTypes; i++)
            if (_geometricTypes[i]==Type)
              return _nodal->getI(_count[i]);
-    } 
-  else 
+    }
+  else
     {
       calculateDescendingConnectivity();
       if (Type==MED_ALL_ELEMENTS)
        return _descending->getValue();
-      for (med_int i=0; i<_numberOfTypes; i++)
+      for (int i=0; i<_numberOfTypes; i++)
        if (_geometricTypes[i]==Type)
          return _descending->getI(Type);
     }
@@ -989,15 +930,15 @@ const med_int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
 
 /*! A DOCUMENTER */
 //---------------------------------------------------------------------//
-const med_int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity) 
+const int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
 //---------------------------------------------------------------------//
 {
-  if (TypeConnectivity == MED_NODAL) 
+  if (TypeConnectivity == MED_NODAL)
     {
       calculateNodalConnectivity();
       return _nodal->getIndex();
     }
-  else 
+  else
     {
       calculateDescendingConnectivity();
       return _descending->getIndex();
@@ -1006,7 +947,7 @@ const med_int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
 
 /*! Not yet implemented */
 //----------------------------------------------//
-const med_int* CONNECTIVITY:: getNeighbourhood() const
+const int* CONNECTIVITY:: getNeighbourhood() const
 //----------------------------------------------//
 {
   throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
@@ -1015,7 +956,7 @@ const med_int* CONNECTIVITY:: getNeighbourhood() const
 /*! Returns an array which contains, for each node, all cells
     arround it. */
 //-------------------------------------------------//
-const med_int* CONNECTIVITY::getReverseNodalConnectivity() 
+const int* CONNECTIVITY::getReverseNodalConnectivity()
 //-------------------------------------------------//
 {
   calculateReverseNodalConnectivity();
@@ -1025,7 +966,7 @@ const med_int* CONNECTIVITY::getReverseNodalConnectivity()
 /*!  Give index array to use with getReverseConnectivity(MED_NODAL).
      It is unusefull with MED_DESCENDING mode, because we have allways two cells.  */
 //-------------------------------------------------------//
-const med_int* CONNECTIVITY::getReverseNodalConnectivityIndex() 
+const int* CONNECTIVITY::getReverseNodalConnectivityIndex()
 //-------------------------------------------------------//
 {
   calculateReverseNodalConnectivity();
@@ -1036,30 +977,30 @@ const med_int* CONNECTIVITY::getReverseNodalConnectivityIndex()
     the 2 cells of each side. First is cell which face normal is outgoing.
     arround it. */
 //------------------------------------------------------//
-const med_int* CONNECTIVITY::getReverseDescendingConnectivity() 
+const int* CONNECTIVITY::getReverseDescendingConnectivity()
 //------------------------------------------------------//
 {
-  // it is in _constituent connectivity only if we are in MED_CELL 
+  // it is in _constituent connectivity only if we are in MED_CELL
   // (we could not for instance calculate face-edge connectivity !)
   if (_entity!=MED_CELL)
     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
 
-  // we need descending connectivity 
+  // we need descending connectivity
   calculateDescendingConnectivity();
   return _reverseDescendingConnectivity->getValue();
 }
 
-/*! calculate the reverse descending Connectivity 
+/*! calculate the reverse descending Connectivity
     and returns the index  ( A DOCUMENTER MIEUX)*/
 //-----------------------------------------------------------//
-const med_int* CONNECTIVITY::getReverseDescendingConnectivityIndex() 
+const int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
 //-----------------------------------------------------------//
 {
   // it is in _constituent connectivity only if we are in MED_CELL
-  if (_entity!=MED_CELL) 
+  if (_entity!=MED_CELL)
     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
 
-  // we need descending connectivity 
+  // we need descending connectivity
   calculateDescendingConnectivity();
   return _reverseDescendingConnectivity->getIndex();
 }
@@ -1069,15 +1010,15 @@ const med_int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
 void CONNECTIVITY::calculateNodalConnectivity()
 //--------------------------------------------//
 {
-  if (_nodal==NULL
+  if (_nodal==NULL && _polygonsNodal==NULL && _polyhedronNodal==NULL)
     {
-      if (_descending==NULL)
+      if (_descending==NULL && _polygonsDescending==NULL && _polyhedronDescending==NULL)
        throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
-      // calculate _nodal from _descending
+      // calculate _nodal, _polygonsNodal and _polyhedronNodal from _descending, _polygonsDescending and _polyhedronDescending
     }
 }
 
-/*! If not yet done, calculate the nodal Connectivity 
+/*! If not yet done, calculate the nodal Connectivity
     and the reverse nodal Connectivity*/
 //---------------------------------------------------//
 void CONNECTIVITY::calculateReverseNodalConnectivity()
@@ -1086,80 +1027,107 @@ void CONNECTIVITY::calculateReverseNodalConnectivity()
   const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
   BEGIN_OF(LOC);
 
-  if (_nodal==NULL) 
+  SCRUTE(_nodal);
+  SCRUTE(_reverseNodalConnectivity);
+
+  if (_nodal==NULL)
     calculateNodalConnectivity();
-  
+
   if(_reverseNodalConnectivity==NULL) {
 
-    med_int node_number = 0;
-    vector <vector <med_int> > reverse_connectivity; 
+    int node_number = 0;
+    vector <vector <int> > reverse_connectivity;
     reverse_connectivity.resize(_numberOfNodes+1);
-  
+
   // Treat all cells types
-  
-    for (med_int j = 0; j < _numberOfTypes; j++)
+
+    for (int j = 0; j < _numberOfTypes; j++)
       {
        // node number of the cell type
        node_number = _type[j].getNumberOfNodes();
        // treat all cells of a particular type
-       for (med_int k = _count[j]; k < _count[j+1]; k++)
+       for (int k = _count[j]; k < _count[j+1]; k++)
          // treat all nodes of the cell type
-         for (med_int local_node_number = 1; local_node_number < node_number+1; local_node_number++)
+         for (int local_node_number = 1; local_node_number < node_number+1; local_node_number++)
            {
-             med_int global_node = _nodal->getIJ(k,local_node_number);
+             int global_node = _nodal->getIJ(k,local_node_number);
              reverse_connectivity[global_node].push_back(k);
            }
       }
-  
+
+    if(_polygonsNodal != NULL && _entityDimension==2)
+      {
+       int nbOfPolygons=_polygonsNodal->getNumberOf();
+       int offset=getNumberOf(_entity,MED_ALL_ELEMENTS);
+       const int *index=_polygonsNodal->getIndex();
+       const int *value=_polygonsNodal->getValue();
+       for (int local_polyg_number = 0; local_polyg_number < nbOfPolygons; local_polyg_number++)
+         {
+           for(int i=index[local_polyg_number];i<index[local_polyg_number+1];i++)
+             reverse_connectivity[value[i-1]].push_back(offset+local_polyg_number+1);
+         }
+      }
+
+    if(_polyhedronNodal != NULL && _entityDimension==3)
+      {
+       int nbOfPolyhedra=_polyhedronNodal->getNumberOfPolyhedron();
+       int offset=getNumberOf(_entity,MED_ALL_ELEMENTS);
+       for (int local_polyh_number = 0; local_polyh_number < nbOfPolyhedra; local_polyh_number++)
+           {
+             int nbOfNodes;
+             int global_polyh_number=offset+local_polyh_number+1;
+             int *nodes=getNodesOfPolyhedron(global_polyh_number,nbOfNodes);
+             for(int i=0;i<nbOfNodes;i++)
+               reverse_connectivity[nodes[i]].push_back(global_polyh_number);
+             delete [] nodes;
+           }
+      }
+
     // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
 
   //calculate size of reverse_nodal_connectivity array
-    med_int size_reverse_nodal_connectivity  = 0;
-    for (med_int i = 1; i < _numberOfNodes+1; i++)
+    int size_reverse_nodal_connectivity  = 0;
+    for (int i = 1; i < _numberOfNodes+1; i++)
       size_reverse_nodal_connectivity += reverse_connectivity[i].size();
-  
-    //MEDSKYLINEARRAY * ReverseConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity);
-    med_int * reverse_nodal_connectivity_index = new (med_int)[_numberOfNodes+1];
-    med_int * reverse_nodal_connectivity = new (med_int)[size_reverse_nodal_connectivity];
-    //const med_int * reverse_nodal_connectivity = ReverseConnectivity->getValue();
-    //const med_int * reverse_nodal_connectivity_index = ReverseConnectivity->getIndex();
+
+    int * reverse_nodal_connectivity_index = new int[_numberOfNodes+1];
+    int * reverse_nodal_connectivity = new int[size_reverse_nodal_connectivity];
 
     reverse_nodal_connectivity_index[0] = 1;
-    for (med_int i = 1; i < _numberOfNodes+1; i++)
+    for (int i = 1; i < _numberOfNodes+1; i++)
       {
-       med_int size = reverse_connectivity[i].size(); 
+       int size = reverse_connectivity[i].size();
        reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
-       for (med_int j = 0; j < size; j++)
+       for (int j = 0; j < size; j++)
          reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
       }
-  
-    //_reverseNodalConnectivity = ReverseConnectivity;
+
     _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
                                                    reverse_nodal_connectivity_index,
-                                                   reverse_nodal_connectivity);
-    delete [] reverse_nodal_connectivity_index;
-    delete [] reverse_nodal_connectivity;
- }
+                                                   reverse_nodal_connectivity,true);
+  }
+  END_OF(LOC);
 }
 
 /*! If not yet done, calculate the Descending Connectivity */
 //-------------------------------------------------//
 void CONNECTIVITY::calculateDescendingConnectivity()
-//-------------------------------------------------//
-{
+  //-------------------------------------------------//
+  {
   const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
   BEGIN_OF(LOC);
-  
-  if (_descending==NULL)
+
+  if (_descending==NULL && _polygonsDescending==NULL && _polyhedronDescending==NULL)
     {
-      if (_nodal==NULL) {
-       MESSAGE(LOC<<"No connectivity found !");
-       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
-      }
+      if (_nodal==NULL)
+       {
+         MESSAGE(LOC<<"No connectivity found !");
+         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
+       }
       // calcul _descending from _nodal
       // we need CONNECTIVITY for constituent
-      
-      if (_constituent != NULL) 
+
+      if (_constituent != NULL)
        //      throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
        MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
 
@@ -1167,16 +1135,17 @@ void CONNECTIVITY::calculateDescendingConnectivity()
        _constituent = new CONNECTIVITY(MED_FACE);
       else if (_entityDimension == 2)
        _constituent = new CONNECTIVITY(MED_EDGE);
-      else {
-       MESSAGE(LOC<<"We are in 1D");
-       return;
-      }
+      else
+       {
+         MESSAGE(LOC<<"We are in 1D");
+         return;
+       }
 
       _constituent->_typeConnectivity = MED_NODAL;
       _constituent->_numberOfNodes = _numberOfNodes;
       // foreach cells, we built array of constituent
       int DescendingSize = 0;
-      for(int i=0; i<_numberOfTypes; i++) 
+      for(int i=0; i<_numberOfTypes; i++)
        DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
       //_descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize);
       //const int * descend_connectivity = _descending->getValue();
@@ -1192,44 +1161,49 @@ void CONNECTIVITY::calculateDescendingConnectivity()
       int * NumberOfConstituentsForeachType = new int[2];
       NumberOfConstituentsForeachType[0]=0;
       NumberOfConstituentsForeachType[1]=0;
-      for(int i=0; i<_numberOfTypes; i++) {
-       // initialize descend_connectivity_index array :
-       int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
-       for (int j=_count[i];j<_count[i+1];j++){
-         descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
-         // compute number of constituent of all cell for each type
-         for(int k=1;k<NumberOfConstituents+1;k++) {
-           medGeometryElement MEDType = _type[i].getConstituentType(1,k);
-           if (ConstituentsTypes[0]==MED_NONE) {
-             ConstituentsTypes[0]=MEDType;
-             NumberOfConstituentsForeachType[0]++;
-           } else if (ConstituentsTypes[0]==MEDType)
-             NumberOfConstituentsForeachType[0]++;
-           else if (ConstituentsTypes[1]==MED_NONE) {
-             ConstituentsTypes[1]=MEDType;
-             NumberOfConstituentsForeachType[1]++;
-           } else if (ConstituentsTypes[1]==MEDType)
-             NumberOfConstituentsForeachType[1]++;
-           else
-             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
-         }
+      for(int i=0; i<_numberOfTypes; i++)
+       {
+         // initialize descend_connectivity_index array :
+         int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
+         for (int j=_count[i];j<_count[i+1];j++)
+           {
+             descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
+             // compute number of constituent of all cell for each type
+             for(int k=1;k<NumberOfConstituents+1;k++)
+               {
+                 medGeometryElement MEDType = _type[i].getConstituentType(1,k);
+                 if (ConstituentsTypes[0]==MED_NONE)
+                   {
+                     ConstituentsTypes[0]=MEDType;
+                     NumberOfConstituentsForeachType[0]++;
+                   } else if (ConstituentsTypes[0]==MEDType)
+                     NumberOfConstituentsForeachType[0]++;
+                 else if (ConstituentsTypes[1]==MED_NONE)
+                   {
+                     ConstituentsTypes[1]=MEDType;
+                     NumberOfConstituentsForeachType[1]++;
+                   } else if (ConstituentsTypes[1]==MEDType)
+                     NumberOfConstituentsForeachType[1]++;
+                 else
+                   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
+               }
+           }
        }
-      }
 
       // we could built _constituent !
       int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
       int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
-      
+
       //_constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes);
 
-      // we use _constituent->_nodal 
+      // we use _constituent->_nodal
       //const int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
       //const int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
       int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
       int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
       ConstituentNodalConnectivityIndex[0]=1;
 
-      _constituent->_entityDimension=ConstituentsTypes[0]/100;
+      _constituent->_entityDimension = _entityDimension-1;
       if (ConstituentsTypes[1]==MED_NONE)
        _constituent->_numberOfTypes = 1;
       else
@@ -1251,7 +1225,7 @@ void CONNECTIVITY::calculateDescendingConnectivity()
       }
       delete [] ConstituentsTypes;
       delete [] NumberOfConstituentsForeachType;
-    
+
       // we need reverse nodal connectivity
       if (! _reverseNodalConnectivity)
        calculateReverseNodalConnectivity();
@@ -1260,120 +1234,128 @@ void CONNECTIVITY::calculateDescendingConnectivity()
 
       // array to keep reverse descending connectivity
       int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
-      
+
       int TotalNumberOfSubCell = 0;
-      for (int i=0; i<_numberOfTypes; i++) { // loop on all cell type
-       CELLMODEL Type = _type[i];
-       //      int NumberOfNodesPerCell = Type.getNumberOfNodes();
-       int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
-       for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
-         for (int k=1; k<=NumberOfConstituentPerCell; k++)
-           { // we loop on all sub cell of it
-             if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
-               // it is a new sub cell !
-               //          TotalNumberOfSubCell++;
-               // Which type ?
-               if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
-                 tmp_NumberOfConstituentsForeachType[0]++;
-                 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
-               } else {
-                 tmp_NumberOfConstituentsForeachType[1]++;
-                 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
-               }
-             //we have maximum two types
-               
-               descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
-               ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
-               int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
-           
-               int * NodesLists = new int[NumberOfNodesPerConstituent];
-               for (int l=0; l<NumberOfNodesPerConstituent; l++) {
-                 NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
-                 ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
-               }
-               // we use reverse_nodal_connectivity to find the other element which contain this sub cell
-
-               // all elements which contains first node of sub cell :
-               int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
-               int ReverseNodalConnectivityIndex_1 = ReverseNodalConnectivityIndex[NodesLists[0]];
-               int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
-
-               if (NumberOfCellsInList > 0)  { // we could have no element !
-                 int * CellsList = new int[NumberOfCellsInList];
-                 for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
-                   CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
-                 
-                 // foreach node in sub cell, we search elements which are in common
-                 // at the end, we must have only one !
-
-                 for (int l=1; l<NumberOfNodesPerConstituent; l++) {
-                   int NewNumberOfCellsInList = 0;
-                   int * NewCellsList = new int[NumberOfCellsInList];
-                   for (int l1=0; l1<NumberOfCellsInList; l1++)
-                     for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++) {
-                       if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
-                         // increasing order : CellsList[l1] are not in elements list of node l
-                         break;
-                       if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
-                         // we have found one
-                         NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
-                         NewNumberOfCellsInList++;
-                         break;
-                       }
-                     }
-                   NumberOfCellsInList = NewNumberOfCellsInList;
+      for (int i=0; i<_numberOfTypes; i++)
+       { // loop on all cell type
+         CELLMODEL Type = _type[i];
+         //      int NumberOfNodesPerCell = Type.getNumberOfNodes();
+         int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
+         for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
+           for (int k=1; k<=NumberOfConstituentPerCell; k++)
+             { // we loop on all sub cell of it
+               if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
+                 // it is a new sub cell !
+                 //        TotalNumberOfSubCell++;
+                 // Which type ?
+                 if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
+                   tmp_NumberOfConstituentsForeachType[0]++;
+                   TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
+                 } else {
+                   tmp_NumberOfConstituentsForeachType[1]++;
+                   TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
+                 }
+                 //we have maximum two types
+
+                 descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
+                 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
+                 int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
 
-                   delete [] CellsList;
-                   CellsList = NewCellsList;
+                 int * NodesLists = new int[NumberOfNodesPerConstituent];
+                 for (int l=0; l<NumberOfNodesPerConstituent; l++) {
+                   NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
+                   ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
                  }
-           
-                 if (NumberOfCellsInList > 0) { // We have found some elements !
-                   int CellNumber = CellsList[0];
-                   //delete [] CellsList;
-                   if (NumberOfCellsInList>1)
-                     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
-
-                   if (NumberOfCellsInList==1) {
-                     ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
-                     
-                     // we search sub cell number in this cell to not calculate it another time
-                     // which type ?
-                     CELLMODEL Type2;
-                     for (int l=0; l<_numberOfTypes; l++)
-                       if (CellNumber < _count[l+1]) {
-                         Type2=_type[l];
-                         break;
-                       }
-                     //int sub_cell_count2 = Type2.get_entities_count(1);
-                     //int nodes_cell_count2 = Type2.get_nodes_count();
-                     bool find2 = false;
-                     for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) { // on all sub cell
-                       int counter = 0;
-                       for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
-                         for (int n=1; n<=Type.getConstituentType(1,k)%100; n++) { 
-                           if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
-                             counter++;
+                 // we use reverse_nodal_connectivity to find the other element which contain this sub cell
+
+                 // all elements which contains first node of sub cell :
+                 int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
+                 int ReverseNodalConnectivityIndex_1 = getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[0]-1);
+                 //ReverseNodalConnectivityIndex[NodesLists[0]];
+                 int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
+
+                 if (NumberOfCellsInList > 0)
+                   { // we could have no element !
+                     int * CellsList = new int[NumberOfCellsInList];
+                     for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
+                       CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
+
+                     // foreach node in sub cell, we search elements which are in common
+                     // at the end, we must have only one !
+
+                     for (int l=1; l<NumberOfNodesPerConstituent; l++) {
+                       int NewNumberOfCellsInList = 0;
+                       int * NewCellsList = new int[NumberOfCellsInList];
+                       for (int l1=0; l1<NumberOfCellsInList; l1++)
+                         for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[l]-1); l2++)
+                         //for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++)
+                           {
+                             if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
+                               // increasing order : CellsList[l1] are not in elements list of node l
+                               break;
+                             if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
+                               // we have found one
+                               NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
+                               NewNumberOfCellsInList++;
+                               break;
+                             }
+                           }
+                       NumberOfCellsInList = NewNumberOfCellsInList;
+
+                       delete [] CellsList;
+                       CellsList = NewCellsList;
+                     }
+
+                     if (NumberOfCellsInList > 0) { // We have found some elements !
+                       int CellNumber = CellsList[0];
+                       //delete [] CellsList;
+                       if (NumberOfCellsInList>1)
+                         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
+
+                       if (NumberOfCellsInList==1)
+                         {
+                           ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
+
+                           // we search sub cell number in this cell to not calculate it another time
+                           // which type ?
+                           CELLMODEL Type2;
+                           for (int l=0; l<_numberOfTypes; l++)
+                             if (CellNumber < _count[l+1]) {
+                               Type2=_type[l];
+                               break;
+                             }
+                           //int sub_cell_count2 = Type2.get_entities_count(1);
+                           //int nodes_cell_count2 = Type2.get_nodes_count();
+                           bool find2 = false;
+                           for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) { // on all sub cell
+                             int counter = 0;
+                             for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
+                               for (int n=1; n<=Type.getConstituentType(1,k)%100; n++)
+                                 {
+                                   if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
+                                     counter++;
+                                 }
+                             if (counter==Type.getConstituentType(1,k)%100)
+                               {
+                                 descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
+                                 find2 = true;
+                               }
+                             if (find2)
+                               break;
+                           }
+                           if (!find2)
+                             MESSAGE(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
                          }
-                       if (counter==Type.getConstituentType(1,k)%100) {
-                         descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
-                         find2 = true;
-                       }
-                       if (find2)
-                         break;
+                     } else {
+                       ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
                      }
-                     if (!find2)
-                       MESSAGE(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
-                   } 
-                 } else {
-                   ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
-                 }
-                 delete [] CellsList;
-               }
+                     delete [] CellsList;
+                   }
 
-               delete [] NodesLists;
+                 delete [] NodesLists;
+               }
              }
-           }
-      }
+       }
       // we adjust _constituent
       int NumberOfConstituent=0;
       int SizeOfConstituentNodal=0;
@@ -1388,7 +1370,7 @@ void CONNECTIVITY::calculateDescendingConnectivity()
       int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
       int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
       ConstituentNodalIndex[0]=1;
-      // we build _reverseDescendingConnectivity 
+      // we build _reverseDescendingConnectivity
       //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
       //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
       //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
@@ -1432,6 +1414,9 @@ void CONNECTIVITY::calculateDescendingConnectivity()
 
       delete [] ConstituentNodalConnectivityIndex;
       delete [] ConstituentNodalConnectivity;
+      delete [] ReverseDescendingConnectivityValue;
+      _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
+      delete [] tmp_NumberOfConstituentsForeachType;
 
       _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
                                        DescendingSize,
@@ -1439,42 +1424,185 @@ void CONNECTIVITY::calculateDescendingConnectivity()
                                        descend_connectivity);
       delete [] descend_connectivity_index;
       delete [] descend_connectivity;
-      _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,
-                                                          2*NumberOfConstituent,
-                                                          reverseDescendingConnectivityIndex,
-                                                          reverseDescendingConnectivityValue);
-      delete [] reverseDescendingConnectivityIndex;
+
+      ////
+      vector<int> PolyDescending;
+      vector<int> Reversedescendingconnectivityvalue(reverseDescendingConnectivityValue,reverseDescendingConnectivityValue + 2*NumberOfConstituent);
+      vector<int> Reversedescendingconnectivityindex(reverseDescendingConnectivityIndex,reverseDescendingConnectivityIndex + NumberOfConstituent);
       delete [] reverseDescendingConnectivityValue;
+      delete [] reverseDescendingConnectivityIndex;
 
-      _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
-      delete [] tmp_NumberOfConstituentsForeachType;
-      
-      //delete _constituent->_nodal;
-      _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
-                                                SizeOfConstituentNodal,
-                                                ConstituentNodalIndex,
-                                                ConstituentNodalValue);
-    
-      delete [] ConstituentNodalIndex;
+
+      // polygons (2D mesh)
+
+      vector<int> Constituentnodalvalue(ConstituentNodalValue,ConstituentNodalValue + SizeOfConstituentNodal);
+      vector<int> Constituentnodalindex(ConstituentNodalIndex,ConstituentNodalIndex + NumberOfConstituent+1);
       delete [] ConstituentNodalValue;
+      delete [] ConstituentNodalIndex;
+      int NumberOfNewSeg = 0;
 
-      delete [] ReverseDescendingConnectivityValue;
+      for (int i=0; i <getNumberOfPolygons(); i++) // for each polygon
+       {
+         const int * vector_begin = &_polygonsNodal->getValue()[_polygonsNodal->getIndex()[i]-1];
+         int vector_size = _polygonsNodal->getIndex()[i+1]-_polygonsNodal->getIndex()[i]+1;
+         vector<int> myPolygon(vector_begin, vector_begin+vector_size);
+         myPolygon[myPolygon.size()-1] = myPolygon[0]; // because first and last point make a segment
+
+         for (int j=0; j<myPolygon.size()-1; j++) // for each segment of polygon
+           {
+             MEDMODULUSARRAY segment_poly(2,&myPolygon[j]);
+             int ret_compare = 0;
+
+             // we search it in existing segments
+
+             for (int k=0; k<Constituentnodalindex.size()-1; k++)
+               {
+                 MEDMODULUSARRAY segment(2,&Constituentnodalvalue[0] + Constituentnodalindex[k]-1);
+                 ret_compare = segment_poly.compare(segment);
+                 if (ret_compare) // segment_poly already exists
+                   {
+                     PolyDescending.push_back(ret_compare*(k+1)); // we had it to the connectivity
+                     insert_vector(Reversedescendingconnectivityvalue, 2*k+1, i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polygon i to reverse descending connectivity for segment_poly (in 2sd place)
+                     break;
+                   }
+               }
+
+             // segment_poly must be created
+
+             if (!ret_compare)
+               {
+                 NumberOfNewSeg++;
+                 PolyDescending.push_back(NumberOfConstituent+NumberOfNewSeg); // we had it to the connectivity
+                 Constituentnodalvalue.push_back(segment_poly[0]);
+                 Constituentnodalvalue.push_back(segment_poly[1]);
+                 Constituentnodalindex.push_back(Constituentnodalindex[Constituentnodalindex.size()-1] + 2); // we have only segments
+                 insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+NumberOfNewSeg-1), i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polygon i to reverse descending connectivity for segment_poly (in 1st place)
+                 insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewSeg-1, 2*(NumberOfConstituent+NumberOfNewSeg-1)+1); // idem with index
+               }
+           }
+       }
+
+      if (getNumberOfPolygons() > 0)
+       {
+         _polygonsDescending = new MEDSKYLINEARRAY(getNumberOfPolygons(),_polygonsNodal->getLength(),_polygonsNodal->getIndex(),&PolyDescending[0]); // index are the same for polygons nodal and descending connectivities
+
+         NumberOfConstituent += NumberOfNewSeg;
+         SizeOfConstituentNodal += 2*NumberOfNewSeg;
+         _constituent->_count[1] = NumberOfConstituent+1;
+       }
+
+
+      // polyhedron (3D mesh)
+
+      vector<int> Constituentpolygonsnodalvalue;
+      vector<int> Constituentpolygonsnodalindex(1,1);
+      int NumberOfNewFaces = 0; // by convention new faces are polygons
+
+      for (int i=0; i<getNumberOfPolyhedron(); i++) // for each polyhedron
+       {
+         // we create a POLYHEDRONARRAY containing only this polyhedra
+         int myNumberOfFaces = _polyhedronNodal->getPolyhedronIndex()[i+1]-_polyhedronNodal->getPolyhedronIndex()[i];
+         int myNumberOfNodes = _polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i+1]-1]-_polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i]-1];
+         POLYHEDRONARRAY myPolyhedra(1,myNumberOfFaces,myNumberOfNodes);
+         vector<int> myFacesIndex(_polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1, _polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1 + myNumberOfFaces+1);
+         for (int j=myNumberOfFaces; j>=0; j--)
+           myFacesIndex[j] -= myFacesIndex[0]-1;
+         myPolyhedra.setFacesIndex(&myFacesIndex[0]);
+         myPolyhedra.setNodes(_polyhedronNodal->getNodes() + _polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i]-1]-1);
+
+         for (int j=0; j<myPolyhedra.getNumberOfFaces(); j++) // for each face of polyhedra
+           {
+             int myFaceNumberOfNodes = myPolyhedra.getFacesIndex()[j+1]-myPolyhedra.getFacesIndex()[j];
+             MEDMODULUSARRAY face_poly(myFaceNumberOfNodes,myPolyhedra.getNodes() + myPolyhedra.getFacesIndex()[j]-1);
+             int ret_compare = 0;
+
+             // we search it in existing faces
+
+             // we search first in TRIA3 and QUAD4
+             if (myFaceNumberOfNodes == 3 || myFaceNumberOfNodes == 4)
+               {
+                 int Begin = -1; // first TRIA3 or QUAD4
+                 int Number = 0; // number of TRIA3 or QUAD4
+                 for (int k=0; k<Constituentnodalindex.size()-1; k++)
+                   if (Constituentnodalindex[k+1]-Constituentnodalindex[k] == myFaceNumberOfNodes)
+                     {
+                       if (Begin == -1)
+                         Begin = k;
+                       Number++;
+                     }
+
+                 for (int k=0; k<Number; k++)
+                   {
+                     MEDMODULUSARRAY face(myFaceNumberOfNodes,&Constituentnodalvalue[0] + Constituentnodalindex[Begin+k]-1);
+                     ret_compare = face_poly.compare(face);
+                     if (ret_compare)
+                       {
+                         PolyDescending.push_back(ret_compare*(Begin+k+1)); // we had it to the connectivity
+                         insert_vector(Reversedescendingconnectivityvalue, 2*(Begin+k)+1, i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polyhedra i to reverse descending connectivity for face_poly (in 2sd place)
+                         break;
+                       }
+                   }
+               }
+
+             // we search last in POLYGONS
+             if (!ret_compare)
+               {
+                 for (int k=0; k<static_cast<int>(Constituentpolygonsnodalindex.size())-1; k++) // we must cast the unsigned int into int before doing -1
+                   {
+                     if (Constituentpolygonsnodalindex[k+1]-Constituentpolygonsnodalindex[k] == myFaceNumberOfNodes)
+                       {
+                         MEDMODULUSARRAY face(myFaceNumberOfNodes,&Constituentpolygonsnodalvalue[0] + Constituentpolygonsnodalindex[k]-1);
+                         ret_compare = face_poly.compare(face);
+                         if (ret_compare)
+                           {
+                             PolyDescending.push_back(ret_compare*(NumberOfConstituent+k+1)); // we had it to the connectivity
+                             insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+k)+1, i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polyhedra i to reverse descending connectivity for face_poly (in 2sd place)
+                             break;
+                           }
+                       }
+                   }
+               }
+
+             // if not found, face_poly must be created
+
+             if (!ret_compare)
+               {
+                 NumberOfNewFaces++;
+                 PolyDescending.push_back(NumberOfConstituent+NumberOfNewFaces); // we had it to the connectivity
+                 for (int k=0; k<myFaceNumberOfNodes; k++)
+                   Constituentpolygonsnodalvalue.push_back(face_poly[k]);
+                 Constituentpolygonsnodalindex.push_back(Constituentpolygonsnodalindex[Constituentpolygonsnodalindex.size()-1]+myFaceNumberOfNodes);
+                 insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+NumberOfNewFaces-1), i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polyhedra i to reverse descending connectivity for face_poly (in 1st place)
+                 insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewFaces-1, 2*(NumberOfConstituent+NumberOfNewFaces-1)+1); // idem with index
+               }
+           }
+       }
+
+      if (getNumberOfPolyhedron() > 0)
+       {
+         _polyhedronDescending = new MEDSKYLINEARRAY(getNumberOfPolyhedron(),_polyhedronNodal->getNumberOfFaces(),_polyhedronNodal->getPolyhedronIndex(),&PolyDescending[0]); // polyhedron index are the same for nodal and descending connectivities
+
+         if (_constituent->_polygonsNodal != NULL)
+           delete [] _constituent->_polygonsNodal;
+         _constituent->_polygonsNodal = new MEDSKYLINEARRAY(Constituentpolygonsnodalindex.size()-1,Constituentpolygonsnodalvalue.size(),&Constituentpolygonsnodalindex[0],&Constituentpolygonsnodalvalue[0]);
+       }
+
+      // delete _constituent->_nodal
+      _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
+                                                SizeOfConstituentNodal,
+                                                &Constituentnodalindex[0],
+                                                &Constituentnodalvalue[0]);
+
+      int NumberOfConstituentWithPolygons = NumberOfConstituent + NumberOfNewFaces;
+      Reversedescendingconnectivityindex.push_back(Reversedescendingconnectivityindex[Reversedescendingconnectivityindex.size()-1]+2); // we complete the index
+      _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituentWithPolygons+1,
+                                                          2*NumberOfConstituentWithPolygons,
+                                                          &Reversedescendingconnectivityindex[0],
+                                                          &Reversedescendingconnectivityvalue[0]);
+      ////
     }
   END_OF(LOC);
-}
-
-//-----------------------------------------------------------------------------------------//
-//  void CONNECTIVITY::calculateReverseDescendingConnectivity(CONNECTIVITY *myConnectivity)
-//-----------------------------------------------------------------------------------------//
-//  {
-//    if (_entity==MED_CELL)
-//      throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
-//        Connectivity : we could not do that with MED_CELL entity !");
-//    if (myConnectivity->getEntity()!=MED_CELL)
-//      throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
-//             Connectivity : we need MED_CELL connectivity !");
-//    return;
-//  }
+  }
 
 /*! Not implemented yet */
 //--------------------------------------------------------------------//
@@ -1492,7 +1620,7 @@ void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
 }
 
 
-/*! 
+/*!
   Returns the geometry of an element given by its entity type & its global number.
 
   Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
@@ -1503,11 +1631,19 @@ medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalN
 {
   const char * LOC = "medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity, int globalNumber) const : ";
   BEGIN_OF(LOC);
+  int globalNumberMin = 1;
+  int globalNumberMax ;
+
+  if (_entity==Entity) globalNumberMax = _count[_numberOfTypes];
+  else if (_constituent!=NULL) globalNumberMax = _constituent->_count[_constituent->_numberOfTypes];
+  else
+    throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
 
   // The globalNumber must verify  : 1 <=  globalNumber < _count[_numberOfTypes] (== totalNumberOfElement+1)
-  if ( (globalNumber < 1) || (globalNumber >  _count[_numberOfTypes]-1 )  )
-    throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |" 
-                                 << 1 <<"| and <= |" << _count[_numberOfTypes]-1 ));
+
+  if ( (globalNumber < globalNumberMin) || (globalNumber >  globalNumberMax-1 )  )
+    throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |"
+                                 << globalNumberMin <<"| and <= |" << globalNumberMax-1 ));
 
   if (_entity==Entity) {
     for(int i=1; i<=_numberOfTypes;i++)
@@ -1523,6 +1659,47 @@ medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalN
   END_OF(LOC);
 }
 
+/*!
+  Idem as getElementType method except that it manages polygon and polyhedron.
+*/
+//--------------------------------------------------------------------//
+medGeometryElement CONNECTIVITY::getElementTypeWithPoly(medEntityMesh Entity,int globalNumber) const
+//--------------------------------------------------------------------//
+{
+  int globalNumberMaxOfClassicType;
+  if(_entity==Entity)
+    {
+      globalNumberMaxOfClassicType=_count[_numberOfTypes];
+      if(globalNumber>=1)
+       {
+         if(globalNumber<globalNumberMaxOfClassicType)
+           {
+             for(int i=1; i<=_numberOfTypes;i++)
+               if (globalNumber<_count[i])
+                 return _geometricTypes[i-1];
+             throw MEDEXCEPTION("never happens just for compilo");
+           }
+         else
+           {
+             int localNumberInPolyArray=globalNumber-globalNumberMaxOfClassicType;
+             int nbOfPol=getNumberOfElementOfPolyType(_entity);
+             if(localNumberInPolyArray<nbOfPol)
+               return getPolyTypeRelativeTo();
+             else
+               throw MEDEXCEPTION("getElementTypeWithPoly : unexisting element type");
+           }
+       }
+      else
+       throw MEDEXCEPTION("getElementTypeWithPoly : globalNumber < 1");
+    }
+  else
+    {
+      if(_constituent!=NULL)
+       return _constituent->getElementTypeWithPoly(Entity,globalNumber);
+      else
+       throw MEDEXCEPTION("getElementTypeWithPoly : Entity not defined !");
+    }
+}
 
 ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
 {
@@ -1547,7 +1724,7 @@ ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
        default:
            os << "Unknown";
            break;
-    } 
+    }
     os  << " -------------\n\nMedConnectivity : ";
     switch (co._typeConnectivity)
     {
@@ -1559,21 +1736,21 @@ ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
            break;
        default:
            break;
-    } 
+    }
     os << "Entity dimension : " << co._entityDimension << endl;
     os << "Number of nodes : " << co._numberOfNodes << endl;
     os << "Number of types : " << co._numberOfTypes << endl;
     for (int i=0; i!=co._numberOfTypes ; ++i)
-       os << "  -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : " 
+       os << "  -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : "
           << co._count[i+1]-co._count[i] << " elements" << endl;
 
     if (co._typeConnectivity == MED_NODAL )
     {
-       for (int i=0; i<co._numberOfTypes; i++) 
+       for (int i=0; i<co._numberOfTypes; i++)
        {
            os << endl << co._type[i].getName() << " : " << endl;
            int numberofelements = co._count[i+1]-co._count[i];
-           const med_int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, co._geometricTypes[i]);
+           const int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, co._geometricTypes[i]);
            int numberofnodespercell = co._geometricTypes[i]%100;
            for (int j=0;j<numberofelements;j++)
            {
@@ -1583,20 +1760,91 @@ ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
                os << endl;
            }
        }
+
+       os << endl << "MED_POLYGON : " << co.getNumberOfPolygons() << " polygons" << endl;
+       if (co.getNumberOfPolygons() > 0)
+         {
+           const int* polygonsconnectivity = co.getPolygonsConnectivity(MED_NODAL,co.getEntity());
+           const int* polygonsconnectivityindex = co.getPolygonsConnectivityIndex(MED_NODAL,co.getEntity());
+
+           for (int i=0; i<co.getNumberOfPolygons(); i++)
+             {
+               int numberofnodesperpolygon = polygonsconnectivityindex[i+1]-polygonsconnectivityindex[i];
+               for (int j=0; j<numberofnodesperpolygon; j++)
+                 os << polygonsconnectivity[polygonsconnectivityindex[i]-1+j] << " ";
+               os << endl;
+             }
+         }
+
+       os << endl << "MED_POLYHEDRA : " << co.getNumberOfPolyhedron() << " polyhedron" << endl;
+       if (co.getNumberOfPolyhedron() > 0)
+         {
+           const int* polyhedronconnectivity = co.getPolyhedronConnectivity(MED_NODAL);
+           const int* polyhedronfacesindex = co.getPolyhedronFacesIndex();
+           const int* polyhedronindex = co.getPolyhedronIndex(MED_NODAL);
+
+           for (int i=0; i<co.getNumberOfPolyhedron(); i++)
+             {
+               int numberoffacesperpolyhedra = polyhedronindex[i+1]-polyhedronindex[i];
+               for (int j=0; j<numberoffacesperpolyhedra; j++)
+                 {
+                   int numberofnodesperface = polyhedronfacesindex[polyhedronindex[i]-1+j+1]-polyhedronfacesindex[polyhedronindex[i]-1+j];
+                   for (int k=0; k<numberofnodesperface; k++)
+                     os << polyhedronconnectivity[polyhedronfacesindex[polyhedronindex[i]-1+j]-1+k] << " ";
+                   if (j != numberoffacesperpolyhedra-1) os << "| ";
+                 }
+               os << endl;
+             }
+         }
+
     }
     else if (co._typeConnectivity == MED_DESCENDING)
     {
        int numberofelements = co.getNumberOf( co._entity , MED_ALL_ELEMENTS);
-       const med_int *connectivity =  co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
-       const int *connectivity_index =  co.getConnectivityIndex( co._typeConnectivity, co._entity );
-       
-       for ( int j=0; j!=numberofelements; j++ ) 
-       {
-           os << "Element " << j+1 << " : ";
-           for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
-               os << connectivity[k-1] << " ";
-           os << endl;
-       }
+       if (numberofelements > 0)
+         {
+           const int *connectivity =  co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
+           const int *connectivity_index =  co.getConnectivityIndex( co._typeConnectivity, co._entity );
+
+           for ( int j=0; j!=numberofelements; j++ )
+             {
+               os << "Element " << j+1 << " : ";
+               for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
+                 os << connectivity[k-1] << " ";
+               os << endl;
+             }
+         }
+
+       os << endl << "MED_POLYGON : " << co.getNumberOfPolygons() << " polygons" << endl;
+       if (co.getNumberOfPolygons() > 0)
+         {
+           const int* polygonsconnectivity = co.getPolygonsConnectivity(MED_DESCENDING,co.getEntity());
+           const int* polygonsconnectivityindex = co.getPolygonsConnectivityIndex(MED_DESCENDING,co.getEntity());
+
+           for (int i=0; i<co.getNumberOfPolygons(); i++)
+             {
+               int numberofedgesperpolygon = polygonsconnectivityindex[i+1]-polygonsconnectivityindex[i];
+               for (int j=0; j<numberofedgesperpolygon; j++)
+                 os << polygonsconnectivity[polygonsconnectivityindex[i]-1+j] << " ";
+               os << endl;
+             }
+         }
+
+       os << endl << "MED_POLYHEDRA : " << co.getNumberOfPolyhedron() << " polyhedron" << endl;
+       if (co.getNumberOfPolyhedron() > 0)
+         {
+           const int* polyhedronconnectivity = co.getPolyhedronConnectivity(MED_DESCENDING);
+           const int* polyhedronindex = co.getPolyhedronIndex(MED_DESCENDING);
+
+           for (int i=0; i<co.getNumberOfPolyhedron(); i++)
+             {
+               int numberoffacesperpolyhedra = polyhedronindex[i+1]-polyhedronindex[i];
+               for (int j=0; j<numberoffacesperpolyhedra; j++)
+                 os << polyhedronconnectivity[polyhedronindex[i]-1+j] << " ";
+               os << endl;
+             }
+         }
+
     }
 
     if (co._constituent)
@@ -1604,3 +1852,307 @@ ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
 
     return os;
 }
+
+/*!
+  method that adds to vector 'nodes' all the nodes of polyhedron with id 'polyhedronId'.
+  WARNING the returned pointer should be deallocated. Returned nodes and polyhedronId are in form [1,...]
+ */
+int *CONNECTIVITY::getNodesOfPolyhedron(int polyhedronId, int& lgthOfTab) const
+{
+  const int *nodes=getPolyhedronConnectivity(MED_EN::MED_NODAL);
+  const int *faces=getPolyhedronFacesIndex();
+  const int *glob=getPolyhedronIndex(MED_EN::MED_NODAL);
+  int offsetWithClassicType=getNumberOf(_entity,MED_ALL_ELEMENTS);
+  set<int> retInSet;
+  int startFace=glob[polyhedronId-offsetWithClassicType-1]-1;
+  int endFace=glob[polyhedronId-offsetWithClassicType]-1;
+  int i;
+  for(i=startFace;i!=endFace;i++)
+    {
+      for(int j=faces[i];j<faces[i+1];j++)
+       retInSet.insert(nodes[j-1]);
+    }
+  lgthOfTab=retInSet.size();
+  int *ret=new int[lgthOfTab];
+  set<int>::iterator iter=retInSet.begin();
+  i=0;
+  for(iter=retInSet.begin();iter!=retInSet.end();iter++)
+    ret[i++]=*iter;
+  return ret;
+}
+
+/*!
+  Idem as MESH::getNodesOfPolyhedron except that returned nodes are sorted by face. 'nbOfNodesPerFaces' is an array of size 'nbOfFaces'.
+  Returned int** has a size of 'nbOfNodesPerFaces' too, and for each element j in int** the size is  nbOfNodesPerFaces[j].
+  Warning both returned 'nbOfNodesPerFaces' and returned value should be deallocated. Returned nodes and 'polyhedronId' are in form [1,...]
+ */
+int **CONNECTIVITY::getNodesPerFaceOfPolyhedron(int polyhedronId, int& nbOfFaces, int* & nbOfNodesPerFaces) const
+{
+  int i;
+  const int *nodes=getPolyhedronConnectivity(MED_EN::MED_NODAL);
+  const int *faces=getPolyhedronFacesIndex();
+  const int *glob=getPolyhedronIndex(MED_EN::MED_NODAL);
+  int offsetWithClassicType=getNumberOf(_entity,MED_ALL_ELEMENTS);
+
+  int startFace=glob[polyhedronId-offsetWithClassicType-1]-1;
+  nbOfFaces=glob[polyhedronId-offsetWithClassicType]-startFace-1;
+  nbOfNodesPerFaces=new int[nbOfFaces];
+  int **ret=new int* [nbOfFaces];
+  for(i=0;i<nbOfFaces;i++)
+    {
+      int startNode=faces[startFace+i]-1;
+      nbOfNodesPerFaces[i]=faces[startFace+i+1]-startNode-1;
+      ret[i]=(int *)(nodes)+startNode;
+    }
+  return ret;
+}
+
+int CONNECTIVITY::getNumberOfTypesWithPoly(MED_EN::medEntityMesh Entity) const
+{
+  if (_entity==Entity)
+    return _numberOfTypes+getNumberOfPolyType();
+  else if (_constituent!=NULL)
+    return _constituent->getNumberOfTypesWithPoly(Entity);
+  else
+    return 0;
+}
+
+int CONNECTIVITY::getNumberOfPolyType()  const
+{
+  if (_entity==MED_CELL && _entityDimension==3)
+    {
+      if(getNumberOfPolyhedron()>0)
+       return 1;
+    }
+  else if ((_entity==MED_CELL && _entityDimension==2) || (_entity==MED_FACE && _entityDimension==2))
+    if(getNumberOfPolygons()>0)
+      return 1;
+  return 0;
+}
+
+int CONNECTIVITY::getNumberOfElementOfPolyType(MED_EN::medEntityMesh Entity)  const
+{
+  if(_entity==Entity)
+    {
+      if (_entity==MED_CELL && _entityDimension==3)
+       return getNumberOfPolyhedron();
+      else if ((_entity==MED_CELL && _entityDimension==2) || (_entity==MED_FACE && _entityDimension==2))
+       return getNumberOfPolygons();
+      return 0;
+    }
+  else
+    {
+      if(_constituent!=NULL)
+       return _constituent->getNumberOfElementOfPolyType(Entity);
+      else
+       throw MEDEXCEPTION("_constituent required to evaluate getNumberOfElementOfPolyType");
+    }
+}
+
+/*
+  Method equivalent to getGeometricTypes except that it includes not only classical Types but polygons/polyhedra also.
+  WARNING the returned array MUST be deallocated.
+ */
+MED_EN::medGeometryElement * CONNECTIVITY::getGeometricTypesWithPoly(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
+{
+  if(_entity==Entity)
+    {
+      int nbOfTypesTotal=getNumberOfTypesWithPoly(Entity);
+      int nbOfTypesWithoutPoly=getNumberOfTypes(Entity);
+      medGeometryElement *ret=new medGeometryElement[nbOfTypesTotal];
+      memcpy(ret,getGeometricTypes(Entity),nbOfTypesWithoutPoly*sizeof(medGeometryElement));
+      if(nbOfTypesTotal>nbOfTypesWithoutPoly)
+       {
+         if (Entity==MED_CELL && _entityDimension==3)
+           ret[nbOfTypesWithoutPoly]=MED_EN::MED_POLYHEDRA;
+         else if((Entity==MED_CELL && _entityDimension==2) || (Entity==MED_FACE && _entityDimension==2))
+           ret[nbOfTypesWithoutPoly]=MED_EN::MED_POLYGON;
+       }
+      return ret;
+    }
+  else
+    {
+      if(_constituent)
+       return _constituent->getGeometricTypesWithPoly(Entity);
+      throw MEDEXCEPTION("_constituent required to evaluate getGeometricTypesWithPoly");
+    }
+}
+
+/*
+  Method used in CalculateDescendingConnectivity. So it's typically a private method.
+  The aim of this method is to hide to CalculateDescendingConnectivity algorithm the fact that in reverse connectivity polygons and polyhedrons
+  are not in separate data structure, contrary to not reverse connectivities.
+ */
+int CONNECTIVITY::getIndexOfEndClassicElementInReverseNodal(const int *reverseNodalValue, const int *reverseNodalIndex, int rk)  const
+{
+  int nbOfLastElt=getNumberOf(_entity,MED_ALL_ELEMENTS);
+  int ret=reverseNodalIndex[rk];
+  for(int i=reverseNodalIndex[rk];i<reverseNodalIndex[rk+1];i++)
+    {
+      if(reverseNodalValue[i-1]<=nbOfLastElt)
+       ret++;
+    }
+  return ret;
+}
+
+/*
+  Method that inverts the face with faceId 'faceId' in the data structure. If it's a polygon face 'faceId' is a value between 1 and nbOfPolygons.
+  This method has to be applied on a instance of MEDMEM::CONNECTIVITY with _entityDimension==3.
+  WARNING calculateDescendingConnectivity must have been called before.
+ */
+void CONNECTIVITY::invertConnectivityForAFace(int faceId, const int *nodalConnForFace, bool polygonFace)
+{
+  // we have always 2 neighbourings
+  int cell1 = _reverseDescendingConnectivity->getIJ(faceId,1);
+  int cell2 = _reverseDescendingConnectivity->getIJ(faceId,2);
+
+  if (cell2 != 0) { // we are not on border, update compulsary. If we aren't on border no update necessary so WARNING because user specified a bad oriented face
+    // Updating _reverseDescendingConnectivity
+    _reverseDescendingConnectivity->setIJ(faceId,1,cell2);
+    _reverseDescendingConnectivity->setIJ(faceId,2,cell1);
+    // Updating _constituent->_nodal because of reversity
+    MEDSKYLINEARRAY *currentNodal=(!polygonFace)?_constituent->_nodal:_constituent->_polygonsNodal;
+    MEDSKYLINEARRAY *currentDescending=(!polygonFace)?_descending:_polygonsDescending;
+    const int *descendingNodalIndex=(!polygonFace)?_constituent->_nodal->getIndex():_constituent->_polygonsNodal->getIndex();
+    const int *newDescendingIndex=(!polygonFace)?_descending->getIndex():_polygonsDescending->getIndex();
+    for(int iarray=1;iarray<=(descendingNodalIndex[faceId]-descendingNodalIndex[faceId-1]);iarray++)
+      currentNodal->setIJ(faceId,iarray,nodalConnForFace[iarray-1]);
+
+    // Updating _descending for cell1 and cell2
+    for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
+      if (currentDescending->getIndexValue(iface)==faceId)
+       currentDescending->setIndexValue(iface,-faceId);
+      else if (currentDescending->getIndexValue(iface)==-faceId)
+       currentDescending->setIndexValue(iface,faceId);
+
+    for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
+      if (currentDescending->getIndexValue(iface)==faceId)
+       currentDescending->setIndexValue(iface,-faceId);
+      else if (_descending->getIndexValue(iface)==-faceId)
+       currentDescending->setIndexValue(iface,faceId);
+  }
+}
+
+/*
+  Method with 2 output : the connectivity required and its length 'lgth'. This method gives the connectivity independently it is a polygons/polyhedrons or normal
+  element.
+ */
+const int * CONNECTIVITY::getConnectivityOfAnElementWithPoly(MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity, int Number, int &lgth)
+{
+  if(Entity==MED_EN::MED_NODE)
+    throw  MEDEXCEPTION("No connectivity attached to a node entity");
+  if(Entity==_entity)
+    {
+      if(_entity==MED_EDGE && _entityDimension==1)
+       {
+         const int * newConstituentValue = 0;
+         const int * newConstituentIndex = 0;
+         newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
+         newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
+         lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
+         return newConstituentValue+newConstituentIndex[Number-1]-1;
+       }
+      int nbOfClassicalElements=getNumberOf(_entity,MED_ALL_ELEMENTS);
+      if((_entity==MED_FACE && _entityDimension==2) || (_entity==MED_CELL && _entityDimension==2))
+       {
+         const int * newConstituentValue = 0;
+         const int * newConstituentIndex = 0;
+         if(Number<=nbOfClassicalElements)
+           {
+             newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
+             newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
+             lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
+             return newConstituentValue+newConstituentIndex[Number-1]-1;
+           }
+         else
+           {
+             int localNumber=Number-nbOfClassicalElements-1;
+             if(localNumber<getNumberOfPolygons())
+               {
+                 newConstituentValue = getPolygonsConnectivity(ConnectivityType,Entity);
+                 newConstituentIndex = getPolygonsConnectivityIndex(ConnectivityType,Entity);
+                 lgth=newConstituentIndex[localNumber+1]-newConstituentIndex[localNumber];
+                 return newConstituentValue+newConstituentIndex[localNumber]-1;
+               }
+             else
+               throw  MEDEXCEPTION("Unknown number");
+           }
+       }
+      else if(_entity==MED_CELL && _entityDimension==3)
+       {
+         const int * newConstituentValue = 0;
+         const int * newConstituentIndex = 0;
+         if(Number<=nbOfClassicalElements)
+           {
+             newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
+             newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
+             lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
+             return newConstituentValue+newConstituentIndex[Number-1]-1;
+           }
+         else
+           {
+             int localNumber=Number-nbOfClassicalElements-1;
+             if(localNumber<getNumberOfPolyhedron())
+               {
+                 if(ConnectivityType==MED_NODAL)
+                   throw  MEDEXCEPTION("NODAL Connectivity required for a polyhedron");
+                 newConstituentValue = _polyhedronDescending->getValue();
+                 newConstituentIndex = _polyhedronDescending->getIndex();
+                 lgth=newConstituentIndex[localNumber+1]-newConstituentIndex[localNumber];
+                 return newConstituentValue+newConstituentIndex[localNumber]-1;
+               }
+             else
+               throw  MEDEXCEPTION("Unknown number");
+           }
+       }
+      else
+       throw MEDEXCEPTION("No connectivity available");
+    }
+  else
+    {
+      if(_constituent==NULL)
+       calculateDescendingConnectivity();
+      return _constituent->getConnectivityOfAnElementWithPoly(ConnectivityType,Entity,Number,lgth);
+    }
+}
+
+int CONNECTIVITY::getNumberOfElementsWithPoly(MED_EN::medEntityMesh Entity, MED_EN::medGeometryElement Type) const
+{
+  if(Entity==MED_EN::MED_NODE)
+    return _numberOfNodes;
+  if(Entity==_entity)
+    {
+      if(Type==MED_POLYGON || Type==MED_POLYHEDRA)
+       return getNumberOfElementOfPolyType(_entity);
+      else
+       return getNumberOf(_entity,Type);
+    }
+  else
+    {
+      if(_constituent)
+       return _constituent->getNumberOfElementsWithPoly(Entity,Type);
+      else
+       throw MEDEXCEPTION("CONNECTIVITY::getNumberOfElementsWithPoly : _constituent needed");
+    }
+}
+
+/*!
+  Perform a deep comparison of the 2 connectivities in NODAL mode and on all elements.
+*/
+bool CONNECTIVITY::deepCompare(const CONNECTIVITY& other) const
+{
+  CONNECTIVITY* temp=(CONNECTIVITY* )this;
+  const int *conn1=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
+  int size1=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
+  temp=(CONNECTIVITY* )(&other);
+  const int *conn2=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
+  int size2=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
+  if(size1!=size2)
+    return false;
+  bool ret=true;
+  for(int i=0;i<size1 && ret;i++)
+    {
+      ret=(conn1[i]==conn2[i]);
+    }
+  return ret;
+}