src/MEDMEM/MEDMEM_Connectivity.cxx

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