src/MEDMEM/MEDMEM_Connectivity.cxx

   1 using namespace std;
   2 #include "MEDMEM_Connectivity.hxx"
   3 #include "MEDMEM_Family.hxx"
   4 #include "MEDMEM_Group.hxx"
   5 #include "MEDMEM_CellModel.hxx"
   6
   7 #include "MEDMEM_SkyLineArray.hxx"
   8 #include "MEDMEM_ModulusArray.hxx"
   9
  10 #include "MEDMEM_STRING.hxx"
  11 #include <iomanip>
  12
  13 using namespace MEDMEM;
  14 /*!
  15    Default Constructor. /n
  16    Default for Entity is MED_CELL and type of Connectivity is MED_NODAL */
  17 //--------------------------------------------------------------//
  18 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity /* =MED_CELL */) :
  19 //--------------------------------------------------------------//
  20                                 _entity(Entity),
  21                                 _typeConnectivity(MED_NODAL),
  22                                 _numberOfTypes(0),
  23                                 _geometricTypes((medGeometryElement*)NULL),
  24                                 _type((CELLMODEL*)NULL),
  25                                 _entityDimension(0),
  26                                 _count((int*)NULL),
  27                                 _nodal((MEDSKYLINEARRAY*)NULL),
  28                                 _descending((MEDSKYLINEARRAY*)NULL),
  29                                 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
  30                                 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
  31                                 _neighbourhood((MEDSKYLINEARRAY*)NULL),
  32                                 _constituent((CONNECTIVITY*)NULL)
  33 {
  34    MESSAGE("CONNECTIVITY(medEntityMesh Entity=MED_CELL)");
  35 }
  36
  37 /*!
  38    Constructor. /n
  39    Default for Entity is MED_CELL */
  40 //------------------------------------------------------------------------------//
  41 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity /* =MED_CELL */):
  42 //------------------------------------------------------------------------------//
  43                                 _entity(Entity),
  44                                 _typeConnectivity(MED_NODAL),
  45                                 _numberOfTypes(numberOfTypes),
  46                                 _entityDimension(0),
  47                                 _nodal((MEDSKYLINEARRAY*)NULL),
  48                                 _descending((MEDSKYLINEARRAY*)NULL),
  49                                 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
  50                                 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
  51                                 _neighbourhood((MEDSKYLINEARRAY*)NULL),
  52                                 _constituent((CONNECTIVITY*)NULL)
  53 {
  54   MESSAGE("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)");
  55   _geometricTypes = new medGeometryElement[numberOfTypes];
  56   _type = new CELLMODEL[numberOfTypes];
  57   _count = new int[numberOfTypes+1];
  58 }
  59
  60 /*!
  61   Copy Constructor.
  62 */
  63 //------------------------------------------------------//
  64 CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
  65 //----------------------------------------------------//
  66                                 _entity           (m._entity),
  67                                 _typeConnectivity (m._typeConnectivity),
  68                                 _numberOfTypes    (m._numberOfTypes),
  69                                 _entityDimension  (m._entityDimension),
  70                                 _numberOfNodes    (m._numberOfNodes)
  71 {
  72  if (m._geometricTypes != NULL)
  73  {
  74     _geometricTypes = new medGeometryElement[_numberOfTypes];
  75     memcpy(_geometricTypes,m._geometricTypes,_numberOfTypes*sizeof(medGeometryElement));
  76  }
  77  else
  78     _geometricTypes = (medGeometryElement *) NULL;
  79
  80  if (m._type != NULL)
  81  {
  82     _type = new CELLMODEL[_numberOfTypes];
  83     for (int i=0;i<_numberOfTypes;i++)
  84         _type[i] = CELLMODEL(m._type[i]);
  85  }
  86  else
  87     _type = (CELLMODEL *) NULL;
  88
  89  if (m._count != NULL)
  90  {
  91       _count = new med_int[_numberOfTypes+1];
  92       memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(med_int));
  93  }
  94  else
  95     _count = (med_int *) NULL;
  96
  97  if (m._nodal != NULL)
  98     _nodal = new MEDSKYLINEARRAY(* m._nodal);
  99  else
 100     _nodal = (MEDSKYLINEARRAY *) NULL;
 101
 102  if (m._descending != NULL)
 103     _descending = new MEDSKYLINEARRAY(* m._descending);
 104  else
 105     _descending = (MEDSKYLINEARRAY *) NULL;
 106
 107  if (m._reverseNodalConnectivity != NULL)
 108     _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
 109  else
 110     _reverseNodalConnectivity = (MEDSKYLINEARRAY *) NULL;
 111
 112  if (m._reverseDescendingConnectivity != NULL)
 113     _reverseDescendingConnectivity = new MEDSKYLINEARRAY(* m._reverseDescendingConnectivity);
 114  else
 115     _reverseDescendingConnectivity = (MEDSKYLINEARRAY *) NULL;
 116
 117  if (m._neighbourhood != NULL)
 118     _neighbourhood = new MEDSKYLINEARRAY(* m._neighbourhood);
 119  else
 120     _neighbourhood = (MEDSKYLINEARRAY *) NULL;
 121
 122  if (m._constituent != NULL)
 123     _constituent = new CONNECTIVITY(* m._constituent);
 124  else
 125     _constituent = (CONNECTIVITY *) NULL;
 126 }
 127
 128 /*!
 129    Destructor./n
 130    desallocates existing pointers */
 131 //----------------------------//
 132 CONNECTIVITY::~CONNECTIVITY()
 133 //----------------------------//
 134 {
 135   MESSAGE("Destructeur de CONNECTIVITY()");
 136
 137   if (_geometricTypes != NULL)
 138      delete [] _geometricTypes;
 139   if (_type != NULL)
 140      delete [] _type;
 141   if (_count != NULL)
 142      delete [] _count;
 143   if (_nodal != NULL)
 144      delete _nodal;
 145   if (_descending != NULL)
 146      delete _descending;
 147   if (_reverseNodalConnectivity != NULL)
 148      delete _reverseNodalConnectivity;
 149   if (_reverseDescendingConnectivity != NULL)
 150      delete _reverseDescendingConnectivity;
 151   if (_constituent != NULL)
 152      delete _constituent;
 153 }
 154
 155 /*!
 156    set _constituent to Constituent
 157    be aware desallocation of _constituent is done by CONNECTIVITY:~CONNECTIVITY
 158    throws an exception if Constituent = MED_CELL
 159    A DOCUMENTER
 160     */
 161 //----------------------------------------------------------//
 162 void CONNECTIVITY::setConstituent(CONNECTIVITY * Constituent)
 163                     throw (MEDEXCEPTION)
 164 //----------------------------------------------------------//
 165 {
 166   medEntityMesh Entity = Constituent->getEntity();
 167   if (Entity == MED_CELL)
 168     throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : could not set constituent on MED_CELL !"));
 169
 170   if ((Entity == MED_EDGE)&(_entityDimension == 3))
 171   {
 172     if (_constituent == NULL)
 173       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : Entity not found !"));
 174     _constituent->setConstituent(Constituent);
 175   }
 176   else
 177     _constituent = Constituent;
 178 }
 179
 180 /*! Duplicated Types array in CONNECTIVITY object. */
 181 //--------------------------------------------------------------------//
 182 void CONNECTIVITY::setGeometricTypes(const medGeometryElement * Types,
 183                                      const medEntityMesh Entity)
 184                                      throw (MEDEXCEPTION)
 185 //--------------------------------------------------------------------//
 186 {
 187   if (Entity == _entity)
 188     for (int i=0; i<_numberOfTypes; i++)
 189     {
 190       _geometricTypes[i] = Types[i];
 191       _type[i] = CELLMODEL(_geometricTypes[i]);
 192     }
 193   else
 194     {
 195     if (_constituent == NULL)
 196         throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setGeometricTypes : Entity not found !"));
 197     _constituent->setGeometricTypes(Types,Entity);
 198   }
 199 }
 200
 201 /*! A DOCUMENTER */
 202 //--------------------------------------------------------------------//
 203 void CONNECTIVITY::setCount(const int * Count, const medEntityMesh Entity)
 204                             throw (MEDEXCEPTION)
 205 //--------------------------------------------------------------------//
 206 {
 207   if (Entity == _entity)
 208   {
 209     int * index = new int[Count[_numberOfTypes]];
 210     index[0]=1;
 211     _count[0]=1;
 212     for (int i=0; i<_numberOfTypes; i++) {
 213       _count[i+1] = Count[i+1];
 214       int NumberOfNodesPerElement = _type[i].getNumberOfNodes() ;
 215       for (int j=_count[i]; j<_count[i+1]; j++)
 216         index[j] = index[j-1]+NumberOfNodesPerElement;
 217     }
 218     // allocate _nodal
 219     if (_nodal != NULL) delete _nodal;
 220     _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,index[_count[_numberOfTypes]-1]-1);
 221     _nodal->setIndex(index);
 222     delete [] index;
 223   }
 224   else
 225   {
 226     if (_constituent == NULL)
 227       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setCount : Entity not found !"));
 228     _constituent->setCount(Count,Entity);
 229   }
 230 }
 231
 232 //--------------------------------------------------------------------//
 233 void CONNECTIVITY::setNodal(const int * Connectivity,
 234                             const medEntityMesh Entity,
 235                             const medGeometryElement Type)
 236                             throw (MEDEXCEPTION)
 237 //--------------------------------------------------------------------//
 238 {
 239   if (_entity == Entity)
 240   {
 241     // find geometric type
 242     bool find = false;
 243     for (int i=0; i<_numberOfTypes; i++)
 244       if (_geometricTypes[i] == Type) {
 245         find = true;
 246         int NumberOfNodePerElements = _type[i].getNumberOfNodes() ;
 247         //_nodal->setI(i+1,Connectivity);
 248         for( int j=_count[i];j<_count[i+1]; j++)
 249           _nodal->setI(j,Connectivity+(j-_count[i])*NumberOfNodePerElements);
 250       }
 251     if (!find)
 252       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : geometric type not found"));
 253   } else
 254   {
 255     if (_constituent == NULL)
 256       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : Entity not found !"));
 257       _constituent->setNodal(Connectivity,Entity,Type);
 258   }
 259 }
 260
 261 /*! A DOCUMENTER */
 262 //------------------------------------------------------------------------------------------//
 263 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
 264 //------------------------------------------------------------------------------------------//
 265 {
 266   MESSAGE("CONNECTIVITY::calculateConnectivity");
 267
 268   // a temporary limitation due to calculteDescendingConnectivity function !!!
 269   if ((_entityDimension==3) & (Entity==MED_EDGE))
 270     throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
 271
 272   if (Entity==_entity)
 273     if (ConnectivityType==MED_NODAL)
 274       calculateNodalConnectivity();
 275     else
 276       if (Entity==MED_CELL)
 277         calculateDescendingConnectivity();
 278       else
 279         throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
 280   if (Entity!=_entity) {
 281     calculateDescendingConnectivity();
 282     if (_entityDimension == 2 || _entityDimension == 3)
 283       _constituent->calculateConnectivity(ConnectivityType,Entity);
 284   }
 285 }
 286
 287 /*!  Give, in full or no interlace mode (for nodal connectivity),
 288      descending or nodal connectivity.
 289
 290       You must give a <medEntityMesh> (ie:MED_EDGE)
 291       and a <medGeometryElement> (ie:MED_SEG3).  */
 292
 293 //------------------------------------------------------------//
 294 void CONNECTIVITY::updateFamily(vector<FAMILY*> myFamilies)
 295 //------------------------------------------------------------//
 296 {
 297   const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
 298   BEGIN_OF(LOC);
 299
 300   int numberOfFamilies = myFamilies.size();
 301   if (numberOfFamilies == 0 ) {
 302     MESSAGE(LOC<<"No family");
 303     return;
 304   }
 305   // does we do an update ?
 306   if ((_constituent != NULL)&(_descending != NULL)) {
 307     MESSAGE(LOC<<"Constituent is already defined");
 308     return;
 309   }
 310
 311   if ((_constituent != NULL)&(_descending == NULL)) {
 312     if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
 313       MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
 314       return;
 315     }
 316
 317     // well we could go !
 318     CONNECTIVITY * oldConstituent = _constituent;
 319
 320 //      for(int i=0; i<numberOfFamilies; i++) {
 321 //        FAMILY * myFamily = myFamilies[i];
 322 //        MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
 323 //      }
 324
 325     _constituent = (CONNECTIVITY *)NULL;
 326     // for instance we must have nodal connectivity in constituent :
 327     if (oldConstituent->_nodal == NULL)
 328       {
 329         MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
 330         throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
 331       }
 332     int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
 333     const int * oldConstituentValue = oldConstituent->_nodal->getValue();
 334     const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
 335
 336     calculateDescendingConnectivity();
 337
 338     //    int newNumberOfFace = _constituent->_nodal->getNumberOf();
 339     const int * newConstituentValue = _constituent->_nodal->getValue();
 340     const int * newConstituentIndex = _constituent->_nodal->getIndex();
 341
 342     const int * newReverseDescendingIndex =
 343       _reverseDescendingConnectivity->getIndex();
 344
 345     const int * newDescendingIndex = _descending->getIndex();
 346     //    const int * newDescendingValue = _descending->getValue();
 347
 348     // loop on all family,
 349     //   for all constituent in family, we get it's old connectivity
 350     //   (with oldCconstituentValue and oldConstituentIndex)
 351     //   and search the constituent in newConstituentValue with class
 352     //   ModulusArry
 353     //
 354     //   when a new face is found, replace old constituent
 355     //   number in family with new one
 356     //   If face is not rigth oriented, we must change _descending attribute
 357     //   and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
 358
 359     // Voila a toi de jouer Nadir :-)
 360
 361     // First we get the renumbering from the oldCconstituentValue and
 362     // oldConstituentIndex in the the new one, newConstituentValue and
 363     // newConstituentIndex with a negative sign if the face is not
 364     // right orented
 365
 366     int * renumberingFromOldToNew = new int [oldNumberOfFace];
 367     int index1 = 0;
 368     int indexm1 = 0;
 369
 370     _constituent->calculateReverseNodalConnectivity();
 371
 372     for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
 373       {
 374         int index = 0;
 375
 376         renumberingFromOldToNew[iOldFace] = iOldFace+1;
 377         //        renumberingFromOldToNew[iOldFace] = 999999;
 378
 379         int face_it_beginOld = oldConstituentIndex[iOldFace];
 380         int face_it_endOld = oldConstituentIndex[iOldFace+1];
 381         int face_size_itOld = face_it_endOld - face_it_beginOld;
 382
 383         const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
 384         int face_size_itNew;
 385
 386         const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
 387         const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
 388
 389         // set an array wich contains faces numbers arround first node
 390         int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
 391         int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
 392         int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
 393
 394         int * FacesList = new int[NumberOfFacesInList];
 395
 396         for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
 397           FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
 398         }
 399         // foreach node in sub cell, we search elements which are in common
 400         // at the end, we must have only one !
 401
 402         for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
 403           {
 404             int NewNumberOfFacesInList = 0;
 405             int * NewFacesList = new int[NumberOfFacesInList];
 406
 407             for (int l1=0; l1<NumberOfFacesInList; l1++) {
 408               for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
 409                 if (FacesList[l1]<reverseFaceNodal[l2-1])
 410                   // increasing order : FacesList[l1] are not in elements list of node l
 411                   break;
 412                 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
 413                   // we have found one
 414                   NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
 415                   NewNumberOfFacesInList++;
 416                   break;
 417                 }
 418               }
 419             }
 420             NumberOfFacesInList = NewNumberOfFacesInList;
 421             delete [] FacesList;
 422             FacesList = NewFacesList;
 423           }
 424
 425         if (!NumberOfFacesInList==0)
 426           {
 427             if (NumberOfFacesInList>1)
 428               throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
 429
 430             MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
 431
 432             int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
 433             int face_it_endNew = newConstituentIndex[FacesList[0]];
 434             face_size_itNew = face_it_endNew - face_it_beginNew;
 435
 436             const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
 437             MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
 438
 439             int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
 440
 441             //SCRUTE(retCompareNewOld);
 442
 443             // Real new face found
 444
 445             if(retCompareNewOld == 1)
 446               {
 447                 renumberingFromOldToNew[iOldFace] = FacesList[0];
 448                 index = 1;
 449                 index1++;
 450               }
 451
 452             // Reverse new face found
 453
 454             if(retCompareNewOld == -1)
 455               {
 456                 renumberingFromOldToNew[iOldFace] = FacesList[0];
 457                 index = 1;
 458                 indexm1++;
 459
 460                 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
 461                 int face_it_end = newReverseDescendingIndex[FacesList[0]];
 462                 int face_size_it = face_it_end - face_it_begin;
 463
 464                 if (face_size_it == 1)
 465                   throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
 466
 467                 if (face_size_it > 2)
 468                   throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
 469
 470                 // we have always 2 neighbourings
 471                 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
 472                 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
 473                 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
 474                 //                  if (cell2 == 0)
 475                 //                    throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
 476
 477                 if (cell2 != 0) { // we are not on border !!!!
 478
 479                   _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
 480                   // Updating _constituent->_nodal because of reversity
 481                   const int * oldArray = oldConstituentValue+face_it_beginOld-1;
 482                   for(int iarray=1;iarray<=face_size_itNew;iarray++){
 483                     _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
 484                   }
 485
 486                   // Updating _reverseDescendingConnectivity
 487
 488
 489                   _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
 490
 491                   // Updating _descending for cell1 and cell2
 492                   for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
 493                     if (_descending->getIndexValue(iface)==FacesList[0])
 494                       _descending->setIndexValue(iface,-FacesList[0]);
 495                     else if (_descending->getIndexValue(iface)==-FacesList[0])
 496                       _descending->setIndexValue(iface,FacesList[0]);
 497                   // else nothing to do
 498
 499                   for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
 500                     if (_descending->getIndexValue(iface)==FacesList[0])
 501                       _descending->setIndexValue(iface,-FacesList[0]);
 502                     else if (_descending->getIndexValue(iface)==-FacesList[0])
 503                       _descending->setIndexValue(iface,FacesList[0]);
 504                   // else nothing to do
 505
 506                 } else {// else we are on border and we do nothing !!!!!!!!
 507                   MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
 508                   MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
 509                   MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
 510                 }
 511               }
 512
 513             if(index == 0)
 514               {
 515                 MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
 516                 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
 517               }
 518           }
 519         delete [] FacesList;
 520       }
 521
 522     MESSAGE(LOC<<"The Renumbering is finished and the status is");
 523
 524     // Updating the Family
 525
 526     for(int i=0; i<numberOfFamilies; i++) {
 527       FAMILY * myFamily = myFamilies[i];
 528
 529       MEDSKYLINEARRAY * number = myFamily->getnumber();
 530       int numberOfLines_skyline = number->getNumberOf();
 531       const int * index_skyline = number->getIndex();
 532
 533       for (int i=0;i<numberOfLines_skyline;i++) {
 534         for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
 535           number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
 536         }
 537       }
 538       MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
 539     }
 540
 541     delete oldConstituent ;
 542     delete [] renumberingFromOldToNew;
 543   }
 544
 545
 546   END_OF(LOC);
 547   return;
 548 }
 549
 550 // meme methode que updateFamily, mais avec des groupes. Il n'est pas possible d'utiliser
 551 // l'heritage car les pointeurs sont dans un conteneur.
 552 void CONNECTIVITY::updateGroup(vector<GROUP*> myFamilies)
 553 //------------------------------------------------------------//
 554 {
 555   const char * LOC = "CONNECTIVITY::updateFamily(vector<GROUP*>) ";
 556   BEGIN_OF(LOC);
 557
 558   int numberOfFamilies = myFamilies.size();
 559   if (numberOfFamilies == 0 ) {
 560     MESSAGE(LOC<<"No family");
 561     return;
 562   }
 563   // does we do an update ?
 564   if ((_constituent != NULL)&(_descending != NULL)) {
 565     MESSAGE(LOC<<"Constituent is already defined");
 566     return;
 567   }
 568
 569   if ((_constituent != NULL)&(_descending == NULL)) {
 570     if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
 571       MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
 572       return;
 573     }
 574
 575     // well we could go !
 576     CONNECTIVITY * oldConstituent = _constituent;
 577
 578 //      for(int i=0; i<numberOfFamilies; i++) {
 579 //        FAMILY * myFamily = myFamilies[i];
 580 //        MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
 581 //      }
 582
 583     _constituent = (CONNECTIVITY *)NULL;
 584     // for instance we must have nodal connectivity in constituent :
 585     if (oldConstituent->_nodal == NULL)
 586       {
 587         MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
 588         throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
 589       }
 590     int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
 591     const int * oldConstituentValue = oldConstituent->_nodal->getValue();
 592     const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
 593
 594     calculateDescendingConnectivity();
 595
 596     //    int newNumberOfFace = _constituent->_nodal->getNumberOf();
 597     const int * newConstituentValue = _constituent->_nodal->getValue();
 598     const int * newConstituentIndex = _constituent->_nodal->getIndex();
 599
 600     const int * newReverseDescendingIndex =
 601       _reverseDescendingConnectivity->getIndex();
 602
 603     const int * newDescendingIndex = _descending->getIndex();
 604     //    const int * newDescendingValue = _descending->getValue();
 605
 606     // loop on all family,
 607     //   for all constituent in family, we get it's old connectivity
 608     //   (with oldCconstituentValue and oldConstituentIndex)
 609     //   and search the constituent in newConstituentValue with class
 610     //   ModulusArry
 611     //
 612     //   when a new face is found, replace old constituent
 613     //   number in family with new one
 614     //   If face is not rigth oriented, we must change _descending attribute
 615     //   and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
 616
 617     // Voila a toi de jouer Nadir :-)
 618
 619     // First we get the renumbering from the oldCconstituentValue and
 620     // oldConstituentIndex in the the new one, newConstituentValue and
 621     // newConstituentIndex with a negative sign if the face is not
 622     // right orented
 623
 624     int * renumberingFromOldToNew = new int [oldNumberOfFace];
 625     int index1 = 0;
 626     int indexm1 = 0;
 627
 628     _constituent->calculateReverseNodalConnectivity();
 629
 630     for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
 631       {
 632         int index = 0;
 633
 634         renumberingFromOldToNew[iOldFace] = iOldFace+1;
 635         //        renumberingFromOldToNew[iOldFace] = 999999;
 636
 637         int face_it_beginOld = oldConstituentIndex[iOldFace];
 638         int face_it_endOld = oldConstituentIndex[iOldFace+1];
 639         int face_size_itOld = face_it_endOld - face_it_beginOld;
 640
 641         const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
 642         int face_size_itNew;
 643
 644         const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
 645         const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
 646
 647         // set an array wich contains faces numbers arround first node
 648         int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
 649         int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
 650         int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
 651
 652         int * FacesList = new int[NumberOfFacesInList];
 653
 654         for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
 655           FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
 656         }
 657         // foreach node in sub cell, we search elements which are in common
 658         // at the end, we must have only one !
 659
 660         for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
 661           {
 662             int NewNumberOfFacesInList = 0;
 663             int * NewFacesList = new int[NumberOfFacesInList];
 664
 665             for (int l1=0; l1<NumberOfFacesInList; l1++) {
 666               for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
 667                 if (FacesList[l1]<reverseFaceNodal[l2-1])
 668                   // increasing order : FacesList[l1] are not in elements list of node l
 669                   break;
 670                 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
 671                   // we have found one
 672                   NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
 673                   NewNumberOfFacesInList++;
 674                   break;
 675                 }
 676               }
 677             }
 678             NumberOfFacesInList = NewNumberOfFacesInList;
 679             delete [] FacesList;
 680             FacesList = NewFacesList;
 681           }
 682
 683         if (!NumberOfFacesInList==0)
 684           {
 685             if (NumberOfFacesInList>1)
 686               throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
 687
 688             MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
 689
 690             int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
 691             int face_it_endNew = newConstituentIndex[FacesList[0]];
 692             face_size_itNew = face_it_endNew - face_it_beginNew;
 693
 694             const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
 695             MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
 696
 697             int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
 698
 699             //SCRUTE(retCompareNewOld);
 700
 701             // Real new face found
 702
 703             if(retCompareNewOld == 1)
 704               {
 705                 renumberingFromOldToNew[iOldFace] = FacesList[0];
 706                 index = 1;
 707                 index1++;
 708               }
 709
 710             // Reverse new face found
 711
 712             if(retCompareNewOld == -1)
 713               {
 714                 renumberingFromOldToNew[iOldFace] = FacesList[0];
 715                 index = 1;
 716                 indexm1++;
 717
 718                 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
 719                 int face_it_end = newReverseDescendingIndex[FacesList[0]];
 720                 int face_size_it = face_it_end - face_it_begin;
 721
 722                 if (face_size_it == 1)
 723                   throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
 724
 725                 if (face_size_it > 2)
 726                   throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
 727
 728                 // we have always 2 neighbourings
 729                 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
 730                 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
 731                 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
 732                 //                  if (cell2 == 0)
 733                 //                    throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
 734
 735                 if (cell2 != 0) { // we are not on border !!!!
 736
 737                   _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
 738                   // Updating _constituent->_nodal because of reversity
 739                   const int * oldArray = oldConstituentValue+face_it_beginOld-1;
 740                   for(int iarray=1;iarray<=face_size_itNew;iarray++){
 741                     _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
 742                   }
 743
 744                   // Updating _reverseDescendingConnectivity
 745
 746
 747                   _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
 748
 749                   // Updating _descending for cell1 and cell2
 750                   for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
 751                     if (_descending->getIndexValue(iface)==FacesList[0])
 752                       _descending->setIndexValue(iface,-FacesList[0]);
 753                     else if (_descending->getIndexValue(iface)==-FacesList[0])
 754                       _descending->setIndexValue(iface,FacesList[0]);
 755                   // else nothing to do
 756
 757                   for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
 758                     if (_descending->getIndexValue(iface)==FacesList[0])
 759                       _descending->setIndexValue(iface,-FacesList[0]);
 760                     else if (_descending->getIndexValue(iface)==-FacesList[0])
 761                       _descending->setIndexValue(iface,FacesList[0]);
 762                   // else nothing to do
 763
 764                 } else {// else we are on border and we do nothing !!!!!!!!
 765                   MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
 766                   MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
 767                   MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
 768                 }
 769               }
 770
 771             if(index == 0)
 772               {
 773                 MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
 774                 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
 775               }
 776           }
 777         delete [] FacesList;
 778       }
 779
 780     MESSAGE(LOC<<"The Renumbering is finished and the status is");
 781
 782     // Updating the Family
 783
 784     for(int i=0; i<numberOfFamilies; i++) {
 785       GROUP * myFamily = myFamilies[i];
 786
 787       MEDSKYLINEARRAY * number = myFamily->getnumber();
 788       int numberOfLines_skyline = number->getNumberOf();
 789       const int * index_skyline = number->getIndex();
 790
 791       for (int i=0;i<numberOfLines_skyline;i++) {
 792         for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
 793           number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
 794         }
 795       }
 796       MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
 797     }
 798
 799     delete oldConstituent ;
 800     delete [] renumberingFromOldToNew;
 801   }
 802
 803
 804   END_OF(LOC);
 805   return;
 806 }
 807
 808 //------------------------------------------------------------------------------------------------------------------//
 809 const med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
 810 //------------------------------------------------------------------------------------------------------------------//
 811 {
 812   const char * LOC = "CONNECTIVITY::getConnectivity";
 813   BEGIN_OF(LOC);
 814
 815   MEDSKYLINEARRAY * Connectivity;
 816   if (Entity==_entity) {
 817
 818     if (ConnectivityType==MED_NODAL)
 819       {
 820         calculateNodalConnectivity();
 821         Connectivity=_nodal;
 822       }
 823     else
 824       {
 825         calculateDescendingConnectivity();
 826         Connectivity=_descending;
 827       }
 828
 829     if (Connectivity!=NULL)
 830       if (Type==MED_ALL_ELEMENTS)
 831         return Connectivity->getValue();
 832       else {
 833         for (med_int i=0; i<_numberOfTypes; i++)
 834           if (_geometricTypes[i]==Type)
 835             //return Connectivity->getI(i+1);
 836             return Connectivity->getI(_count[i]);
 837         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
 838       }
 839     else
 840       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 841   } else
 842     if (_constituent != NULL)
 843       return _constituent->getConnectivity(ConnectivityType,Entity,Type);
 844
 845   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 846 }
 847
 848 /*!  Give morse index array to use with
 849      getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
 850
 851       Each value give start index for corresponding entity in connectivity array.
 852
 853       Example : i-th element, j-th node of it :
 854       - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
 855       - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
 856 //-----------------------------------------------------------------------------------------------//
 857 const med_int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
 858 //----0000000--------------------------------------------------------------------------------------------//
 859 {
 860   const char * LOC = "CONNECTIVITY::getConnectivityIndex";
 861   BEGIN_OF(LOC);
 862
 863   MEDSKYLINEARRAY * Connectivity;
 864   if (Entity==_entity) {
 865
 866     if (ConnectivityType==MED_NODAL)
 867       Connectivity=_nodal;
 868     else
 869       Connectivity=_descending;
 870
 871     if (Connectivity!=NULL)
 872       return Connectivity->getIndex();
 873     else
 874       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 875
 876   } else
 877     if (_constituent != NULL)
 878       return _constituent->getConnectivityIndex(ConnectivityType,Entity);
 879
 880   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 881 }
 882
 883 /*! A DOCUMENTER */
 884 //--------------------------------------------------------------//
 885 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
 886 //--------------------------------------------------------------//
 887 {
 888   const char * LOC = "CONNECTIVITY::getType";
 889   BEGIN_OF(LOC);
 890
 891   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 892     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !"));
 893   for (med_int i=0; i<_numberOfTypes; i++)
 894     if (_geometricTypes[i]==Type)
 895       return _type[i];
 896   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 897 }
 898
 899 /*!  Returns the number of elements of type <med geometrie element>.
 900      Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
 901 //------------------------------------------------------------------------//
 902 med_int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
 903 //------------------------------------------------------------------------//
 904 {
 905   const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
 906   BEGIN_OF(LOC);
 907
 908   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 909     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
 910   for (med_int i=0; i<_numberOfTypes; i++)
 911     if (_geometricTypes[i]==Type)
 912       return _type[i].getNumberOfNodes();
 913   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 914 }
 915
 916 /*!  Returns the number of geometric sub cells of <med geometrie element> type.
 917 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
 918 //------------------------------------------------------------------------//
 919 med_int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
 920 //------------------------------------------------------------------------//
 921 {
 922   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 923     throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
 924   for (med_int i=0; i<_numberOfTypes; i++)
 925     if (_geometricTypes[i]==Type)
 926       return _type[i].getNumberOfConstituents(1);
 927   throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
 928 }
 929
 930 /*!  Returns the number of elements of type <med geometrie element>.
 931
 932      Note :
 933       - Implemented for MED_ALL_ELEMENTS
 934       - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
 935       - Not implemented for MED_NONE */
 936 //-----------------------------------------------------------------------------------//
 937 med_int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
 938 //-----------------------------------------------------------------------------------//
 939 {
 940   const char * LOC = "CONNECTIVITY::getNumberOf";
 941   BEGIN_OF(LOC);
 942
 943   MESSAGE(LOC<<" Entity = "<< Entity << ", _entity = "<<_entity);
 944
 945   if (Entity==_entity) {
 946     if (Type==MED_NONE)
 947       return 0; // not defined !
 948     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
 949     if (Type==MED_ALL_ELEMENTS)
 950       return _count[_numberOfTypes]-1;
 951     for (med_int i=0; i<_numberOfTypes; i++)
 952       if (_geometricTypes[i]==Type)
 953         return (_count[i+1] - _count[i]);
 954     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 955   } else
 956     if (_constituent != NULL)
 957       return _constituent->getNumberOf(Entity,Type);
 958
 959   return 0; // valid if they are nothing else !
 960   //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 961 }
 962
 963 /*! A DOCUMENTER */
 964 //--------------------------------------------------------------//
 965 const med_int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
 966                                 medGeometryElement Type)
 967 //--------------------------------------------------------------//
 968 {
 969   if (TypeConnectivity == MED_NODAL)
 970     {
 971           calculateNodalConnectivity();
 972           if (Type==MED_ALL_ELEMENTS)
 973             return _nodal->getValue();
 974           for (med_int i=0; i<_numberOfTypes; i++)
 975             if (_geometricTypes[i]==Type)
 976               return _nodal->getI(_count[i]);
 977     }
 978   else
 979     {
 980       calculateDescendingConnectivity();
 981       if (Type==MED_ALL_ELEMENTS)
 982         return _descending->getValue();
 983       for (med_int i=0; i<_numberOfTypes; i++)
 984         if (_geometricTypes[i]==Type)
 985           return _descending->getI(Type);
 986     }
 987   throw MEDEXCEPTION("Not found");
 988 }
 989
 990 /*! A DOCUMENTER */
 991 //---------------------------------------------------------------------//
 992 const med_int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
 993 //---------------------------------------------------------------------//
 994 {
 995   if (TypeConnectivity == MED_NODAL)
 996     {
 997       calculateNodalConnectivity();
 998       return _nodal->getIndex();
 999     }
1000   else
1001     {
1002       calculateDescendingConnectivity();
1003       return _descending->getIndex();
1004     }
1005 }
1006
1007 /*! Not yet implemented */
1008 //----------------------------------------------//
1009 const med_int* CONNECTIVITY:: getNeighbourhood() const
1010 //----------------------------------------------//
1011 {
1012   throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
1013 }
1014
1015 /*! Returns an array which contains, for each node, all cells
1016     arround it. */
1017 //-------------------------------------------------//
1018 const med_int* CONNECTIVITY::getReverseNodalConnectivity()
1019 //-------------------------------------------------//
1020 {
1021   calculateReverseNodalConnectivity();
1022   return _reverseNodalConnectivity->getValue();
1023 }
1024
1025 /*!  Give index array to use with getReverseConnectivity(MED_NODAL).
1026      It is unusefull with MED_DESCENDING mode, because we have allways two cells.  */
1027 //-------------------------------------------------------//
1028 const med_int* CONNECTIVITY::getReverseNodalConnectivityIndex()
1029 //-------------------------------------------------------//
1030 {
1031   calculateReverseNodalConnectivity();
1032   return _reverseNodalConnectivity->getIndex();
1033 }
1034
1035 /*! Returns an array which contains, for each face (or edge),
1036     the 2 cells of each side. First is cell which face normal is outgoing.
1037     arround it. */
1038 //------------------------------------------------------//
1039 const med_int* CONNECTIVITY::getReverseDescendingConnectivity()
1040 //------------------------------------------------------//
1041 {
1042   // it is in _constituent connectivity only if we are in MED_CELL
1043   // (we could not for instance calculate face-edge connectivity !)
1044   if (_entity!=MED_CELL)
1045     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
1046
1047   // we need descending connectivity
1048   calculateDescendingConnectivity();
1049   return _reverseDescendingConnectivity->getValue();
1050 }
1051
1052 /*! calculate the reverse descending Connectivity
1053     and returns the index  ( A DOCUMENTER MIEUX)*/
1054 //-----------------------------------------------------------//
1055 const med_int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
1056 //-----------------------------------------------------------//
1057 {
1058   // it is in _constituent connectivity only if we are in MED_CELL
1059   if (_entity!=MED_CELL)
1060     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
1061
1062   // we need descending connectivity
1063   calculateDescendingConnectivity();
1064   return _reverseDescendingConnectivity->getIndex();
1065 }
1066
1067 /*! A DOCUMENTER (et a finir ???) */
1068 //--------------------------------------------//
1069 void CONNECTIVITY::calculateNodalConnectivity()
1070 //--------------------------------------------//
1071 {
1072   if (_nodal==NULL)
1073     {
1074       if (_descending==NULL)
1075         throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
1076       // calculate _nodal from _descending
1077     }
1078 }
1079
1080 /*! If not yet done, calculate the nodal Connectivity
1081     and the reverse nodal Connectivity*/
1082 //---------------------------------------------------//
1083 void CONNECTIVITY::calculateReverseNodalConnectivity()
1084 //---------------------------------------------------//
1085 {
1086   const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
1087   BEGIN_OF(LOC);
1088
1089   if (_nodal==NULL)
1090     calculateNodalConnectivity();
1091
1092   if(_reverseNodalConnectivity==NULL) {
1093
1094     med_int node_number = 0;
1095     vector <vector <med_int> > reverse_connectivity;
1096     reverse_connectivity.resize(_numberOfNodes+1);
1097
1098   // Treat all cells types
1099
1100     for (med_int j = 0; j < _numberOfTypes; j++)
1101       {
1102         // node number of the cell type
1103         node_number = _type[j].getNumberOfNodes();
1104         // treat all cells of a particular type
1105         for (med_int k = _count[j]; k < _count[j+1]; k++)
1106           // treat all nodes of the cell type
1107           for (med_int local_node_number = 1; local_node_number < node_number+1; local_node_number++)
1108             {
1109               med_int global_node = _nodal->getIJ(k,local_node_number);
1110               reverse_connectivity[global_node].push_back(k);
1111             }
1112       }
1113
1114     // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
1115
1116   //calculate size of reverse_nodal_connectivity array
1117     med_int size_reverse_nodal_connectivity  = 0;
1118     for (med_int i = 1; i < _numberOfNodes+1; i++)
1119       size_reverse_nodal_connectivity += reverse_connectivity[i].size();
1120
1121     //MEDSKYLINEARRAY * ReverseConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity);
1122     med_int * reverse_nodal_connectivity_index = new (med_int)[_numberOfNodes+1];
1123     med_int * reverse_nodal_connectivity = new (med_int)[size_reverse_nodal_connectivity];
1124     //const med_int * reverse_nodal_connectivity = ReverseConnectivity->getValue();
1125     //const med_int * reverse_nodal_connectivity_index = ReverseConnectivity->getIndex();
1126
1127     reverse_nodal_connectivity_index[0] = 1;
1128     for (med_int i = 1; i < _numberOfNodes+1; i++)
1129       {
1130         med_int size = reverse_connectivity[i].size();
1131         reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
1132         for (med_int j = 0; j < size; j++)
1133           reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
1134       }
1135
1136     //_reverseNodalConnectivity = ReverseConnectivity;
1137     _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
1138                                                     reverse_nodal_connectivity_index,
1139                                                     reverse_nodal_connectivity);
1140     delete [] reverse_nodal_connectivity_index;
1141     delete [] reverse_nodal_connectivity;
1142  }
1143 }
1144
1145 /*! If not yet done, calculate the Descending Connectivity */
1146 //-------------------------------------------------//
1147 void CONNECTIVITY::calculateDescendingConnectivity()
1148 //-------------------------------------------------//
1149 {
1150   const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
1151   BEGIN_OF(LOC);
1152
1153   if (_descending==NULL)
1154     {
1155       if (_nodal==NULL) {
1156         MESSAGE(LOC<<"No connectivity found !");
1157         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
1158       }
1159       // calcul _descending from _nodal
1160       // we need CONNECTIVITY for constituent
1161
1162       if (_constituent != NULL)
1163         //      throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
1164         MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
1165
1166       if (_entityDimension == 3)
1167         _constituent = new CONNECTIVITY(MED_FACE);
1168       else if (_entityDimension == 2)
1169         _constituent = new CONNECTIVITY(MED_EDGE);
1170       else {
1171         MESSAGE(LOC<<"We are in 1D");
1172         return;
1173       }
1174
1175       _constituent->_typeConnectivity = MED_NODAL;
1176       _constituent->_numberOfNodes = _numberOfNodes;
1177       // foreach cells, we built array of constituent
1178       int DescendingSize = 0;
1179       for(int i=0; i<_numberOfTypes; i++)
1180         DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
1181       //_descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize);
1182       //const int * descend_connectivity = _descending->getValue();
1183       int * descend_connectivity = new int[DescendingSize];
1184       for (int i=0; i<DescendingSize; i++)
1185         descend_connectivity[i]=0;
1186       //const int * descend_connectivity_index = _descending->getIndex();
1187       int * descend_connectivity_index = new int[_count[_numberOfTypes]];
1188       descend_connectivity_index[0]=1;
1189       medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
1190       ConstituentsTypes[0]=MED_NONE;
1191       ConstituentsTypes[1]=MED_NONE;
1192       int * NumberOfConstituentsForeachType = new int[2];
1193       NumberOfConstituentsForeachType[0]=0;
1194       NumberOfConstituentsForeachType[1]=0;
1195       for(int i=0; i<_numberOfTypes; i++) {
1196         // initialize descend_connectivity_index array :
1197         int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
1198         for (int j=_count[i];j<_count[i+1];j++){
1199           descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
1200           // compute number of constituent of all cell for each type
1201           for(int k=1;k<NumberOfConstituents+1;k++) {
1202             medGeometryElement MEDType = _type[i].getConstituentType(1,k);
1203             if (ConstituentsTypes[0]==MED_NONE) {
1204               ConstituentsTypes[0]=MEDType;
1205               NumberOfConstituentsForeachType[0]++;
1206             } else if (ConstituentsTypes[0]==MEDType)
1207               NumberOfConstituentsForeachType[0]++;
1208             else if (ConstituentsTypes[1]==MED_NONE) {
1209               ConstituentsTypes[1]=MEDType;
1210               NumberOfConstituentsForeachType[1]++;
1211             } else if (ConstituentsTypes[1]==MEDType)
1212               NumberOfConstituentsForeachType[1]++;
1213             else
1214               throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
1215           }
1216         }
1217       }
1218
1219       // we could built _constituent !
1220       int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
1221       int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
1222
1223       //_constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes);
1224
1225       // we use _constituent->_nodal
1226       //const int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
1227       //const int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
1228       int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
1229       int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
1230       ConstituentNodalConnectivityIndex[0]=1;
1231
1232       _constituent->_entityDimension=ConstituentsTypes[0]/100;
1233       if (ConstituentsTypes[1]==MED_NONE)
1234         _constituent->_numberOfTypes = 1;
1235       else
1236         _constituent->_numberOfTypes = 2;
1237       _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
1238       _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes];
1239       _constituent->_count = new med_int[_constituent->_numberOfTypes+1];
1240       _constituent->_count[0]=1;
1241       int* tmp_NumberOfConstituentsForeachType = new med_int[_constituent->_numberOfTypes+1];
1242       tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
1243       for (int i=0; i<_constituent->_numberOfTypes;i++) {
1244         _constituent->_geometricTypes[i]=ConstituentsTypes[i];
1245         _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
1246         CELLMODEL Type(ConstituentsTypes[i]);
1247         _constituent->_type[i]=Type;
1248         tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
1249         for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
1250           ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
1251       }
1252       delete [] ConstituentsTypes;
1253       delete [] NumberOfConstituentsForeachType;
1254
1255       // we need reverse nodal connectivity
1256       if (! _reverseNodalConnectivity)
1257         calculateReverseNodalConnectivity();
1258       const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
1259       const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
1260
1261       // array to keep reverse descending connectivity
1262       int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
1263
1264       int TotalNumberOfSubCell = 0;
1265       for (int i=0; i<_numberOfTypes; i++) { // loop on all cell type
1266         CELLMODEL Type = _type[i];
1267         //      int NumberOfNodesPerCell = Type.getNumberOfNodes();
1268         int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
1269         for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
1270           for (int k=1; k<=NumberOfConstituentPerCell; k++)
1271             { // we loop on all sub cell of it
1272               if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
1273                 // it is a new sub cell !
1274                 //          TotalNumberOfSubCell++;
1275                 // Which type ?
1276                 if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
1277                   tmp_NumberOfConstituentsForeachType[0]++;
1278                   TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
1279                 } else {
1280                   tmp_NumberOfConstituentsForeachType[1]++;
1281                   TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
1282                 }
1283               //we have maximum two types
1284
1285                 descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
1286                 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
1287                 int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
1288
1289                 int * NodesLists = new int[NumberOfNodesPerConstituent];
1290                 for (int l=0; l<NumberOfNodesPerConstituent; l++) {
1291                   NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
1292                   ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
1293                 }
1294                 // we use reverse_nodal_connectivity to find the other element which contain this sub cell
1295
1296                 // all elements which contains first node of sub cell :
1297                 int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
1298                 int ReverseNodalConnectivityIndex_1 = ReverseNodalConnectivityIndex[NodesLists[0]];
1299                 int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
1300
1301                 if (NumberOfCellsInList > 0)  { // we could have no element !
1302                   int * CellsList = new int[NumberOfCellsInList];
1303                   for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
1304                     CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
1305
1306                   // foreach node in sub cell, we search elements which are in common
1307                   // at the end, we must have only one !
1308
1309                   for (int l=1; l<NumberOfNodesPerConstituent; l++) {
1310                     int NewNumberOfCellsInList = 0;
1311                     int * NewCellsList = new int[NumberOfCellsInList];
1312                     for (int l1=0; l1<NumberOfCellsInList; l1++)
1313                       for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++) {
1314                         if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
1315                           // increasing order : CellsList[l1] are not in elements list of node l
1316                           break;
1317                         if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
1318                           // we have found one
1319                           NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
1320                           NewNumberOfCellsInList++;
1321                           break;
1322                         }
1323                       }
1324                     NumberOfCellsInList = NewNumberOfCellsInList;
1325
1326                     delete [] CellsList;
1327                     CellsList = NewCellsList;
1328                   }
1329
1330                   if (NumberOfCellsInList > 0) { // We have found some elements !
1331                     int CellNumber = CellsList[0];
1332                     //delete [] CellsList;
1333                     if (NumberOfCellsInList>1)
1334                       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
1335
1336                     if (NumberOfCellsInList==1) {
1337                       ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
1338
1339                       // we search sub cell number in this cell to not calculate it another time
1340                       // which type ?
1341                       CELLMODEL Type2;
1342                       for (int l=0; l<_numberOfTypes; l++)
1343                         if (CellNumber < _count[l+1]) {
1344                           Type2=_type[l];
1345                           break;
1346                         }
1347                       //int sub_cell_count2 = Type2.get_entities_count(1);
1348                       //int nodes_cell_count2 = Type2.get_nodes_count();
1349                       bool find2 = false;
1350                       for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) { // on all sub cell
1351                         int counter = 0;
1352                         for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
1353                           for (int n=1; n<=Type.getConstituentType(1,k)%100; n++) {
1354                             if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
1355                               counter++;
1356                           }
1357                         if (counter==Type.getConstituentType(1,k)%100) {
1358                           descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
1359                           find2 = true;
1360                         }
1361                         if (find2)
1362                           break;
1363                       }
1364                       if (!find2)
1365                         MESSAGE(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
1366                     }
1367                   } else {
1368                     ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
1369                   }
1370                   delete [] CellsList;
1371                 }
1372
1373                 delete [] NodesLists;
1374               }
1375             }
1376       }
1377       // we adjust _constituent
1378       int NumberOfConstituent=0;
1379       int SizeOfConstituentNodal=0;
1380       for (int i=0;i<_constituent->_numberOfTypes; i++) {
1381         NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
1382         SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
1383       }
1384       // we built new _nodal attribute in _constituent
1385       //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
1386       //const int *ConstituentNodalValue = ConstituentNodal->getValue();
1387       //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
1388       int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
1389       int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
1390       ConstituentNodalIndex[0]=1;
1391       // we build _reverseDescendingConnectivity
1392       //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
1393       //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
1394       //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
1395       int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
1396       int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
1397       reverseDescendingConnectivityIndex[0]=1;
1398
1399       // first constituent type
1400       for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++) {
1401         ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1402         for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1403           ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1404         }
1405         reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1406         for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1407           reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
1408         }
1409       }
1410       // second type if any
1411       if (_constituent->_numberOfTypes==2) {
1412         int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
1413         int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1414         int offset2=offset*2;
1415         int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1416         for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++) {
1417           ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1418           for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1419             ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1420           }
1421           reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1422           for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1423             reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
1424           }
1425         }
1426         _constituent->_count[2]=NumberOfConstituent+1;
1427         // we correct _descending to adjust face number
1428         for(int j=0;j<DescendingSize;j++)
1429           if (descend_connectivity[j]>tmp_NumberOfConstituentsForeachType[0])
1430             descend_connectivity[j]-=offset;
1431       }
1432
1433       delete [] ConstituentNodalConnectivityIndex;
1434       delete [] ConstituentNodalConnectivity;
1435
1436       _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
1437                                         DescendingSize,
1438                                         descend_connectivity_index,
1439                                         descend_connectivity);
1440       delete [] descend_connectivity_index;
1441       delete [] descend_connectivity;
1442       _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,
1443                                                            2*NumberOfConstituent,
1444                                                            reverseDescendingConnectivityIndex,
1445                                                            reverseDescendingConnectivityValue);
1446       delete [] reverseDescendingConnectivityIndex;
1447       delete [] reverseDescendingConnectivityValue;
1448
1449       _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1450       delete [] tmp_NumberOfConstituentsForeachType;
1451
1452       //delete _constituent->_nodal;
1453       _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
1454                                                  SizeOfConstituentNodal,
1455                                                  ConstituentNodalIndex,
1456                                                  ConstituentNodalValue);
1457
1458       delete [] ConstituentNodalIndex;
1459       delete [] ConstituentNodalValue;
1460
1461       delete [] ReverseDescendingConnectivityValue;
1462     }
1463   END_OF(LOC);
1464 }
1465
1466 //-----------------------------------------------------------------------------------------//
1467 //  void CONNECTIVITY::calculateReverseDescendingConnectivity(CONNECTIVITY *myConnectivity)
1468 //-----------------------------------------------------------------------------------------//
1469 //  {
1470 //    if (_entity==MED_CELL)
1471 //      throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1472 //        Connectivity : we could not do that with MED_CELL entity !");
1473 //    if (myConnectivity->getEntity()!=MED_CELL)
1474 //      throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1475 //              Connectivity : we need MED_CELL connectivity !");
1476 //    return;
1477 //  }
1478
1479 /*! Not implemented yet */
1480 //--------------------------------------------------------------------//
1481 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
1482 //--------------------------------------------------------------------//
1483 {
1484   const char * LOC = "CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity) : ";
1485   BEGIN_OF(LOC);
1486
1487   MESSAGE(LOC<<"method not yet implemented " << myConnectivity._entity);
1488   // Mesh dimension !
1489
1490   END_OF(LOC);
1491   return;
1492 }
1493
1494
1495 /*!
1496   Returns the geometry of an element given by its entity type & its global number.
1497
1498   Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
1499 */
1500 //--------------------------------------------------------------------//
1501 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalNumber) const
1502 //--------------------------------------------------------------------//
1503 {
1504   const char * LOC = "medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity, int globalNumber) const : ";
1505   BEGIN_OF(LOC);
1506
1507   // The globalNumber must verify  : 1 <=  globalNumber < _count[_numberOfTypes] (== totalNumberOfElement+1)
1508   if ( (globalNumber < 1) || (globalNumber >  _count[_numberOfTypes]-1 )  )
1509     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |"
1510                                  << 1 <<"| and <= |" << _count[_numberOfTypes]-1 ));
1511
1512   if (_entity==Entity) {
1513     for(int i=1; i<=_numberOfTypes;i++)
1514       if (globalNumber<_count[i])
1515         return _geometricTypes[i-1];
1516   }
1517   else if (_constituent!=NULL)
1518     return _constituent->getElementType(Entity,globalNumber);
1519   else
1520     throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1521   throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
1522
1523   END_OF(LOC);
1524 }
1525
1526
1527 ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
1528 {
1529     os << endl << "------------- Entity = ";
1530     switch (co._entity)
1531     {
1532         case MED_CELL:
1533             os << "MED_CELL";
1534             break;
1535         case MED_FACE:
1536             os << "MED_FACE";
1537             break;
1538         case MED_EDGE:
1539             os << "MED_EDGE";
1540             break;
1541         case MED_NODE:
1542             os << "MED_NODE";
1543             break;
1544         case MED_ALL_ENTITIES:
1545             os << "MED_ALL_ENTITIES";
1546             break;
1547         default:
1548             os << "Unknown";
1549             break;
1550     }
1551     os  << " -------------\n\nMedConnectivity : ";
1552     switch (co._typeConnectivity)
1553     {
1554         case MED_NODAL:
1555             os << "MED_NODAL\n";
1556             break;
1557         case MED_DESCENDING:
1558             os << "MED_DESCENDING\n";
1559             break;
1560         default:
1561             break;
1562     }
1563     os << "Entity dimension : " << co._entityDimension << endl;
1564     os << "Number of nodes : " << co._numberOfNodes << endl;
1565     os << "Number of types : " << co._numberOfTypes << endl;
1566     for (int i=0; i!=co._numberOfTypes ; ++i)
1567         os << "  -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : "
1568            << co._count[i+1]-co._count[i] << " elements" << endl;
1569
1570     if (co._typeConnectivity == MED_NODAL )
1571     {
1572         for (int i=0; i<co._numberOfTypes; i++)
1573         {
1574             os << endl << co._type[i].getName() << " : " << endl;
1575             int numberofelements = co._count[i+1]-co._count[i];
1576             const med_int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, co._geometricTypes[i]);
1577             int numberofnodespercell = co._geometricTypes[i]%100;
1578             for (int j=0;j<numberofelements;j++)
1579             {
1580                 os << setw(6) << j+1 << " : ";
1581                 for (int k=0;k<numberofnodespercell;k++)
1582                     os << connectivity[j*numberofnodespercell+k]<<" ";
1583                 os << endl;
1584             }
1585         }
1586     }
1587     else if (co._typeConnectivity == MED_DESCENDING)
1588     {
1589         int numberofelements = co.getNumberOf( co._entity , MED_ALL_ELEMENTS);
1590         const med_int *connectivity =  co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
1591         const int *connectivity_index =  co.getConnectivityIndex( co._typeConnectivity, co._entity );
1592
1593         for ( int j=0; j!=numberofelements; j++ )
1594         {
1595             os << "Element " << j+1 << " : ";
1596             for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
1597                 os << connectivity[k-1] << " ";
1598             os << endl;
1599         }
1600     }
1601
1602     if (co._constituent)
1603         os << endl << *co._constituent << endl;
1604
1605     return os;
1606 }