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