src/MEDMEM/MEDMEM_Connectivity.cxx

   1 // Copyright (C) 2005  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 // This library is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU Lesser General Public
   6 // License as published by the Free Software Foundation; either
   7 // version 2.1 of the License.
   8 //
   9 // This library is distributed in the hope that it will be useful
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 // Lesser General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU Lesser General Public
  15 // License along with this library; if not, write to the Free Software
  16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 //
  20 #include "MEDMEM_Connectivity.hxx"
  21 #include "MEDMEM_Family.hxx"
  22 #include "MEDMEM_Group.hxx"
  23 #include "MEDMEM_CellModel.hxx"
  24
  25 #include "MEDMEM_SkyLineArray.hxx"
  26 #include "MEDMEM_ModulusArray.hxx"
  27
  28 #include "MEDMEM_STRING.hxx"
  29 #include <iomanip>
  30
  31 using namespace std;
  32 using namespace MEDMEM;
  33 using namespace MED_EN;
  34
  35 // Enlarge the vector if necessary to insert the element
  36 static inline void insert_vector(vector<int> &Vect, int Indice, int Element)
  37 {
  38   if (Indice >= Vect.capacity())
  39     Vect.reserve(Indice + 1000);
  40
  41   if (Indice >= Vect.size())
  42     Vect.resize(Indice+1);
  43
  44   Vect[Indice] = Element;
  45 }
  46
  47 /*!
  48    Default Constructor. /n
  49    Default for Entity is MED_CELL and type of Connectivity is MED_NODAL */
  50 //--------------------------------------------------------------//
  51 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity /* =MED_CELL */) :
  52 //--------------------------------------------------------------//
  53                                 _entity(Entity),
  54                                 _typeConnectivity(MED_NODAL),
  55                                 _numberOfTypes(0),
  56                                 _geometricTypes((medGeometryElement*)NULL),
  57                                 _type((CELLMODEL*)NULL),
  58                                 _entityDimension(0),
  59                                 _numberOfNodes(0),
  60                                 _count((int*)NULL),
  61                                 _nodal((MEDSKYLINEARRAY*)NULL),
  62                                 _polygonsNodal((MEDSKYLINEARRAY*)NULL),
  63                                 _polyhedronNodal((POLYHEDRONARRAY*)NULL),
  64                                 _descending((MEDSKYLINEARRAY*)NULL),
  65                                 _polygonsDescending((MEDSKYLINEARRAY*)NULL),
  66                                 _polyhedronDescending((MEDSKYLINEARRAY*)NULL),
  67                                 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
  68                                 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
  69                                 _neighbourhood((MEDSKYLINEARRAY*)NULL),
  70                                 _constituent((CONNECTIVITY*)NULL)
  71 {
  72    MESSAGE("CONNECTIVITY(medEntityMesh Entity=MED_CELL)");
  73 }
  74
  75 /*!
  76    Constructor. /n
  77    Default for Entity is MED_CELL */
  78 //------------------------------------------------------------------------------//
  79 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity /* =MED_CELL */):
  80 //------------------------------------------------------------------------------//
  81                                 _entity(Entity),
  82                                 _typeConnectivity(MED_NODAL),
  83                                 _numberOfTypes(numberOfTypes),
  84                                 _entityDimension(0),
  85                                 _numberOfNodes(0),
  86                                 _nodal((MEDSKYLINEARRAY*)NULL),
  87                                 _polygonsNodal((MEDSKYLINEARRAY*)NULL),
  88                                 _polyhedronNodal((POLYHEDRONARRAY*)NULL),
  89                                 _descending((MEDSKYLINEARRAY*)NULL),
  90                                 _polygonsDescending((MEDSKYLINEARRAY*)NULL),
  91                                 _polyhedronDescending((MEDSKYLINEARRAY*)NULL),
  92                                 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
  93                                 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
  94                                 _neighbourhood((MEDSKYLINEARRAY*)NULL),
  95                                 _constituent((CONNECTIVITY*)NULL)
  96 {
  97   MESSAGE("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)");
  98   _geometricTypes = new medGeometryElement[numberOfTypes];
  99   _type = new CELLMODEL[numberOfTypes];
 100   _count = new int[numberOfTypes+1];
 101 }
 102
 103 /*!
 104   Copy Constructor.
 105 */
 106 //------------------------------------------------------//
 107 CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
 108 //----------------------------------------------------//
 109                                 _entity           (m._entity),
 110                                 _typeConnectivity (m._typeConnectivity),
 111                                 _numberOfTypes    (m._numberOfTypes),
 112                                 _entityDimension  (m._entityDimension),
 113                                 _numberOfNodes    (m._numberOfNodes)
 114 {
 115  if (m._geometricTypes != NULL)
 116  {
 117     _geometricTypes = new medGeometryElement[_numberOfTypes];
 118     memcpy(_geometricTypes,m._geometricTypes,_numberOfTypes*sizeof(medGeometryElement));
 119  }
 120  else
 121     _geometricTypes = (medGeometryElement *) NULL;
 122
 123  if (m._type != NULL)
 124  {
 125     _type = new CELLMODEL[_numberOfTypes];
 126     for (int i=0;i<_numberOfTypes;i++)
 127         _type[i] = CELLMODEL(m._type[i]);
 128  }
 129  else
 130     _type = (CELLMODEL *) NULL;
 131
 132  if (m._count != NULL)
 133  {
 134       _count = new int[_numberOfTypes+1];
 135       memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(int));
 136  }
 137  else
 138     _count = (int *) NULL;
 139
 140  if (m._nodal != NULL)
 141     _nodal = new MEDSKYLINEARRAY(* m._nodal);
 142  else
 143     _nodal = (MEDSKYLINEARRAY *) NULL;
 144
 145  if (m._polygonsNodal != NULL)
 146    _polygonsNodal = new MEDSKYLINEARRAY(* m._polygonsNodal);
 147  else
 148    _polygonsNodal = (MEDSKYLINEARRAY *) NULL;
 149
 150  if (m._polyhedronNodal != NULL)
 151    _polyhedronNodal = new POLYHEDRONARRAY(* m._polyhedronNodal);
 152  else
 153    _polyhedronNodal = (POLYHEDRONARRAY *) NULL;
 154
 155  if (m._descending != NULL)
 156     _descending = new MEDSKYLINEARRAY(* m._descending);
 157  else
 158     _descending = (MEDSKYLINEARRAY *) NULL;
 159
 160  if (m._polygonsDescending != NULL)
 161    _polygonsDescending = new MEDSKYLINEARRAY(* m._polygonsDescending);
 162  else
 163    _polygonsDescending = (MEDSKYLINEARRAY *) NULL;
 164
 165  if (m._polyhedronDescending != NULL)
 166    _polyhedronDescending = new MEDSKYLINEARRAY(* m._polyhedronDescending);
 167  else
 168    _polyhedronDescending = (MEDSKYLINEARRAY *) NULL;
 169
 170  if (m._reverseNodalConnectivity != NULL)
 171     _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
 172  else
 173     _reverseNodalConnectivity = (MEDSKYLINEARRAY *) NULL;
 174
 175  if (m._reverseDescendingConnectivity != NULL)
 176     _reverseDescendingConnectivity = new MEDSKYLINEARRAY(* m._reverseDescendingConnectivity);
 177  else
 178     _reverseDescendingConnectivity = (MEDSKYLINEARRAY *) NULL;
 179
 180  if (m._neighbourhood != NULL)
 181     _neighbourhood = new MEDSKYLINEARRAY(* m._neighbourhood);
 182  else
 183     _neighbourhood = (MEDSKYLINEARRAY *) NULL;
 184
 185  if (m._constituent != NULL)
 186     _constituent = new CONNECTIVITY(* m._constituent);
 187  else
 188     _constituent = (CONNECTIVITY *) NULL;
 189 }
 190
 191 /*!
 192    Destructor./n
 193    desallocates existing pointers */
 194 //----------------------------//
 195 CONNECTIVITY::~CONNECTIVITY()
 196 //----------------------------//
 197 {
 198   MESSAGE("Destructeur de CONNECTIVITY()");
 199
 200   if (_geometricTypes != NULL)
 201      delete [] _geometricTypes;
 202   if (_type != NULL)
 203      delete [] _type;
 204   if (_count != NULL)
 205      delete [] _count;
 206   if (_nodal != NULL)
 207      delete _nodal;
 208   if (_polygonsNodal != NULL)
 209      delete _polygonsNodal;
 210   if (_polyhedronNodal != NULL)
 211      delete _polyhedronNodal;
 212   if (_descending != NULL)
 213      delete _descending;
 214   if (_polygonsDescending != NULL)
 215      delete _polygonsDescending;
 216   if (_polyhedronDescending != NULL)
 217     delete _polyhedronDescending;
 218   if (_reverseNodalConnectivity != NULL)
 219      delete _reverseNodalConnectivity;
 220   if (_reverseDescendingConnectivity != NULL)
 221      delete _reverseDescendingConnectivity;
 222   if (_constituent != NULL)
 223      delete _constituent;
 224 }
 225
 226 /*!
 227    set _constituent to Constituent
 228    be aware desallocation of _constituent is done by CONNECTIVITY:~CONNECTIVITY
 229    throws an exception if Constituent = MED_CELL
 230    A DOCUMENTER
 231     */
 232 //----------------------------------------------------------//
 233 void CONNECTIVITY::setConstituent(CONNECTIVITY * Constituent)
 234                     throw (MEDEXCEPTION)
 235 //----------------------------------------------------------//
 236 {
 237   medEntityMesh Entity = Constituent->getEntity();
 238   if (Entity == MED_CELL)
 239     throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : could not set constituent on MED_CELL !"));
 240
 241   if ((Entity == MED_EDGE)&(_entityDimension == 3))
 242   {
 243     if (_constituent == NULL)
 244       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : Entity not found !"));
 245     _constituent->setConstituent(Constituent);
 246   }
 247   else
 248     _constituent = Constituent;
 249 }
 250
 251 /*! Duplicated Types array in CONNECTIVITY object. */
 252 //--------------------------------------------------------------------//
 253 void CONNECTIVITY::setGeometricTypes(const medGeometryElement * Types,
 254                                      const medEntityMesh Entity)
 255                                      throw (MEDEXCEPTION)
 256 //--------------------------------------------------------------------//
 257 {
 258   if (Entity == _entity)
 259     for (int i=0; i<_numberOfTypes; i++)
 260     {
 261       _geometricTypes[i] = Types[i];
 262       _type[i] = CELLMODEL(_geometricTypes[i]);
 263     }
 264   else
 265     {
 266     if (_constituent == NULL)
 267         throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setGeometricTypes : Entity not found !"));
 268     _constituent->setGeometricTypes(Types,Entity);
 269   }
 270 }
 271
 272 /*! A DOCUMENTER */
 273 //--------------------------------------------------------------------//
 274 void CONNECTIVITY::setCount(const int * Count, const medEntityMesh Entity)
 275                             throw (MEDEXCEPTION)
 276 //--------------------------------------------------------------------//
 277 {
 278   if (Entity == _entity)
 279   {
 280     int * index = new int[Count[_numberOfTypes]];
 281     index[0]=1;
 282     _count[0]=1;
 283     for (int i=0; i<_numberOfTypes; i++) {
 284       _count[i+1] = Count[i+1];
 285       int NumberOfNodesPerElement = _type[i].getNumberOfNodes() ;
 286       for (int j=_count[i]; j<_count[i+1]; j++)
 287         index[j] = index[j-1]+NumberOfNodesPerElement;
 288     }
 289     // allocate _nodal
 290     if (_nodal != NULL) delete _nodal;
 291     _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,index[_count[_numberOfTypes]-1]-1);
 292     _nodal->setIndex(index);
 293     delete [] index;
 294   }
 295   else
 296   {
 297     if (_constituent == NULL)
 298       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setCount : Entity not found !"));
 299     _constituent->setCount(Count,Entity);
 300   }
 301 }
 302
 303 //--------------------------------------------------------------------//
 304 void CONNECTIVITY::setNodal(const int * Connectivity,
 305                             const medEntityMesh Entity,
 306                             const medGeometryElement Type)
 307                             throw (MEDEXCEPTION)
 308 //--------------------------------------------------------------------//
 309 {
 310   if (_entity == Entity)
 311   {
 312     // find geometric type
 313     bool find = false;
 314     for (int i=0; i<_numberOfTypes; i++)
 315       if (_geometricTypes[i] == Type) {
 316         find = true;
 317         int NumberOfNodePerElements = _type[i].getNumberOfNodes() ;
 318         //_nodal->setI(i+1,Connectivity);
 319         for( int j=_count[i];j<_count[i+1]; j++)
 320           _nodal->setI(j,Connectivity+(j-_count[i])*NumberOfNodePerElements);
 321       }
 322     if (!find)
 323       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : geometric type not found"));
 324   } else
 325   {
 326     if (_constituent == NULL)
 327       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : Entity not found !"));
 328       _constituent->setNodal(Connectivity,Entity,Type);
 329   }
 330 }
 331
 332
 333 //--------------------------------------------------------------------//
 334 void CONNECTIVITY::setPolygonsConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, const int* PolygonsConnectivity, const int* PolygonsConnectivityIndex, int ConnectivitySize, int NumberOfPolygons)
 335 //--------------------------------------------------------------------//
 336 {
 337   const char* LOC = "CONNECTIVITY::setPolygonsConnectivity";
 338   BEGIN_OF(LOC);
 339
 340   if (_entity == Entity)
 341     {
 342       MEDSKYLINEARRAY* Connectivity = new MEDSKYLINEARRAY(NumberOfPolygons,ConnectivitySize,PolygonsConnectivityIndex,PolygonsConnectivity);
 343
 344       if (ConnectivityType == MED_NODAL)
 345         {
 346           if (_polygonsNodal != (MEDSKYLINEARRAY*) NULL)
 347             delete _polygonsNodal;
 348           _polygonsNodal = Connectivity;
 349         }
 350       else
 351         {
 352           if (_typeConnectivity != MED_DESCENDING)
 353             _typeConnectivity = MED_DESCENDING; //by default it is set to MED_NODAL
 354           if (_polygonsDescending != (MEDSKYLINEARRAY*) NULL)
 355             delete _polygonsDescending;
 356           _polygonsDescending = Connectivity;
 357         }
 358     }
 359   else
 360     {
 361       if (_constituent == (CONNECTIVITY*) NULL)
 362         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not found !"));
 363       _constituent->setPolygonsConnectivity(ConnectivityType, Entity, PolygonsConnectivity, PolygonsConnectivityIndex, ConnectivitySize, NumberOfPolygons);
 364     }
 365 }
 366
 367
 368 //--------------------------------------------------------------------//
 369 void CONNECTIVITY::setPolyhedronConnectivity(medConnectivity ConnectivityType, const int* PolyhedronConnectivity, const int* PolyhedronIndex, int ConnectivitySize, int NumberOfPolyhedron, const int* PolyhedronFacesIndex /* =(int*)NULL */, int NumberOfFaces /* =0 */)
 370 //--------------------------------------------------------------------//
 371 {
 372   const char* LOC = "CONNECTIVITY::setPolyhedronConnectivity";
 373   BEGIN_OF(LOC);
 374
 375   if (_entity == MED_CELL)
 376     {
 377       if (ConnectivityType == MED_NODAL)
 378         {
 379           if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
 380             delete _polyhedronNodal;
 381           _polyhedronNodal = new POLYHEDRONARRAY(NumberOfPolyhedron,NumberOfFaces,ConnectivitySize);
 382 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000)
 383           int i ;
 384           MED_EN::med_int * tmp_PolyhedronIndex = new med_int[NumberOfPolyhedron+1] ;
 385           for ( i = 0 ; i < NumberOfPolyhedron+1 ; i++ )
 386              tmp_PolyhedronIndex[i] = PolyhedronIndex[i] ;
 387           _polyhedronNodal->setPolyhedronIndex(tmp_PolyhedronIndex);
 388           delete [] tmp_PolyhedronIndex ;
 389           MED_EN::med_int * tmp_PolyhedronFacesIndex = new med_int[NumberOfFaces+1] ;
 390           for ( i = 0 ; i < NumberOfFaces+1 ; i++ )
 391              tmp_PolyhedronFacesIndex[i] = PolyhedronFacesIndex[i] ;
 392           _polyhedronNodal->setFacesIndex(tmp_PolyhedronFacesIndex);
 393           delete [] tmp_PolyhedronFacesIndex ;
 394           MED_EN::med_int * tmp_PolyhedronConnectivity = new med_int[ConnectivitySize] ;
 395           for ( i = 0 ; i < ConnectivitySize ; i++ )
 396              tmp_PolyhedronConnectivity[i] = PolyhedronConnectivity[i] ;
 397           _polyhedronNodal->setNodes(tmp_PolyhedronConnectivity);
 398           delete [] tmp_PolyhedronConnectivity ;
 399 #else
 400           _polyhedronNodal->setPolyhedronIndex(PolyhedronIndex);
 401           _polyhedronNodal->setFacesIndex(PolyhedronFacesIndex);
 402           _polyhedronNodal->setNodes(PolyhedronConnectivity);
 403 #endif
 404         }
 405       else
 406         {
 407           if (_typeConnectivity != MED_DESCENDING)
 408             _typeConnectivity = MED_DESCENDING; //by default it is set to MED_NODAL
 409           if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
 410             delete _polyhedronDescending;
 411           _polyhedronDescending = new MEDSKYLINEARRAY(NumberOfPolyhedron,ConnectivitySize,PolyhedronIndex,PolyhedronConnectivity);
 412         }
 413     }
 414   else
 415     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : _entity must be MED_CELL to set polyhedron !"));
 416 }
 417
 418
 419 /*! A DOCUMENTER */
 420 //------------------------------------------------------------------------------------------//
 421 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
 422 //------------------------------------------------------------------------------------------//
 423 {
 424   MESSAGE("CONNECTIVITY::calculateConnectivity");
 425
 426   // a temporary limitation due to calculteDescendingConnectivity function !!!
 427   if ((_entityDimension==3) & (Entity==MED_EDGE))
 428     throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
 429
 430   if (Entity==_entity)
 431     if (ConnectivityType==MED_NODAL)
 432       calculateNodalConnectivity();
 433     else
 434       if (Entity==MED_CELL)
 435         calculateDescendingConnectivity();
 436       else
 437         throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
 438   if (Entity!=_entity) {
 439     calculateDescendingConnectivity();
 440     if (_entityDimension == 2 || _entityDimension == 3)
 441       _constituent->calculateConnectivity(ConnectivityType,Entity);
 442   }
 443 }
 444
 445 /*!  Give, in full or no interlace mode (for nodal connectivity),
 446      descending or nodal connectivity.
 447
 448       You must give a %medEntityMesh (ie:MED_EDGE)
 449       and a %medGeometryElement (ie:MED_SEG3).  */
 450
 451 //------------------------------------------------------------//
 452 void CONNECTIVITY::updateFamily(const vector<FAMILY*>& myFamilies)
 453 //------------------------------------------------------------//
 454 {
 455   const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
 456   int numberOfFamilies = myFamilies.size();
 457   if (numberOfFamilies == 0 || _constituent == NULL)
 458     return;
 459   // does we do an update ?
 460   if ((_constituent != NULL) && (_descending != NULL))
 461     return;
 462   if (myFamilies[0]->getEntity() != _constituent->getEntity())
 463     return;
 464   CONNECTIVITY * oldConstituent = _constituent;
 465   _constituent = (CONNECTIVITY *)NULL;
 466   if (oldConstituent->_nodal==NULL)
 467     throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
 468
 469   //Loc vars defined to treat polygons exactly the same as classic types. Not nice but necessary.
 470   int oldNumberOfFaceTab[2];
 471   const int * oldConstituentValueTab[2];
 472   const int * oldConstituentIndexTab[2];
 473   int * renumberingFromOldToNewTab[2];//Final mapping array between old numbers and new numbers.;
 474
 475   int oldNumberOfFace = oldConstituent->_nodal->getNumberOf(); oldNumberOfFaceTab[0]=oldNumberOfFace;
 476   const int * oldConstituentValue = oldConstituent->_nodal->getValue(); oldConstituentValueTab[0]=oldConstituentValue;
 477   const int * oldConstituentIndex = oldConstituent->_nodal->getIndex(); oldConstituentIndexTab[0]=oldConstituentIndex;
 478   int * renumberingFromOldToNew= new int [oldNumberOfFace]; renumberingFromOldToNewTab[0]=renumberingFromOldToNew;
 479
 480   int oldNumberOfFacePoly = oldConstituent->getNumberOfPolygons();
 481   const int * oldConstituentValuePoly=0;
 482   const int * oldConstituentIndexPoly=0;
 483   int * renumberingFromOldToNewPoly=0;
 484
 485   int nbOfTurnInGlobalLoop=1;//Defined to know if a second search on polygons is needed.
 486   if(oldNumberOfFacePoly>0)
 487     {
 488       oldNumberOfFaceTab[1]=oldNumberOfFacePoly;
 489       oldConstituentValuePoly = oldConstituent->_polygonsNodal->getValue(); oldConstituentValueTab[1]=oldConstituentValuePoly;
 490       oldConstituentIndexPoly = oldConstituent->_polygonsNodal->getIndex(); oldConstituentIndexTab[1]=oldConstituentIndexPoly;
 491       renumberingFromOldToNewPoly=new int[oldNumberOfFacePoly]; renumberingFromOldToNewTab[1]=renumberingFromOldToNewPoly;
 492       nbOfTurnInGlobalLoop++;
 493     }
 494
 495   calculateDescendingConnectivity();//perform calculation of descending connectivity to update all connectivities after having taken into account of direction of faces...
 496   _constituent->calculateReverseNodalConnectivity();//getting d-1 nodal connectivity to find get new face numbers from nodes numbers...
 497
 498   const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity(); //Common to polygons and classic geometric types
 499   const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex(); //Common to polygons and classic geometric types
 500
 501   for(int loop=0;loop<nbOfTurnInGlobalLoop;loop++)
 502     {
 503       int oldNumberOfFaceLoop=oldNumberOfFaceTab[loop];
 504       const int * oldConstituentValueLoop=oldConstituentValueTab[loop];
 505       const int * oldConstituentIndexLoop= oldConstituentIndexTab[loop];
 506       int * renumberingFromOldToNewLoop=renumberingFromOldToNewTab[loop];
 507       for(int iOldFace=0;iOldFace<oldNumberOfFaceLoop;iOldFace++)
 508         {
 509           const int *nodesOfCurrentFaceOld=oldConstituentValueLoop+oldConstituentIndexLoop[iOldFace]-1;
 510           int nbOfNodesOfCurrentFaceOld=oldConstituentIndexLoop[iOldFace+1]-oldConstituentIndexLoop[iOldFace];
 511
 512           //retrieving new number of polygon face...
 513           int sizeOfNewFaceNb=reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1];
 514           int *newFaceNb=new int[sizeOfNewFaceNb];
 515           memcpy(newFaceNb,reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1]-1,sizeOfNewFaceNb*sizeof(int));//all faces candidates are in newFaceNb !!!
 516           for(int curNode=1;curNode<nbOfNodesOfCurrentFaceOld && sizeOfNewFaceNb>1;curNode++)
 517             {
 518               const int *newFacesNbForCurNode=reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1]-1;
 519               int sizeOfNewFacesNbForCurNode=reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1];
 520               for(int i=0;i<sizeOfNewFaceNb;)
 521                 {
 522                   bool found=false;
 523                   for(int j=0;j<sizeOfNewFacesNbForCurNode && !found;j++)
 524                     {
 525                       if(newFacesNbForCurNode[j]==newFaceNb[i])
 526                         found=true;
 527                     }
 528                   if(found)
 529                     i++;
 530                   else
 531                     newFaceNb[i]=newFaceNb[--sizeOfNewFaceNb];
 532                 }
 533             }
 534           if(sizeOfNewFaceNb!=1)
 535             throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"unexisting face specified!!!"));
 536           renumberingFromOldToNewLoop[iOldFace]=newFaceNb[0];
 537           delete [] newFaceNb;
 538           //end of retrieving new number of polygon face...
 539           int nbOfNodesOfCurrentFaceNew;
 540           const int *nodesOfCurrentFaceNew=_constituent->getConnectivityOfAnElementWithPoly(MED_NODAL,_constituent->getEntity(),
 541                                                                                             renumberingFromOldToNewLoop[iOldFace],nbOfNodesOfCurrentFaceNew);
 542           MEDMODULUSARRAY modulusArrayOld(nbOfNodesOfCurrentFaceOld,nodesOfCurrentFaceOld);
 543           MEDMODULUSARRAY modulusArrayNew(nbOfNodesOfCurrentFaceNew,nodesOfCurrentFaceNew);
 544           int retCompareNewOld=modulusArrayNew.compare(modulusArrayOld);
 545           if(retCompareNewOld==0)
 546             throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"Uncompatible given user face with calculated existing faces"));
 547           if(retCompareNewOld==-1)
 548             invertConnectivityForAFace(renumberingFromOldToNewLoop[iOldFace],nodesOfCurrentFaceOld,loop==1);
 549         }
 550     }
 551   // Updating the Family
 552   for(vector<FAMILY*>::const_iterator iter=myFamilies.begin();iter!=myFamilies.end();iter++)
 553     (*iter)->changeElementsNbs(_constituent->getEntity(),renumberingFromOldToNew,oldNumberOfFace,renumberingFromOldToNewPoly);
 554   delete oldConstituent ;
 555   delete [] renumberingFromOldToNew;
 556   if(oldNumberOfFacePoly>0)
 557     delete [] renumberingFromOldToNewPoly;
 558   return;
 559 }
 560
 561 //------------------------------------------------------------------------------------------------------------------//
 562 const int * MEDMEM::CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
 563 //------------------------------------------------------------------------------------------------------------------//
 564 {
 565   const char * LOC = "CONNECTIVITY::getConnectivity";
 566   //  BEGIN_OF(LOC);
 567
 568   MEDSKYLINEARRAY * Connectivity;
 569   if (Entity==_entity) {
 570
 571     if (ConnectivityType==MED_NODAL)
 572       {
 573         calculateNodalConnectivity();
 574         Connectivity=_nodal;
 575       }
 576     else
 577       {
 578         calculateDescendingConnectivity();
 579         Connectivity=_descending;
 580       }
 581
 582     if (Connectivity!=NULL)
 583       if (Type==MED_ALL_ELEMENTS)
 584         return Connectivity->getValue();
 585       else {
 586         for (int i=0; i<_numberOfTypes; i++)
 587           if (_geometricTypes[i]==Type)
 588             //return Connectivity->getI(i+1);
 589             return Connectivity->getI(_count[i]);
 590         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
 591       }
 592     else
 593       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 594   } else
 595     if (_constituent != NULL)
 596       return _constituent->getConnectivity(ConnectivityType,Entity,Type);
 597
 598   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 599 }
 600
 601 //------------------------------------------------------------------------------------------------------------------//
 602 int CONNECTIVITY::getConnectivityLength(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
 603 //------------------------------------------------------------------------------------------------------------------//
 604 {
 605   const char * LOC = "CONNECTIVITY::getConnectivity";
 606   BEGIN_OF(LOC);
 607
 608   MEDSKYLINEARRAY * Connectivity;
 609   if (Entity==_entity) {
 610
 611     if (ConnectivityType==MED_NODAL)
 612       {
 613         calculateNodalConnectivity();
 614         Connectivity=_nodal;
 615       }
 616     else
 617       {
 618         calculateDescendingConnectivity();
 619         Connectivity=_descending;
 620       }
 621
 622     if (Connectivity!=NULL)
 623       if (Type==MED_ALL_ELEMENTS)
 624         return Connectivity->getLength();
 625       else {
 626         for (int i=0; i<_numberOfTypes; i++)
 627           if (_geometricTypes[i]==Type)
 628             return Connectivity->getNumberOfI(_count[i]);
 629         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
 630       }
 631     else
 632       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 633   }
 634   else
 635     if (_constituent != NULL)
 636       return _constituent->getConnectivityLength(ConnectivityType,Entity,Type);
 637
 638   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 639 }
 640
 641 /*!  Give morse index array to use with
 642      getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
 643
 644       Each value give start index for corresponding entity in connectivity array.
 645
 646       Example : i-th element, j-th node of it :
 647       - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
 648       - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
 649 //-----------------------------------------------------------------------------------------------//
 650 const int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
 651 //-----------------------------------------------------------------------------------------------//
 652 {
 653   const char * LOC = "CONNECTIVITY::getConnectivityIndex";
 654
 655   MEDSKYLINEARRAY * Connectivity;
 656   if (Entity==_entity) {
 657
 658     if (ConnectivityType==MED_NODAL)
 659       Connectivity=_nodal;
 660     else
 661       Connectivity=_descending;
 662
 663     if (Connectivity!=NULL)
 664       return Connectivity->getIndex();
 665     else
 666       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 667
 668   } else
 669     if (_constituent != NULL)
 670       return _constituent->getConnectivityIndex(ConnectivityType,Entity);
 671
 672   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 673 }
 674
 675
 676 //-------------------------------------------------------------//
 677 const int* CONNECTIVITY::getPolygonsConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
 678 //-------------------------------------------------------------//
 679 {
 680   const char* LOC = "CONNECTIVITY::getPolygonsConnectivity";
 681   BEGIN_OF(LOC);
 682
 683   MEDSKYLINEARRAY* Connectivity;
 684   if (Entity == _entity)
 685     {
 686       if (ConnectivityType == MED_NODAL)
 687         {
 688           calculateNodalConnectivity();
 689           Connectivity = _polygonsNodal;
 690         }
 691       else
 692         {
 693           calculateDescendingConnectivity();
 694           Connectivity = _polygonsDescending;
 695         }
 696       if (Connectivity != NULL)
 697         return Connectivity->getValue();
 698       else
 699         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polygons Connectivity not defined !"));
 700     }
 701   else
 702     if (_constituent != NULL)
 703       return _constituent->getPolygonsConnectivity(ConnectivityType, Entity);
 704   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 705 }
 706
 707
 708 //-------------------------------------------------------------//
 709 const int* CONNECTIVITY::getPolygonsConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
 710 //-------------------------------------------------------------//
 711 {
 712   const char* LOC = "CONNECTIVITY::getPolygonsConnectivityIndex";
 713   BEGIN_OF(LOC);
 714
 715   MEDSKYLINEARRAY* Connectivity;
 716   if (Entity == _entity)
 717     {
 718       if (ConnectivityType == MED_NODAL)
 719         {
 720           //      calculateNodalConnectivity();
 721           Connectivity = _polygonsNodal;
 722         }
 723       else
 724         {
 725           //      calculateDescendingConnectivity();
 726           Connectivity = _polygonsDescending;
 727         }
 728       if (Connectivity != NULL)
 729         return Connectivity->getIndex();
 730       else
 731         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polygons Connectivity not defined !"));
 732     }
 733   else
 734     if (_constituent != NULL)
 735       return _constituent->getPolygonsConnectivityIndex(ConnectivityType, Entity);
 736   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 737 }
 738
 739
 740 /*! We suppose in this method that nodal and descending connectivities
 741   are coherent.*/
 742 //-------------------------------------------------------------//
 743 int CONNECTIVITY::getNumberOfPolygons(MED_EN::medEntityMesh Entity) const
 744 //-------------------------------------------------------------//
 745 {
 746   if(Entity==MED_ALL_ENTITIES || Entity==_entity )
 747     {
 748       if (_polygonsNodal != (MEDSKYLINEARRAY*) NULL)
 749         return _polygonsNodal->getNumberOf();
 750       else if (_polygonsDescending != (MEDSKYLINEARRAY*) NULL)
 751         return _polygonsDescending->getNumberOf();
 752       else
 753         return 0;
 754     }
 755   else
 756     {
 757       if (_constituent == (CONNECTIVITY*) NULL)
 758         throw MEDEXCEPTION("getNumberOfPolygons : Entity not found !");
 759       return _constituent->getNumberOfPolygons(Entity);
 760     }
 761 }
 762
 763
 764 //--------------------------------------------------------------//
 765 const int* CONNECTIVITY::getPolyhedronConnectivity(medConnectivity ConnectivityType) const
 766 //--------------------------------------------------------------//
 767 {
 768   const char* LOC = "CONNECTIVITY::getPolyhedronConnectivity";
 769   BEGIN_OF(LOC);
 770
 771   if (_entity == MED_CELL) //polyhedron can only be MED_CELL
 772     {
 773       if (ConnectivityType == MED_NODAL)
 774         {
 775           ((CONNECTIVITY *)(this))->calculateNodalConnectivity();
 776           if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL) {
 777 //CCRT
 778 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000)
 779             int i ;
 780             const MED_EN::med_int * tmp_PolyhedronConnectivity = _polyhedronNodal->getNodes();
 781             int * PolyhedronConnectivity = new int[_polyhedronNodal->getNumberOfNodes()] ;
 782             for ( i = 0 ; i < _polyhedronNodal->getNumberOfNodes() ; i++ )
 783                PolyhedronConnectivity[i] = tmp_PolyhedronConnectivity[i] ;
 784 //CCRT : return of a copy of PolyhedronConnectivity
 785             return PolyhedronConnectivity ;
 786 #else
 787             return _polyhedronNodal->getNodes();
 788 #endif
 789           }
 790           else
 791             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Nodal Connectivity not defined !"));
 792         }
 793       else
 794         {
 795           ((CONNECTIVITY *)(this))->calculateDescendingConnectivity();
 796           if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL) {
 797             return _polyhedronDescending->getValue();
 798           }
 799           else
 800             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Descending Connectivity not defined !"));
 801         }
 802     }
 803   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
 804 }
 805
 806
 807 //---------------------------------------------------------------//
 808 const int* CONNECTIVITY::getPolyhedronFacesIndex() const
 809 //---------------------------------------------------------------//
 810 {
 811   const char* LOC = "CONNECTIVITY::getPolyhedronFacesIndex";
 812   BEGIN_OF(LOC);
 813
 814   if (_entity == MED_CELL) //polyhedron can only be MED_CELL
 815     {
 816       //      calculateNodalConnectivity();
 817       if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL) {
 818 //CCRT
 819 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000)
 820             int i ;
 821             const MED_EN::med_int * tmp_PolyhedronFacesIndex = _polyhedronNodal->getFacesIndex();
 822             int * PolyhedronFacesIndex = new int[_polyhedronNodal->getNumberOfFaces()+1] ;
 823             for ( i = 0 ; i < _polyhedronNodal->getNumberOfFaces()+1 ; i++ )
 824                PolyhedronFacesIndex[i] = tmp_PolyhedronFacesIndex[i] ;
 825 //CCRT : return of a copy of PolyhedronFacesIndex
 826             return PolyhedronFacesIndex ;
 827 #else
 828         return _polyhedronNodal->getFacesIndex();
 829 #endif
 830       }
 831       else
 832         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No Polyhedron in that Connectivity !"));
 833     }
 834   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
 835 }
 836
 837
 838 //--------------------------------------------------------------//
 839 const int* CONNECTIVITY::getPolyhedronIndex(medConnectivity ConnectivityType) const
 840 //--------------------------------------------------------------//
 841 {
 842   const char* LOC = "CONNECTIVITY::getPolyhedronIndex";
 843   BEGIN_OF(LOC);
 844
 845   if (_entity == MED_CELL) //polyhedron can only be MED_CELL
 846     {
 847       if (ConnectivityType == MED_NODAL)
 848         {
 849           //      calculateNodalConnectivity();
 850           if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL) {
 851 //CCRT
 852 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000)
 853             int i ;
 854             const MED_EN::med_int * tmp_PolyhedronIndex = _polyhedronNodal->getPolyhedronIndex();
 855             int * PolyhedronIndex = new int[_polyhedronNodal->getNumberOfPolyhedron()+1] ;
 856             for ( i = 0 ; i < _polyhedronNodal->getNumberOfPolyhedron()+1 ; i++ )
 857                PolyhedronIndex[i] = tmp_PolyhedronIndex[i] ;
 858 //CCRT : return of a copy of PolyhedronIndex
 859             return PolyhedronIndex ;
 860 #else
 861             return _polyhedronNodal->getPolyhedronIndex();
 862 #endif
 863           }
 864           else
 865             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Nodal Connectivity not defined !"));
 866         }
 867       else
 868         {
 869           //      calculateDescendingConnectivity();
 870           if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
 871             return _polyhedronDescending->getIndex();
 872           else
 873             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Descending Connectivity not defined !"));
 874         }
 875     }
 876   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
 877 }
 878
 879
 880 /*! We suppose in this method that nodal and descending connectivities
 881   are coherent.*/
 882 //-------------------------------------------------------------//
 883 int CONNECTIVITY::getNumberOfPolyhedronFaces() const
 884 //-------------------------------------------------------------//
 885 {
 886   //  if (_polyhedronNodal == (POLYHEDRONARRAY*) NULL)
 887   //    calculateNodalConnectivity();
 888   if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
 889     return _polyhedronNodal->getNumberOfFaces();
 890   else
 891     return 0;
 892 }
 893
 894
 895 /*! We suppose in this method that nodal and descending connectivities
 896   are coherent.*/
 897 //--------------------------------------------------------------//
 898 int CONNECTIVITY::getNumberOfPolyhedron() const
 899 //--------------------------------------------------------------//
 900 {
 901   if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
 902     return _polyhedronNodal->getNumberOfPolyhedron();
 903   else if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
 904     return _polyhedronDescending->getNumberOf();
 905   else
 906     return 0;
 907 }
 908
 909
 910 /*! A DOCUMENTER */
 911 //--------------------------------------------------------------//
 912 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
 913 //--------------------------------------------------------------//
 914 {
 915   const char * LOC = "CONNECTIVITY::getType";
 916   BEGIN_OF(LOC);
 917
 918   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 919     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !"));
 920   for (int i=0; i<_numberOfTypes; i++)
 921     if (_geometricTypes[i]==Type)
 922       return _type[i];
 923   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 924 }
 925
 926 /*!  Returns the number of elements of type %medGeometryElement.
 927      Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
 928 //------------------------------------------------------------------------//
 929 int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
 930 //------------------------------------------------------------------------//
 931 {
 932   const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
 933   BEGIN_OF(LOC);
 934
 935   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 936     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
 937   for (int i=0; i<_numberOfTypes; i++)
 938     if (_geometricTypes[i]==Type)
 939       return _type[i].getNumberOfNodes();
 940   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 941 }
 942
 943 /*!  Returns the number of geometric sub cells of %medGeometryElement type.
 944 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
 945 //------------------------------------------------------------------------//
 946 int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
 947 //------------------------------------------------------------------------//
 948 {
 949   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 950     throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
 951   for (int i=0; i<_numberOfTypes; i++)
 952     if (_geometricTypes[i]==Type)
 953       return _type[i].getNumberOfConstituents(1);
 954   throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
 955 }
 956
 957 /*!  Returns the number of elements of type %medGeometryElement.
 958
 959      Note :
 960       - Implemented for MED_ALL_ELEMENTS
 961       - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
 962       - Not implemented for MED_NONE */
 963 //-----------------------------------------------------------------------------------//
 964 int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
 965 //-----------------------------------------------------------------------------------//
 966 {
 967   //const char * LOC = "CONNECTIVITY::getNumberOf";
 968   if (Entity==_entity) {
 969     if (Type==MED_NONE)
 970       return 0; // not defined !
 971     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
 972     if (!existConnectivity(MED_NODAL,Entity) && !existConnectivity(MED_DESCENDING,Entity)) return 0; //case with only polygons for example
 973     if (Type==MED_ALL_ELEMENTS)
 974       return _count[_numberOfTypes]-1;
 975     for (int i=0; i<_numberOfTypes; i++)
 976       if (_geometricTypes[i]==Type)
 977         return (_count[i+1] - _count[i]);
 978     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 979   } else
 980     if (_constituent != NULL)
 981       return _constituent->getNumberOf(Entity,Type);
 982
 983   return 0; // valid if they are nothing else !
 984   //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 985 }
 986
 987 /*! A DOCUMENTER */
 988 //--------------------------------------------------------------//
 989 const int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
 990                                 medGeometryElement Type)
 991 //--------------------------------------------------------------//
 992 {
 993   if (TypeConnectivity == MED_NODAL)
 994     {
 995           calculateNodalConnectivity();
 996           if (Type==MED_ALL_ELEMENTS)
 997             return _nodal->getValue();
 998           for (int i=0; i<_numberOfTypes; i++)
 999             if (_geometricTypes[i]==Type)
1000               return _nodal->getI(_count[i]);
1001     }
1002   else
1003     {
1004       calculateDescendingConnectivity();
1005       if (Type==MED_ALL_ELEMENTS)
1006         return _descending->getValue();
1007       for (int i=0; i<_numberOfTypes; i++)
1008         if (_geometricTypes[i]==Type)
1009           return _descending->getI(Type);
1010     }
1011   throw MEDEXCEPTION("Not found");
1012 }
1013
1014 /*! A DOCUMENTER */
1015 //---------------------------------------------------------------------//
1016 const int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
1017 //---------------------------------------------------------------------//
1018 {
1019   if (TypeConnectivity == MED_NODAL)
1020     {
1021       calculateNodalConnectivity();
1022       return _nodal->getIndex();
1023     }
1024   else
1025     {
1026       calculateDescendingConnectivity();
1027       return _descending->getIndex();
1028     }
1029 }
1030
1031 /*! Not yet implemented */
1032 //----------------------------------------------//
1033 const int* CONNECTIVITY:: getNeighbourhood() const
1034 //----------------------------------------------//
1035 {
1036   throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
1037 }
1038
1039 /*! Returns an array which contains, for each node, all cells
1040     arround it. */
1041 //-------------------------------------------------//
1042 const int* CONNECTIVITY::getReverseNodalConnectivity()
1043 //-------------------------------------------------//
1044 {
1045   calculateReverseNodalConnectivity();
1046   return _reverseNodalConnectivity->getValue();
1047 }
1048
1049 /*!  Give index array to use with getReverseConnectivity(MED_NODAL).
1050      It is unusefull with MED_DESCENDING mode, because we have allways two cells.  */
1051 //-------------------------------------------------------//
1052 const int* CONNECTIVITY::getReverseNodalConnectivityIndex()
1053 //-------------------------------------------------------//
1054 {
1055   calculateReverseNodalConnectivity();
1056   return _reverseNodalConnectivity->getIndex();
1057 }
1058
1059 /*! Returns an array which contains, for each face (or edge),
1060     the 2 cells of each side. First is cell which face normal is outgoing.
1061     arround it. */
1062 //------------------------------------------------------//
1063 const int* CONNECTIVITY::getReverseDescendingConnectivity()
1064 //------------------------------------------------------//
1065 {
1066   // it is in _constituent connectivity only if we are in MED_CELL
1067   // (we could not for instance calculate face-edge connectivity !)
1068   if (_entity!=MED_CELL)
1069     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
1070
1071   // we need descending connectivity
1072   calculateDescendingConnectivity();
1073   return _reverseDescendingConnectivity->getValue();
1074 }
1075
1076 /*! calculate the reverse descending Connectivity
1077     and returns the index  ( A DOCUMENTER MIEUX)*/
1078 //-----------------------------------------------------------//
1079 const int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
1080 //-----------------------------------------------------------//
1081 {
1082   // it is in _constituent connectivity only if we are in MED_CELL
1083   if (_entity!=MED_CELL)
1084     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
1085
1086   // we need descending connectivity
1087   calculateDescendingConnectivity();
1088   return _reverseDescendingConnectivity->getIndex();
1089 }
1090
1091 /*! A DOCUMENTER (et a finir ???) */
1092 //--------------------------------------------//
1093 void CONNECTIVITY::calculateNodalConnectivity()
1094 //--------------------------------------------//
1095 {
1096   if (_nodal==NULL && _polygonsNodal==NULL && _polyhedronNodal==NULL)
1097     {
1098       if (_descending==NULL && _polygonsDescending==NULL && _polyhedronDescending==NULL)
1099         throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
1100       // calculate _nodal, _polygonsNodal and _polyhedronNodal from _descending, _polygonsDescending and _polyhedronDescending
1101     }
1102 }
1103
1104 /*! If not yet done, calculate the nodal Connectivity
1105     and the reverse nodal Connectivity*/
1106 //---------------------------------------------------//
1107 void CONNECTIVITY::calculateReverseNodalConnectivity()
1108 //---------------------------------------------------//
1109 {
1110   const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
1111   BEGIN_OF(LOC);
1112
1113   SCRUTE(_nodal);
1114   SCRUTE(_reverseNodalConnectivity);
1115
1116   if (_nodal==NULL)
1117     calculateNodalConnectivity();
1118
1119   if(_reverseNodalConnectivity==NULL) {
1120
1121     int node_number = 0;
1122     vector <vector <int> > reverse_connectivity;
1123     reverse_connectivity.resize(_numberOfNodes+1);
1124
1125   // Treat all cells types
1126
1127     for (int j = 0; j < _numberOfTypes; j++)
1128       {
1129         // node number of the cell type
1130         node_number = _type[j].getNumberOfNodes();
1131         // treat all cells of a particular type
1132         for (int k = _count[j]; k < _count[j+1]; k++)
1133           // treat all nodes of the cell type
1134           for (int local_node_number = 1; local_node_number < node_number+1; local_node_number++)
1135             {
1136               int global_node = _nodal->getIJ(k,local_node_number);
1137               reverse_connectivity[global_node].push_back(k);
1138             }
1139       }
1140
1141     if(_polygonsNodal != NULL && _entityDimension==2)
1142       {
1143         int nbOfPolygons=_polygonsNodal->getNumberOf();
1144         int offset=getNumberOf(_entity,MED_ALL_ELEMENTS);
1145         const int *index=_polygonsNodal->getIndex();
1146         const int *value=_polygonsNodal->getValue();
1147         for (int local_polyg_number = 0; local_polyg_number < nbOfPolygons; local_polyg_number++)
1148           {
1149             for(int i=index[local_polyg_number];i<index[local_polyg_number+1];i++)
1150               reverse_connectivity[value[i-1]].push_back(offset+local_polyg_number+1);
1151           }
1152       }
1153
1154     if(_polyhedronNodal != NULL && _entityDimension==3)
1155       {
1156         int nbOfPolyhedra=_polyhedronNodal->getNumberOfPolyhedron();
1157         int offset=getNumberOf(_entity,MED_ALL_ELEMENTS);
1158         for (int local_polyh_number = 0; local_polyh_number < nbOfPolyhedra; local_polyh_number++)
1159             {
1160               int nbOfNodes;
1161               int global_polyh_number=offset+local_polyh_number+1;
1162               int *nodes=getNodesOfPolyhedron(global_polyh_number,nbOfNodes);
1163               for(int i=0;i<nbOfNodes;i++)
1164                 reverse_connectivity[nodes[i]].push_back(global_polyh_number);
1165               delete [] nodes;
1166             }
1167       }
1168
1169     // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
1170
1171   //calculate size of reverse_nodal_connectivity array
1172     int size_reverse_nodal_connectivity  = 0;
1173     for (int i = 1; i < _numberOfNodes+1; i++)
1174       size_reverse_nodal_connectivity += reverse_connectivity[i].size();
1175
1176     int * reverse_nodal_connectivity_index = new int[_numberOfNodes+1];
1177     int * reverse_nodal_connectivity = new int[size_reverse_nodal_connectivity];
1178
1179     reverse_nodal_connectivity_index[0] = 1;
1180     for (int i = 1; i < _numberOfNodes+1; i++)
1181       {
1182         int size = reverse_connectivity[i].size();
1183         reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
1184         for (int j = 0; j < size; j++)
1185           reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
1186       }
1187
1188     _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
1189                                                     reverse_nodal_connectivity_index,
1190                                                     reverse_nodal_connectivity,true);
1191   }
1192   END_OF(LOC);
1193 }
1194
1195 /*! If not yet done, calculate the Descending Connectivity */
1196 //-------------------------------------------------//
1197 void CONNECTIVITY::calculateDescendingConnectivity()
1198   //-------------------------------------------------//
1199   {
1200   const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
1201   BEGIN_OF(LOC);
1202
1203   if (_descending==NULL && _polygonsDescending==NULL && _polyhedronDescending==NULL)
1204     {
1205       if (_nodal==NULL)
1206         {
1207           MESSAGE(LOC<<"No connectivity found !");
1208           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
1209         }
1210       // calcul _descending from _nodal
1211       // we need CONNECTIVITY for constituent
1212
1213       if (_constituent != NULL)
1214         //      throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
1215         MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
1216
1217       if (_entityDimension == 3)
1218         _constituent = new CONNECTIVITY(MED_FACE);
1219       else if (_entityDimension == 2)
1220         _constituent = new CONNECTIVITY(MED_EDGE);
1221       else
1222         {
1223           MESSAGE(LOC<<"We are in 1D");
1224           return;
1225         }
1226
1227       _constituent->_typeConnectivity = MED_NODAL;
1228       _constituent->_numberOfNodes = _numberOfNodes;
1229       // foreach cells, we built array of constituent
1230       int DescendingSize = 0;
1231       for(int i=0; i<_numberOfTypes; i++)
1232         DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
1233       //_descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize);
1234       //const int * descend_connectivity = _descending->getValue();
1235       int * descend_connectivity = new int[DescendingSize];
1236       for (int i=0; i<DescendingSize; i++)
1237         descend_connectivity[i]=0;
1238       //const int * descend_connectivity_index = _descending->getIndex();
1239       int * descend_connectivity_index = new int[_count[_numberOfTypes]];
1240       descend_connectivity_index[0]=1;
1241       medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
1242       ConstituentsTypes[0]=MED_NONE;
1243       ConstituentsTypes[1]=MED_NONE;
1244       int * NumberOfConstituentsForeachType = new int[2];
1245       NumberOfConstituentsForeachType[0]=0;
1246       NumberOfConstituentsForeachType[1]=0;
1247       for(int i=0; i<_numberOfTypes; i++)
1248         {
1249           // initialize descend_connectivity_index array :
1250           int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
1251           for (int j=_count[i];j<_count[i+1];j++)
1252             {
1253               descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
1254               // compute number of constituent of all cell for each type
1255               for(int k=1;k<NumberOfConstituents+1;k++)
1256                 {
1257                   medGeometryElement MEDType = _type[i].getConstituentType(1,k);
1258                   if (ConstituentsTypes[0]==MED_NONE)
1259                     {
1260                       ConstituentsTypes[0]=MEDType;
1261                       NumberOfConstituentsForeachType[0]++;
1262                     } else if (ConstituentsTypes[0]==MEDType)
1263                       NumberOfConstituentsForeachType[0]++;
1264                   else if (ConstituentsTypes[1]==MED_NONE)
1265                     {
1266                       ConstituentsTypes[1]=MEDType;
1267                       NumberOfConstituentsForeachType[1]++;
1268                     } else if (ConstituentsTypes[1]==MEDType)
1269                       NumberOfConstituentsForeachType[1]++;
1270                   else
1271                     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
1272                 }
1273             }
1274         }
1275
1276       // we could built _constituent !
1277       int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
1278       int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
1279
1280       //_constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes);
1281
1282       // we use _constituent->_nodal
1283       //const int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
1284       //const int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
1285       int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
1286       int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
1287       ConstituentNodalConnectivityIndex[0]=1;
1288
1289       _constituent->_entityDimension = _entityDimension-1;
1290       if (ConstituentsTypes[1]==MED_NONE)
1291         _constituent->_numberOfTypes = 1;
1292       else
1293         _constituent->_numberOfTypes = 2;
1294       _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
1295       _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes];
1296 //CCRT      _constituent->_count = new med_int[_constituent->_numberOfTypes+1];
1297       _constituent->_count = new int[_constituent->_numberOfTypes+1];
1298       _constituent->_count[0]=1;
1299       med_int* tmp_NumberOfConstituentsForeachType = new med_int[_constituent->_numberOfTypes+1];
1300       tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
1301       for (int i=0; i<_constituent->_numberOfTypes;i++) {
1302         _constituent->_geometricTypes[i]=ConstituentsTypes[i];
1303         _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
1304         CELLMODEL Type(ConstituentsTypes[i]);
1305         _constituent->_type[i]=Type;
1306         tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
1307         for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
1308           ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
1309       }
1310       delete [] ConstituentsTypes;
1311       delete [] NumberOfConstituentsForeachType;
1312
1313       // we need reverse nodal connectivity
1314       if (! _reverseNodalConnectivity)
1315         calculateReverseNodalConnectivity();
1316       const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
1317       const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
1318
1319       // array to keep reverse descending connectivity
1320       int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
1321
1322       int TotalNumberOfSubCell = 0;
1323       for (int i=0; i<_numberOfTypes; i++)
1324         { // loop on all cell type
1325           CELLMODEL Type = _type[i];
1326           //      int NumberOfNodesPerCell = Type.getNumberOfNodes();
1327           int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
1328           for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
1329             for (int k=1; k<=NumberOfConstituentPerCell; k++)
1330               { // we loop on all sub cell of it
1331                 if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
1332                   // it is a new sub cell !
1333                   //        TotalNumberOfSubCell++;
1334                   // Which type ?
1335                   if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
1336                     tmp_NumberOfConstituentsForeachType[0]++;
1337                     TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
1338                   } else {
1339                     tmp_NumberOfConstituentsForeachType[1]++;
1340                     TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
1341                   }
1342                   //we have maximum two types
1343
1344                   descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
1345                   ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
1346                   int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
1347
1348                   int * NodesLists = new int[NumberOfNodesPerConstituent];
1349                   for (int l=0; l<NumberOfNodesPerConstituent; l++) {
1350                     NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
1351                     ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
1352                   }
1353                   // we use reverse_nodal_connectivity to find the other element which contain this sub cell
1354
1355                   // all elements which contains first node of sub cell :
1356                   int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
1357                   int ReverseNodalConnectivityIndex_1 = getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[0]-1);
1358                   //ReverseNodalConnectivityIndex[NodesLists[0]];
1359                   int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
1360
1361                   if (NumberOfCellsInList > 0)
1362                     { // we could have no element !
1363                       int * CellsList = new int[NumberOfCellsInList];
1364                       for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
1365                         CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
1366
1367                       // foreach node in sub cell, we search elements which are in common
1368                       // at the end, we must have only one !
1369
1370                       for (int l=1; l<NumberOfNodesPerConstituent; l++) {
1371                         int NewNumberOfCellsInList = 0;
1372                         int * NewCellsList = new int[NumberOfCellsInList];
1373                         for (int l1=0; l1<NumberOfCellsInList; l1++)
1374                           for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[l]-1); l2++)
1375                           //for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++)
1376                             {
1377                               if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
1378                                 // increasing order : CellsList[l1] are not in elements list of node l
1379                                 break;
1380                               if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
1381                                 // we have found one
1382                                 NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
1383                                 NewNumberOfCellsInList++;
1384                                 break;
1385                               }
1386                             }
1387                         NumberOfCellsInList = NewNumberOfCellsInList;
1388
1389                         delete [] CellsList;
1390                         CellsList = NewCellsList;
1391                       }
1392
1393                       if (NumberOfCellsInList > 0) { // We have found some elements !
1394                         int CellNumber = CellsList[0];
1395                         //delete [] CellsList;
1396                         if (NumberOfCellsInList>1)
1397                           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
1398
1399                         if (NumberOfCellsInList==1)
1400                           {
1401                             ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
1402
1403                             // we search sub cell number in this cell to not calculate it another time
1404                             // which type ?
1405                             CELLMODEL Type2;
1406                             for (int l=0; l<_numberOfTypes; l++)
1407                               if (CellNumber < _count[l+1]) {
1408                                 Type2=_type[l];
1409                                 break;
1410                               }
1411                             //int sub_cell_count2 = Type2.get_entities_count(1);
1412                             //int nodes_cell_count2 = Type2.get_nodes_count();
1413                             bool find2 = false;
1414                             for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) { // on all sub cell
1415                               int counter = 0;
1416                               for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
1417                                 for (int n=1; n<=Type.getConstituentType(1,k)%100; n++)
1418                                   {
1419                                     if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
1420                                       counter++;
1421                                   }
1422                               if (counter==Type.getConstituentType(1,k)%100)
1423                                 {
1424                                   descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
1425                                   find2 = true;
1426                                 }
1427                               if (find2)
1428                                 break;
1429                             }
1430                             if (!find2)
1431                               MESSAGE(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
1432                           }
1433                       } else {
1434                         ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
1435                       }
1436                       delete [] CellsList;
1437                     }
1438
1439                   delete [] NodesLists;
1440                 }
1441               }
1442         }
1443       // we adjust _constituent
1444       int NumberOfConstituent=0;
1445       int SizeOfConstituentNodal=0;
1446       for (int i=0;i<_constituent->_numberOfTypes; i++) {
1447         NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
1448         SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
1449       }
1450       // we built new _nodal attribute in _constituent
1451       //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
1452       //const int *ConstituentNodalValue = ConstituentNodal->getValue();
1453       //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
1454       int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
1455       int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
1456       ConstituentNodalIndex[0]=1;
1457       // we build _reverseDescendingConnectivity
1458       //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
1459       //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
1460       //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
1461       int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
1462       int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
1463       reverseDescendingConnectivityIndex[0]=1;
1464
1465       // first constituent type
1466       for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++) {
1467         ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1468         for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1469           ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1470         }
1471         reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1472         for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1473           reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
1474         }
1475       }
1476       // second type if any
1477       if (_constituent->_numberOfTypes==2) {
1478         int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
1479         int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1480         int offset2=offset*2;
1481         int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1482         for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++) {
1483           ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1484           for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1485             ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1486           }
1487           reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1488           for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1489             reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
1490           }
1491         }
1492         _constituent->_count[2]=NumberOfConstituent+1;
1493         // we correct _descending to adjust face number
1494         for(int j=0;j<DescendingSize;j++)
1495           if (abs(descend_connectivity[j])>tmp_NumberOfConstituentsForeachType[0]) {
1496             if ( descend_connectivity[j] > 0 )
1497               descend_connectivity[j]-=offset;
1498             else
1499               descend_connectivity[j]+=offset;
1500           }
1501       }
1502
1503       delete [] ConstituentNodalConnectivityIndex;
1504       delete [] ConstituentNodalConnectivity;
1505       delete [] ReverseDescendingConnectivityValue;
1506       _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1507       delete [] tmp_NumberOfConstituentsForeachType;
1508
1509       _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
1510                                         DescendingSize,
1511                                         descend_connectivity_index,
1512                                         descend_connectivity);
1513       delete [] descend_connectivity_index;
1514       delete [] descend_connectivity;
1515
1516       ////
1517       vector<int> PolyDescending;
1518       vector<int> Reversedescendingconnectivityvalue(reverseDescendingConnectivityValue,reverseDescendingConnectivityValue + 2*NumberOfConstituent);
1519       vector<int> Reversedescendingconnectivityindex(reverseDescendingConnectivityIndex,reverseDescendingConnectivityIndex + NumberOfConstituent);
1520       delete [] reverseDescendingConnectivityValue;
1521       delete [] reverseDescendingConnectivityIndex;
1522
1523
1524       // polygons (2D mesh)
1525
1526       vector<int> Constituentnodalvalue(ConstituentNodalValue,ConstituentNodalValue + SizeOfConstituentNodal);
1527       vector<int> Constituentnodalindex(ConstituentNodalIndex,ConstituentNodalIndex + NumberOfConstituent+1);
1528       delete [] ConstituentNodalValue;
1529       delete [] ConstituentNodalIndex;
1530       int NumberOfNewSeg = 0;
1531
1532       for (int i=0; i <getNumberOfPolygons(); i++) // for each polygon
1533         {
1534           const int * vector_begin = &_polygonsNodal->getValue()[_polygonsNodal->getIndex()[i]-1];
1535           int vector_size = _polygonsNodal->getIndex()[i+1]-_polygonsNodal->getIndex()[i]+1;
1536           vector<int> myPolygon(vector_begin, vector_begin+vector_size);
1537           myPolygon[myPolygon.size()-1] = myPolygon[0]; // because first and last point make a segment
1538
1539           for (int j=0; j<myPolygon.size()-1; j++) // for each segment of polygon
1540             {
1541               MEDMODULUSARRAY segment_poly(2,&myPolygon[j]);
1542               int ret_compare = 0;
1543
1544               // we search it in existing segments
1545
1546               for (int k=0; k<Constituentnodalindex.size()-1; k++)
1547                 {
1548                   MEDMODULUSARRAY segment(2,&Constituentnodalvalue[0] + Constituentnodalindex[k]-1);
1549                   ret_compare = segment_poly.compare(segment);
1550                   if (ret_compare) // segment_poly already exists
1551                     {
1552                       PolyDescending.push_back(ret_compare*(k+1)); // we had it to the connectivity
1553                       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)
1554                       break;
1555                     }
1556                 }
1557
1558               // segment_poly must be created
1559
1560               if (!ret_compare)
1561                 {
1562                   NumberOfNewSeg++;
1563                   PolyDescending.push_back(NumberOfConstituent+NumberOfNewSeg); // we had it to the connectivity
1564                   Constituentnodalvalue.push_back(segment_poly[0]);
1565                   Constituentnodalvalue.push_back(segment_poly[1]);
1566                   Constituentnodalindex.push_back(Constituentnodalindex[Constituentnodalindex.size()-1] + 2); // we have only segments
1567                   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)
1568                   insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewSeg-1, 2*(NumberOfConstituent+NumberOfNewSeg-1)+1); // idem with index
1569                 }
1570             }
1571         }
1572
1573       if (getNumberOfPolygons() > 0)
1574         {
1575           _polygonsDescending = new MEDSKYLINEARRAY(getNumberOfPolygons(),_polygonsNodal->getLength(),_polygonsNodal->getIndex(),&PolyDescending[0]); // index are the same for polygons nodal and descending connectivities
1576
1577           NumberOfConstituent += NumberOfNewSeg;
1578           SizeOfConstituentNodal += 2*NumberOfNewSeg;
1579           _constituent->_count[1] = NumberOfConstituent+1;
1580         }
1581
1582
1583       // polyhedron (3D mesh)
1584
1585       vector<int> Constituentpolygonsnodalvalue;
1586       vector<int> Constituentpolygonsnodalindex(1,1);
1587       int NumberOfNewFaces = 0; // by convention new faces are polygons
1588
1589       for (int i=0; i<getNumberOfPolyhedron(); i++) // for each polyhedron
1590         {
1591           // we create a POLYHEDRONARRAY containing only this polyhedra
1592           int myNumberOfFaces = _polyhedronNodal->getPolyhedronIndex()[i+1]-_polyhedronNodal->getPolyhedronIndex()[i];
1593           int myNumberOfNodes = _polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i+1]-1]-_polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i]-1];
1594           POLYHEDRONARRAY myPolyhedra(1,myNumberOfFaces,myNumberOfNodes);
1595 //CCRT    vector<int> myFacesIndex(_polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1, _polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1 + myNumberOfFaces+1);
1596           vector<med_int> myFacesIndex(_polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1, _polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1 + myNumberOfFaces+1);
1597           for (int j=myNumberOfFaces; j>=0; j--)
1598             myFacesIndex[j] -= myFacesIndex[0]-1;
1599           myPolyhedra.setFacesIndex(&myFacesIndex[0]);
1600           myPolyhedra.setNodes(_polyhedronNodal->getNodes() + _polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i]-1]-1);
1601
1602           for (int j=0; j<myPolyhedra.getNumberOfFaces(); j++) // for each face of polyhedra
1603             {
1604               int myFaceNumberOfNodes = myPolyhedra.getFacesIndex()[j+1]-myPolyhedra.getFacesIndex()[j];
1605 //CCRT
1606 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000)
1607               int ii ;
1608               const med_int * Nodes = myPolyhedra.getNodes() ;
1609               int * tmp_Nodes = new int[myPolyhedra.getNumberOfNodes()] ;
1610               for ( ii = 0 ; ii < myPolyhedra.getNumberOfNodes() ; ii++ )
1611                  tmp_Nodes[ii] = Nodes[ii] ;
1612               const med_int * FacesIndex = myPolyhedra.getFacesIndex() ;
1613               int * tmp_FacesIndex = new int[myPolyhedra.getNumberOfFaces()+1] ;
1614               for ( ii = 0 ; ii < myPolyhedra.getNumberOfNodes() ; ii++ )
1615                  tmp_FacesIndex[ii] = FacesIndex[ii] ;
1616 //CCRT : copy of Nodes
1617               MEDMODULUSARRAY face_poly(myFaceNumberOfNodes,tmp_Nodes + tmp_FacesIndex[j]-1);
1618 //CCRT              delete [] tmp_Nodes ;
1619               delete [] tmp_FacesIndex ;
1620 #else
1621               MEDMODULUSARRAY face_poly(myFaceNumberOfNodes,myPolyhedra.getNodes() + myPolyhedra.getFacesIndex()[j]-1);
1622 #endif
1623               int ret_compare = 0;
1624
1625               // we search it in existing faces
1626
1627               // we search first in TRIA3 and QUAD4
1628               if (myFaceNumberOfNodes == 3 || myFaceNumberOfNodes == 4)
1629                 {
1630                   int Begin = -1; // first TRIA3 or QUAD4
1631                   int Number = 0; // number of TRIA3 or QUAD4
1632                   for (int k=0; k<Constituentnodalindex.size()-1; k++)
1633                     if (Constituentnodalindex[k+1]-Constituentnodalindex[k] == myFaceNumberOfNodes)
1634                       {
1635                         if (Begin == -1)
1636                           Begin = k;
1637                         Number++;
1638                       }
1639
1640                   for (int k=0; k<Number; k++)
1641                     {
1642                       MEDMODULUSARRAY face(myFaceNumberOfNodes,&Constituentnodalvalue[0] + Constituentnodalindex[Begin+k]-1);
1643                       ret_compare = face_poly.compare(face);
1644                       if (ret_compare)
1645                         {
1646                           PolyDescending.push_back(ret_compare*(Begin+k+1)); // we had it to the connectivity
1647                           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)
1648                           break;
1649                         }
1650                     }
1651                 }
1652
1653               // we search last in POLYGONS
1654               if (!ret_compare)
1655                 {
1656                   for (int k=0; k<static_cast<int>(Constituentpolygonsnodalindex.size())-1; k++) // we must cast the unsigned int into int before doing -1
1657                     {
1658                       if (Constituentpolygonsnodalindex[k+1]-Constituentpolygonsnodalindex[k] == myFaceNumberOfNodes)
1659                         {
1660                           MEDMODULUSARRAY face(myFaceNumberOfNodes,&Constituentpolygonsnodalvalue[0] + Constituentpolygonsnodalindex[k]-1);
1661                           ret_compare = face_poly.compare(face);
1662                           if (ret_compare)
1663                             {
1664                               PolyDescending.push_back(ret_compare*(NumberOfConstituent+k+1)); // we had it to the connectivity
1665                               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)
1666                               break;
1667                             }
1668                         }
1669                     }
1670                 }
1671
1672               // if not found, face_poly must be created
1673
1674               if (!ret_compare)
1675                 {
1676                   NumberOfNewFaces++;
1677                   PolyDescending.push_back(NumberOfConstituent+NumberOfNewFaces); // we had it to the connectivity
1678                   for (int k=0; k<myFaceNumberOfNodes; k++)
1679                     Constituentpolygonsnodalvalue.push_back(face_poly[k]);
1680                   Constituentpolygonsnodalindex.push_back(Constituentpolygonsnodalindex[Constituentpolygonsnodalindex.size()-1]+myFaceNumberOfNodes);
1681                   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)
1682                   insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewFaces-1, 2*(NumberOfConstituent+NumberOfNewFaces-1)+1); // idem with index
1683                 }
1684             }
1685         }
1686
1687       if (getNumberOfPolyhedron() > 0)
1688         {
1689 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000)
1690           int * tmp_PolyhedronIndex = new int[getNumberOfPolyhedron()+1] ;
1691           const MED_EN::med_int * PolyhedronIndex = _polyhedronNodal->getPolyhedronIndex() ;
1692           int ii ;
1693           for ( ii = 0 ; ii < getNumberOfPolyhedron()+1 ; ii++ )
1694              tmp_PolyhedronIndex[ii] = PolyhedronIndex[ii] ;
1695 //CCRT : copy of PolyhedronIndex
1696           _polyhedronDescending = new MEDSKYLINEARRAY(getNumberOfPolyhedron(),_polyhedronNodal->getNumberOfFaces(),tmp_PolyhedronIndex,&PolyDescending[0]); // polyhedron index are the same for nodal and descending connectivities
1697 #else
1698           _polyhedronDescending = new MEDSKYLINEARRAY(getNumberOfPolyhedron(),_polyhedronNodal->getNumberOfFaces(),_polyhedronNodal->getPolyhedronIndex(),&PolyDescending[0]); // polyhedron index are the same for nodal and descending connectivities
1699 #endif
1700
1701           if (_constituent->_polygonsNodal != NULL)
1702             delete [] _constituent->_polygonsNodal;
1703           _constituent->_polygonsNodal = new MEDSKYLINEARRAY(Constituentpolygonsnodalindex.size()-1,Constituentpolygonsnodalvalue.size(),&Constituentpolygonsnodalindex[0],&Constituentpolygonsnodalvalue[0]);
1704         }
1705
1706       // delete _constituent->_nodal
1707       _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
1708                                                  SizeOfConstituentNodal,
1709                                                  &Constituentnodalindex[0],
1710                                                  &Constituentnodalvalue[0]);
1711
1712       int NumberOfConstituentWithPolygons = NumberOfConstituent + NumberOfNewFaces;
1713       Reversedescendingconnectivityindex.push_back(Reversedescendingconnectivityindex[Reversedescendingconnectivityindex.size()-1]+2); // we complete the index
1714       _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituentWithPolygons+1,
1715                                                            2*NumberOfConstituentWithPolygons,
1716                                                            &Reversedescendingconnectivityindex[0],
1717                                                            &Reversedescendingconnectivityvalue[0]);
1718       ////
1719     }
1720   END_OF(LOC);
1721   }
1722
1723 /*! Not implemented yet */
1724 //--------------------------------------------------------------------//
1725 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
1726 //--------------------------------------------------------------------//
1727 {
1728   const char * LOC = "CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity) : ";
1729   BEGIN_OF(LOC);
1730
1731   MESSAGE(LOC<<"method not yet implemented " << myConnectivity._entity);
1732   // Mesh dimension !
1733
1734   END_OF(LOC);
1735   return;
1736 }
1737
1738
1739 /*!
1740   Returns the geometry of an element given by its entity type & its global number.
1741
1742   Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
1743 */
1744 //--------------------------------------------------------------------//
1745 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalNumber) const
1746 //--------------------------------------------------------------------//
1747 {
1748   const char * LOC = "medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity, int globalNumber) const : ";
1749   BEGIN_OF(LOC);
1750   int globalNumberMin = 1;
1751   int globalNumberMax ;
1752
1753   if (_entity==Entity) globalNumberMax = _count[_numberOfTypes];
1754   else if (_constituent!=NULL) globalNumberMax = _constituent->_count[_constituent->_numberOfTypes];
1755   else
1756     throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1757
1758   // The globalNumber must verify  : 1 <=  globalNumber < _count[_numberOfTypes] (== totalNumberOfElement+1)
1759
1760   if ( (globalNumber < globalNumberMin) || (globalNumber >  globalNumberMax-1 )  )
1761     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |"
1762                                  << globalNumberMin <<"| and <= |" << globalNumberMax-1 ));
1763
1764   if (_entity==Entity) {
1765     for(int i=1; i<=_numberOfTypes;i++)
1766       if (globalNumber<_count[i])
1767         return _geometricTypes[i-1];
1768   }
1769   else if (_constituent!=NULL)
1770     return _constituent->getElementType(Entity,globalNumber);
1771   else
1772     throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1773   throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
1774
1775   END_OF(LOC);
1776 }
1777
1778 /*!
1779   Idem as getElementType method except that it manages polygon and polyhedron.
1780 */
1781 //--------------------------------------------------------------------//
1782 medGeometryElement CONNECTIVITY::getElementTypeWithPoly(medEntityMesh Entity,int globalNumber) const
1783 //--------------------------------------------------------------------//
1784 {
1785   int globalNumberMaxOfClassicType;
1786   if(_entity==Entity)
1787     {
1788       globalNumberMaxOfClassicType=_count[_numberOfTypes];
1789       if(globalNumber>=1)
1790         {
1791           if(globalNumber<globalNumberMaxOfClassicType)
1792             {
1793               for(int i=1; i<=_numberOfTypes;i++)
1794                 if (globalNumber<_count[i])
1795                   return _geometricTypes[i-1];
1796               throw MEDEXCEPTION("never happens just for compilo");
1797             }
1798           else
1799             {
1800               int localNumberInPolyArray=globalNumber-globalNumberMaxOfClassicType;
1801               int nbOfPol=getNumberOfElementOfPolyType(_entity);
1802               if(localNumberInPolyArray<nbOfPol)
1803                 return getPolyTypeRelativeTo();
1804               else
1805                 throw MEDEXCEPTION("getElementTypeWithPoly : unexisting element type");
1806             }
1807         }
1808       else
1809         throw MEDEXCEPTION("getElementTypeWithPoly : globalNumber < 1");
1810     }
1811   else
1812     {
1813       if(_constituent!=NULL)
1814         return _constituent->getElementTypeWithPoly(Entity,globalNumber);
1815       else
1816         throw MEDEXCEPTION("getElementTypeWithPoly : Entity not defined !");
1817     }
1818 }
1819
1820 ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
1821 {
1822     os << endl << "------------- Entity = ";
1823     switch (co._entity)
1824     {
1825         case MED_CELL:
1826             os << "MED_CELL";
1827             break;
1828         case MED_FACE:
1829             os << "MED_FACE";
1830             break;
1831         case MED_EDGE:
1832             os << "MED_EDGE";
1833             break;
1834         case MED_NODE:
1835             os << "MED_NODE";
1836             break;
1837         case MED_ALL_ENTITIES:
1838             os << "MED_ALL_ENTITIES";
1839             break;
1840         default:
1841             os << "Unknown";
1842             break;
1843     }
1844     os  << " -------------\n\nMedConnectivity : ";
1845     switch (co._typeConnectivity)
1846     {
1847         case MED_NODAL:
1848             os << "MED_NODAL\n";
1849             break;
1850         case MED_DESCENDING:
1851             os << "MED_DESCENDING\n";
1852             break;
1853         default:
1854             break;
1855     }
1856     os << "Entity dimension : " << co._entityDimension << endl;
1857     os << "Number of nodes : " << co._numberOfNodes << endl;
1858     os << "Number of types : " << co._numberOfTypes << endl;
1859     for (int i=0; i!=co._numberOfTypes ; ++i)
1860         os << "  -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : "
1861            << co._count[i+1]-co._count[i] << " elements" << endl;
1862
1863     if (co._typeConnectivity == MED_NODAL )
1864     {
1865         for (int i=0; i<co._numberOfTypes; i++)
1866         {
1867             os << endl << co._type[i].getName() << " : " << endl;
1868             int numberofelements = co._count[i+1]-co._count[i];
1869             const int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, co._geometricTypes[i]);
1870             int numberofnodespercell = co._geometricTypes[i]%100;
1871             for (int j=0;j<numberofelements;j++)
1872             {
1873                 os << setw(6) << j+1 << " : ";
1874                 for (int k=0;k<numberofnodespercell;k++)
1875                     os << connectivity[j*numberofnodespercell+k]<<" ";
1876                 os << endl;
1877             }
1878         }
1879
1880         os << endl << "MED_POLYGON : " << co.getNumberOfPolygons() << " polygons" << endl;
1881         if (co.getNumberOfPolygons() > 0)
1882           {
1883             const int* polygonsconnectivity = co.getPolygonsConnectivity(MED_NODAL,co.getEntity());
1884             const int* polygonsconnectivityindex = co.getPolygonsConnectivityIndex(MED_NODAL,co.getEntity());
1885
1886             for (int i=0; i<co.getNumberOfPolygons(); i++)
1887               {
1888                 int numberofnodesperpolygon = polygonsconnectivityindex[i+1]-polygonsconnectivityindex[i];
1889                 for (int j=0; j<numberofnodesperpolygon; j++)
1890                   os << polygonsconnectivity[polygonsconnectivityindex[i]-1+j] << " ";
1891                 os << endl;
1892               }
1893           }
1894
1895         os << endl << "MED_POLYHEDRA : " << co.getNumberOfPolyhedron() << " polyhedron" << endl;
1896         if (co.getNumberOfPolyhedron() > 0)
1897           {
1898             const int* polyhedronconnectivity = co.getPolyhedronConnectivity(MED_NODAL);
1899             const int* polyhedronfacesindex = co.getPolyhedronFacesIndex();
1900             const int* polyhedronindex = co.getPolyhedronIndex(MED_NODAL);
1901
1902             for (int i=0; i<co.getNumberOfPolyhedron(); i++)
1903               {
1904                 int numberoffacesperpolyhedra = polyhedronindex[i+1]-polyhedronindex[i];
1905                 for (int j=0; j<numberoffacesperpolyhedra; j++)
1906                   {
1907                     int numberofnodesperface = polyhedronfacesindex[polyhedronindex[i]-1+j+1]-polyhedronfacesindex[polyhedronindex[i]-1+j];
1908                     for (int k=0; k<numberofnodesperface; k++)
1909                       os << polyhedronconnectivity[polyhedronfacesindex[polyhedronindex[i]-1+j]-1+k] << " ";
1910                     if (j != numberoffacesperpolyhedra-1) os << "| ";
1911                   }
1912                 os << endl;
1913               }
1914           }
1915
1916     }
1917     else if (co._typeConnectivity == MED_DESCENDING)
1918     {
1919         int numberofelements = co.getNumberOf( co._entity , MED_ALL_ELEMENTS);
1920         if (numberofelements > 0)
1921           {
1922             const int *connectivity =  co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
1923             const int *connectivity_index =  co.getConnectivityIndex( co._typeConnectivity, co._entity );
1924
1925             for ( int j=0; j!=numberofelements; j++ )
1926               {
1927                 os << "Element " << j+1 << " : ";
1928                 for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
1929                   os << connectivity[k-1] << " ";
1930                 os << endl;
1931               }
1932           }
1933
1934         os << endl << "MED_POLYGON : " << co.getNumberOfPolygons() << " polygons" << endl;
1935         if (co.getNumberOfPolygons() > 0)
1936           {
1937             const int* polygonsconnectivity = co.getPolygonsConnectivity(MED_DESCENDING,co.getEntity());
1938             const int* polygonsconnectivityindex = co.getPolygonsConnectivityIndex(MED_DESCENDING,co.getEntity());
1939
1940             for (int i=0; i<co.getNumberOfPolygons(); i++)
1941               {
1942                 int numberofedgesperpolygon = polygonsconnectivityindex[i+1]-polygonsconnectivityindex[i];
1943                 for (int j=0; j<numberofedgesperpolygon; j++)
1944                   os << polygonsconnectivity[polygonsconnectivityindex[i]-1+j] << " ";
1945                 os << endl;
1946               }
1947           }
1948
1949         os << endl << "MED_POLYHEDRA : " << co.getNumberOfPolyhedron() << " polyhedron" << endl;
1950         if (co.getNumberOfPolyhedron() > 0)
1951           {
1952             const int* polyhedronconnectivity = co.getPolyhedronConnectivity(MED_DESCENDING);
1953             const int* polyhedronindex = co.getPolyhedronIndex(MED_DESCENDING);
1954
1955             for (int i=0; i<co.getNumberOfPolyhedron(); i++)
1956               {
1957                 int numberoffacesperpolyhedra = polyhedronindex[i+1]-polyhedronindex[i];
1958                 for (int j=0; j<numberoffacesperpolyhedra; j++)
1959                   os << polyhedronconnectivity[polyhedronindex[i]-1+j] << " ";
1960                 os << endl;
1961               }
1962           }
1963
1964     }
1965
1966     if (co._constituent)
1967         os << endl << *co._constituent << endl;
1968
1969     return os;
1970 }
1971
1972 /*!
1973   method that adds to vector 'nodes' all the nodes of polyhedron with id 'polyhedronId'.
1974   WARNING the returned pointer should be deallocated. Returned nodes and polyhedronId are in form [1,...]
1975  */
1976 int *CONNECTIVITY::getNodesOfPolyhedron(int polyhedronId, int& lgthOfTab) const
1977 {
1978   const int *nodes=getPolyhedronConnectivity(MED_EN::MED_NODAL);
1979   const int *faces=getPolyhedronFacesIndex();
1980   const int *glob=getPolyhedronIndex(MED_EN::MED_NODAL);
1981   int offsetWithClassicType=getNumberOf(_entity,MED_ALL_ELEMENTS);
1982   set<int> retInSet;
1983   int startFace=glob[polyhedronId-offsetWithClassicType-1]-1;
1984   int endFace=glob[polyhedronId-offsetWithClassicType]-1;
1985   int i;
1986   for(i=startFace;i!=endFace;i++)
1987     {
1988       for(int j=faces[i];j<faces[i+1];j++)
1989         retInSet.insert(nodes[j-1]);
1990     }
1991   lgthOfTab=retInSet.size();
1992   int *ret=new int[lgthOfTab];
1993   set<int>::iterator iter=retInSet.begin();
1994   i=0;
1995   for(iter=retInSet.begin();iter!=retInSet.end();iter++)
1996     ret[i++]=*iter;
1997   return ret;
1998 }
1999
2000 /*!
2001   Idem as MESH::getNodesOfPolyhedron except that returned nodes are sorted by face. 'nbOfNodesPerFaces' is an array of size 'nbOfFaces'.
2002   Returned int** has a size of 'nbOfNodesPerFaces' too, and for each element j in int** the size is  nbOfNodesPerFaces[j].
2003   Warning both returned 'nbOfNodesPerFaces' and returned value should be deallocated. Returned nodes and 'polyhedronId' are in form [1,...]
2004  */
2005 int **CONNECTIVITY::getNodesPerFaceOfPolyhedron(int polyhedronId, int& nbOfFaces, int* & nbOfNodesPerFaces) const
2006 {
2007   int i;
2008   const int *nodes=getPolyhedronConnectivity(MED_EN::MED_NODAL);
2009   const int *faces=getPolyhedronFacesIndex();
2010   const int *glob=getPolyhedronIndex(MED_EN::MED_NODAL);
2011   int offsetWithClassicType=getNumberOf(_entity,MED_ALL_ELEMENTS);
2012
2013   int startFace=glob[polyhedronId-offsetWithClassicType-1]-1;
2014   nbOfFaces=glob[polyhedronId-offsetWithClassicType]-startFace-1;
2015   nbOfNodesPerFaces=new int[nbOfFaces];
2016   int **ret=new int* [nbOfFaces];
2017   for(i=0;i<nbOfFaces;i++)
2018     {
2019       int startNode=faces[startFace+i]-1;
2020       nbOfNodesPerFaces[i]=faces[startFace+i+1]-startNode-1;
2021       ret[i]=(int *)(nodes)+startNode;
2022     }
2023   return ret;
2024 }
2025
2026 int CONNECTIVITY::getNumberOfTypesWithPoly(MED_EN::medEntityMesh Entity) const
2027 {
2028   SCRUTE(_entity);
2029   SCRUTE(Entity);
2030   if (_entity==Entity)
2031     {
2032       SCRUTE(_numberOfTypes);
2033       SCRUTE(getNumberOfPolyType());
2034       return _numberOfTypes+getNumberOfPolyType();
2035     }
2036   else if (_constituent!=NULL)
2037     return _constituent->getNumberOfTypesWithPoly(Entity);
2038   else
2039     return 0;
2040 }
2041
2042 int CONNECTIVITY::getNumberOfPolyType()  const
2043 {
2044   if (_entity==MED_CELL && _entityDimension==3)
2045     {
2046       if(getNumberOfPolyhedron()>0)
2047         return 1;
2048     }
2049   else if ((_entity==MED_CELL && _entityDimension==2) || (_entity==MED_FACE && _entityDimension==2))
2050     if(getNumberOfPolygons()>0)
2051       return 1;
2052   return 0;
2053 }
2054
2055 int CONNECTIVITY::getNumberOfElementOfPolyType(MED_EN::medEntityMesh Entity)  const
2056 {
2057   if(_entity==Entity)
2058     {
2059       if (_entity==MED_CELL && _entityDimension==3)
2060         return getNumberOfPolyhedron();
2061       else if ((_entity==MED_CELL && _entityDimension==2) || (_entity==MED_FACE && _entityDimension==2))
2062         return getNumberOfPolygons();
2063       return 0;
2064     }
2065   else
2066     {
2067       if(_constituent!=NULL)
2068         return _constituent->getNumberOfElementOfPolyType(Entity);
2069       else
2070         throw MEDEXCEPTION("_constituent required to evaluate getNumberOfElementOfPolyType");
2071     }
2072 }
2073
2074 /*
2075   Method equivalent to getGeometricTypes except that it includes not only classical Types but polygons/polyhedra also.
2076   WARNING the returned array MUST be deallocated.
2077  */
2078 MED_EN::medGeometryElement * CONNECTIVITY::getGeometricTypesWithPoly(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
2079 {
2080   if(_entity==Entity)
2081     {
2082       int nbOfTypesTotal=getNumberOfTypesWithPoly(Entity);
2083       int nbOfTypesWithoutPoly=getNumberOfTypes(Entity);
2084       medGeometryElement *ret=new medGeometryElement[nbOfTypesTotal];
2085       memcpy(ret,getGeometricTypes(Entity),nbOfTypesWithoutPoly*sizeof(medGeometryElement));
2086       if(nbOfTypesTotal>nbOfTypesWithoutPoly)
2087         {
2088           if (Entity==MED_CELL && _entityDimension==3)
2089             ret[nbOfTypesWithoutPoly]=MED_EN::MED_POLYHEDRA;
2090           else if((Entity==MED_CELL && _entityDimension==2) || (Entity==MED_FACE && _entityDimension==2))
2091             ret[nbOfTypesWithoutPoly]=MED_EN::MED_POLYGON;
2092         }
2093       return ret;
2094     }
2095   else
2096     {
2097       if(_constituent)
2098         return _constituent->getGeometricTypesWithPoly(Entity);
2099       throw MEDEXCEPTION("_constituent required to evaluate getGeometricTypesWithPoly");
2100     }
2101 }
2102
2103 /*
2104   Method used in CalculateDescendingConnectivity. So it's typically a private method.
2105   The aim of this method is to hide to CalculateDescendingConnectivity algorithm the fact that in reverse connectivity polygons and polyhedrons
2106   are not in separate data structure, contrary to not reverse connectivities.
2107  */
2108 int CONNECTIVITY::getIndexOfEndClassicElementInReverseNodal(const int *reverseNodalValue, const int *reverseNodalIndex, int rk)  const
2109 {
2110   int nbOfLastElt=getNumberOf(_entity,MED_ALL_ELEMENTS);
2111   int ret=reverseNodalIndex[rk];
2112   for(int i=reverseNodalIndex[rk];i<reverseNodalIndex[rk+1];i++)
2113     {
2114       if(reverseNodalValue[i-1]<=nbOfLastElt)
2115         ret++;
2116     }
2117   return ret;
2118 }
2119
2120 /*
2121   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.
2122   This method has to be applied on a instance of MEDMEM::CONNECTIVITY with _entityDimension==3.
2123   WARNING calculateDescendingConnectivity must have been called before.
2124  */
2125 void CONNECTIVITY::invertConnectivityForAFace(int faceId, const int *nodalConnForFace, bool polygonFace)
2126 {
2127   // we have always 2 neighbourings
2128   int cell1 = _reverseDescendingConnectivity->getIJ(faceId,1);
2129   int cell2 = _reverseDescendingConnectivity->getIJ(faceId,2);
2130
2131   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
2132     // Updating _reverseDescendingConnectivity
2133     _reverseDescendingConnectivity->setIJ(faceId,1,cell2);
2134     _reverseDescendingConnectivity->setIJ(faceId,2,cell1);
2135     // Updating _constituent->_nodal because of reversity
2136     const int *descendingNodalIndex=(!polygonFace)?_constituent->_nodal->getIndex():_constituent->_polygonsNodal->getIndex();
2137     MEDSKYLINEARRAY *currentNodal=(!polygonFace)?_constituent->_nodal:_constituent->_polygonsNodal;
2138     int faceIdRelative=(!polygonFace)?faceId:faceId-getNumberOf(MED_FACE,MED_ALL_ELEMENTS);
2139     for(int iarray=1;iarray<=(descendingNodalIndex[faceIdRelative]-descendingNodalIndex[faceIdRelative-1]);iarray++)
2140       currentNodal->setIJ(faceIdRelative,iarray,nodalConnForFace[iarray-1]);
2141
2142     // Updating _descending for cell1 and cell2
2143     const int NB_OF_CELLS_SHARING_A_FACE=2;
2144     int cellsToUpdate[NB_OF_CELLS_SHARING_A_FACE]; cellsToUpdate[0]=cell1; cellsToUpdate[1]=cell2;
2145     for(int curCell=0;curCell<NB_OF_CELLS_SHARING_A_FACE;curCell++)
2146       {
2147         int cell=cellsToUpdate[curCell];
2148         bool polyhCell=(getElementTypeWithPoly(MED_CELL,cell)==MED_POLYHEDRA);
2149         if(polyhCell)
2150           cell-=getNumberOf(MED_CELL,MED_ALL_ELEMENTS);
2151         const int *newDescendingIndex=(!polyhCell)?_descending->getIndex():_polyhedronDescending->getIndex();
2152         MEDSKYLINEARRAY *currentDescending=(!polyhCell)?_descending:_polyhedronDescending;
2153         for(int iface=newDescendingIndex[cell-1];iface<newDescendingIndex[cell];iface++)
2154           {
2155             int curValue=currentDescending->getIndexValue(iface);
2156             if (abs(curValue)==faceId)
2157               currentDescending->setIndexValue(iface,-curValue);
2158           }
2159       }
2160   }
2161 }
2162
2163 /*
2164   Method with 2 output : the connectivity required and its length 'lgth'. This method gives the connectivity independently it is a polygons/polyhedrons or normal
2165   element.
2166  */
2167 const int * CONNECTIVITY::getConnectivityOfAnElementWithPoly(MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity, int Number, int &lgth)
2168 {
2169   if(Entity==MED_EN::MED_NODE)
2170     throw  MEDEXCEPTION("No connectivity attached to a node entity");
2171   if(Entity==_entity)
2172     {
2173       if(_entity==MED_EDGE && _entityDimension==1)
2174         {
2175           const int * newConstituentValue = 0;
2176           const int * newConstituentIndex = 0;
2177           newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
2178           newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
2179           lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
2180           return newConstituentValue+newConstituentIndex[Number-1]-1;
2181         }
2182       int nbOfClassicalElements=getNumberOf(_entity,MED_ALL_ELEMENTS);
2183       if((_entity==MED_FACE && _entityDimension==2) || (_entity==MED_CELL && _entityDimension==2))
2184         {
2185           const int * newConstituentValue = 0;
2186           const int * newConstituentIndex = 0;
2187           if(Number<=nbOfClassicalElements)
2188             {
2189               newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
2190               newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
2191               lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
2192               return newConstituentValue+newConstituentIndex[Number-1]-1;
2193             }
2194           else
2195             {
2196               int localNumber=Number-nbOfClassicalElements-1;
2197               if(localNumber<getNumberOfPolygons())
2198                 {
2199                   newConstituentValue = getPolygonsConnectivity(ConnectivityType,Entity);
2200                   newConstituentIndex = getPolygonsConnectivityIndex(ConnectivityType,Entity);
2201                   lgth=newConstituentIndex[localNumber+1]-newConstituentIndex[localNumber];
2202                   return newConstituentValue+newConstituentIndex[localNumber]-1;
2203                 }
2204               else
2205                 throw  MEDEXCEPTION("Unknown number");
2206             }
2207         }
2208       else if(_entity==MED_CELL && _entityDimension==3)
2209         {
2210           const int * newConstituentValue = 0;
2211           const int * newConstituentIndex = 0;
2212           if(Number<=nbOfClassicalElements)
2213             {
2214               newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
2215               newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
2216               lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
2217               return newConstituentValue+newConstituentIndex[Number-1]-1;
2218             }
2219           else
2220             {
2221               int localNumber=Number-nbOfClassicalElements-1;
2222               if(localNumber<getNumberOfPolyhedron())
2223                 {
2224                   if(ConnectivityType==MED_NODAL)
2225                     throw  MEDEXCEPTION("NODAL Connectivity required for a polyhedron");
2226 //                newConstituentValue = _polyhedronDescending->getValue();
2227 //                newConstituentIndex = _polyhedronDescending->getIndex();
2228                   newConstituentValue = getPolyhedronConnectivity( ConnectivityType );
2229                   newConstituentIndex = getPolyhedronIndex( ConnectivityType );
2230                   lgth=newConstituentIndex[localNumber+1]-newConstituentIndex[localNumber];
2231                   return newConstituentValue+newConstituentIndex[localNumber]-1;
2232                 }
2233               else
2234                 throw  MEDEXCEPTION("Unknown number");
2235             }
2236         }
2237       else
2238         throw MEDEXCEPTION("No connectivity available");
2239     }
2240   else
2241     {
2242       if(_constituent==NULL)
2243         calculateDescendingConnectivity();
2244       return _constituent->getConnectivityOfAnElementWithPoly(ConnectivityType,Entity,Number,lgth);
2245     }
2246 }
2247
2248 int CONNECTIVITY::getNumberOfElementsWithPoly(MED_EN::medEntityMesh Entity, MED_EN::medGeometryElement Type) const
2249 {
2250   if(Entity==MED_EN::MED_NODE)
2251     return _numberOfNodes;
2252   if(Entity==_entity)
2253     {
2254       if(Type==MED_POLYGON || Type==MED_POLYHEDRA)
2255         return getNumberOfElementOfPolyType(_entity);
2256       else
2257         return getNumberOf(_entity,Type);
2258     }
2259   else
2260     {
2261       if(_constituent)
2262         return _constituent->getNumberOfElementsWithPoly(Entity,Type);
2263       else
2264         throw MEDEXCEPTION("CONNECTIVITY::getNumberOfElementsWithPoly : _constituent needed");
2265     }
2266 }
2267
2268 /*!
2269   Perform a deep comparison of the 2 connectivities in NODAL mode and on all elements.
2270 */
2271 bool CONNECTIVITY::deepCompare(const CONNECTIVITY& other) const
2272 {
2273   CONNECTIVITY* temp=(CONNECTIVITY* )this;
2274   const int *conn1=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2275   int size1=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2276   temp=(CONNECTIVITY* )(&other);
2277   const int *conn2=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2278   int size2=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2279   if(size1!=size2)
2280     return false;
2281   bool ret=true;
2282   for(int i=0;i<size1 && ret;i++)
2283     {
2284       ret=(conn1[i]==conn2[i]);
2285     }
2286   return ret;
2287 }