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