src/MEDMEM/MEDMEM_Connectivity.cxx

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