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