src/MEDMEM/MEDMEM_Connectivity.cxx

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